W R, jak znaleźć błąd standardowy średniej?


89

Czy jest jakieś polecenie, aby znaleźć błąd standardowy średniej w R?

Odpowiedzi:



165

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

92

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. sdFunkcja na badania już ma jeden pierwiastek kwadratowy (kod sdjest 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 varradzenie sobie z brakującymi wartościami. varW tym sewywołaniu można użyć wszystkiego, co można przekazać jako nazwany argument .


4
Co ciekawe, twoja funkcja i Ian są prawie identycznie szybkie. Przetestowałem je oba 1000 razy przeciwko 10 ^ 6 milionom rnormalnych drawów (za mało mocy, aby je docisnąć mocniej). I odwrotnie, funkcja plotrix była zawsze wolniejsza niż nawet najwolniejsze przebiegi tych dwóch funkcji - ale pod maską dzieje się też znacznie więcej.
Matt Parker

6
Zauważ, że stderrjest to nazwa funkcji w base.
Tom

3
To bardzo dobra uwaga. Zwykle używam se. Zmieniłem tę odpowiedź, aby to odzwierciedlić.
John

5
Tom, NIE stderrNIE oblicza standardowego błędu, który wyświetladisplay aspects. of connection
prognostyk

15
@forecaster Tom nie powiedział, że stderroblicza błąd standardowy, ostrzegał, że ta nazwa jest używana w bazie, a John pierwotnie nazwał swoją funkcję stderr(sprawdź historię edycji ...).
Molx

60

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))
}

Zauważ, że stderrw basepakiecie istnieje funkcja o nazwie , która robi coś innego, więc może lepiej wybrać inną nazwę dla tej funkcji, np.se
sparrow


3

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łem sum(!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)

wprowadź opis obrazu tutaj


0

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)



0

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:

  • Otrzymujesz natychmiastowe przedziały ufności za pomocą confint()
  • Możesz użyć testów dla różnych hipotez dotyczących średniej, używając na przykład car::linear.hypothesis()
  • Możesz użyć bardziej wyrafinowanych szacunków odchylenia standardowego, jeśli masz pewną heteroskedastyczność, dane klastrowe, dane przestrzenne itp., Patrz pakiet 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)


-12
y <- mean(x, na.rm=TRUE)

sd(y)odchylenie standardowe var(y)wariancji.

Oba wyprowadzenia używają n-1mianownika, więc są oparte na danych przykładowych.

Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.