Czy jest jakieś polecenie, aby znaleźć błąd standardowy średniej w R?
źródło
Błąd standardowy to po prostu odchylenie standardowe podzielone przez pierwiastek kwadratowy z wielkości próby. Możesz więc łatwo stworzyć własną funkcję:
> std <- function(x) sd(x)/sqrt(length(x))
> std(c(1,2,3,4))
[1] 0.6454972
Błąd standardowy (SE) to po prostu odchylenie standardowe rozkładu próbkowania. Wariancja rozkładu próbkowania to wariancja danych podzielonych przez N, a SE to pierwiastek kwadratowy z tego. Wychodząc z tego zrozumienia, można zauważyć, że bardziej efektywne jest użycie wariancji w obliczeniach SE. sd
Funkcja na badania już ma jeden pierwiastek kwadratowy (kod sd
jest w R i ujawniana przez wpisanie „SD”). Dlatego poniższe są najbardziej wydajne.
se <- function(x) sqrt(var(x)/length(x))
aby funkcja była tylko trochę bardziej złożona i obsługiwała wszystkie opcje, do których można było przejść var
, można wprowadzić tę modyfikację.
se <- function(x, ...) sqrt(var(x, ...)/length(x))
Używając tej składni, można wykorzystać takie rzeczy, jak var
radzenie sobie z brakującymi wartościami. var
W tym se
wywołaniu można użyć wszystkiego, co można przekazać jako nazwany argument .
stderr
jest to nazwa funkcji wbase
.stderr
NIE oblicza standardowego błędu, który wyświetladisplay aspects. of connection
stderr
oblicza błąd standardowy, ostrzegał, że ta nazwa jest używana w bazie, a John pierwotnie nazwał swoją funkcjęstderr
(sprawdź historię edycji ...).Wersja powyższej odpowiedzi Johna, która usuwa brzydkie NA:
stderr <- function(x, na.rm=FALSE) { if (na.rm) x <- na.omit(x) sqrt(var(x)/length(x)) }
źródło
stderr
wbase
pakiecie istnieje funkcja o nazwie , która robi coś innego, więc może lepiej wybrać inną nazwę dla tej funkcji, np.se
Pakiet sciplot ma wbudowaną funkcję se (x)
źródło
Ponieważ od czasu do czasu wracam do tego pytania i ponieważ jest ono stare, zamieszczam punkt odniesienia dla najczęściej głosowanych odpowiedzi.
Zauważ, że dla odpowiedzi @ Iana i @ Johna stworzyłem inną wersję. Zamiast używać
length(x)
, użyłemsum(!is.na(x))
(aby uniknąć NA). Użyłem wektora 10 ^ 6 z 1000 powtórzeń.library(microbenchmark) set.seed(123) myVec <- rnorm(10^6) IanStd <- function(x) sd(x)/sqrt(length(x)) JohnSe <- function(x) sqrt(var(x)/length(x)) IanStdisNA <- function(x) sd(x)/sqrt(sum(!is.na(x))) JohnSeisNA <- function(x) sqrt(var(x)/sum(!is.na(x))) AranStderr <- function(x, na.rm=FALSE) { if (na.rm) x <- na.omit(x) sqrt(var(x)/length(x)) } mbm <- microbenchmark( "plotrix" = {plotrix::std.error(myVec)}, "IanStd" = {IanStd(myVec)}, "JohnSe" = {JohnSe(myVec)}, "IanStdisNA" = {IanStdisNA(myVec)}, "JohnSeisNA" = {JohnSeisNA(myVec)}, "AranStderr" = {AranStderr(myVec)}, times = 1000) mbm
Wyniki:
Unit: milliseconds expr min lq mean median uq max neval cld plotrix 10.3033 10.89360 13.869947 11.36050 15.89165 125.8733 1000 c IanStd 4.3132 4.41730 4.618690 4.47425 4.63185 8.4388 1000 a JohnSe 4.3324 4.41875 4.640725 4.48330 4.64935 9.4435 1000 a IanStdisNA 8.4976 8.99980 11.278352 9.34315 12.62075 120.8937 1000 b JohnSeisNA 8.5138 8.96600 11.127796 9.35725 12.63630 118.4796 1000 b AranStderr 4.3324 4.41995 4.634949 4.47440 4.62620 14.3511 1000 a library(ggplot2) autoplot(mbm)
źródło
bardziej ogólnie, w przypadku standardowych błędów dowolnego innego parametru można użyć pakietu rozruchowego do symulacji ładowania początkowego (lub napisać je samodzielnie)
źródło
Możesz użyć funkcji stat.desc z pakietu pastec.
library(pastec) stat.desc(x, BASIC =TRUE, NORMAL =TRUE)
więcej na ten temat znajdziesz tutaj: https://www.rdocumentation.org/packages/pastecs/versions/1.3.21/topics/stat.desc
źródło
Pamiętając, że średnią można również uzyskać za pomocą modelu liniowego, regresując zmienną względem pojedynczego punktu przecięcia z osią, można również użyć
lm(x~1)
do tego funkcji!Zalety to:
confint()
car::linear.hypothesis()
sandwich
## generate data x <- rnorm(1000) ## estimate reg reg <- lm(x~1) coef(summary(reg))[,"Std. Error"] #> [1] 0.03237811 ## conpare with simple formula all.equal(sd(x)/sqrt(length(x)), coef(summary(reg))[,"Std. Error"]) #> [1] TRUE ## extract confidence interval confint(reg) #> 2.5 % 97.5 % #> (Intercept) -0.06457031 0.0625035
Utworzono 06.10.2020 r. Przez pakiet reprex (v0.3.0)
źródło
y <- mean(x, na.rm=TRUE)
sd(y)
odchylenie standardowevar(y)
wariancji.Oba wyprowadzenia używają
n-1
mianownika, więc są oparte na danych przykładowych.źródło