Jak symulować niestandardową analizę mocy modelu lm (za pomocą R)

W odpowiedzi na ostatnie pytania, które tu mieliśmy .

Miałem nadzieję, że dowiem się, czy ktoś się zetknął lub może udostępnić kod R do wykonania niestandardowej analizy mocy na podstawie symulacji dla modelu liniowego?

Później oczywiście chciałbym rozszerzyć go na bardziej złożone modele, ale wydaje mi się, że dobrze zacząć. Dzięki.

Nie jestem pewien, czy potrzebujesz symulacji dla prostego modelu regresji. Na przykład zobacz papier Portable Power . W przypadku bardziej złożonych modeli, w szczególności efektów mieszanych, pakiet pamm w R wykonuje analizy mocy poprzez symulacje. Zobacz także post Todda Jobe, który ma kod R do symulacji.

Oto kilka źródeł kodu symulacyjnego w R. Nie jestem pewien, czy któryś konkretnie dotyczy modeli liniowych, ale być może stanowią one wystarczający przykład, aby uzyskać sedno:

• Benjamin Bolker napisał świetną książkę dane ekologiczne i modele R . Wczesny szkic całej książki wraz z kodem Sweave jest dostępny online. Rozdział 5 dotyczy analizy i symulacji mocy.

Istnieje jeszcze kilka przykładów symulacji w następujących witrynach:

• http://www.personality-project.org/R/r.datageneration.html
• http://psy-ed.wikidot.com/simulation

Na podstawie modeli ekologicznych i danych Bolker 2009 w R. Musisz zadeklarować siłę trendu (tj. Nachylenie), który chcesz przetestować. Intuicyjnie silny trend i niska zmienność będą wymagały małej wielkości próby, słaby trend i duża zmienność będą wymagały dużej wielkości próby.

a = 2  #desired slope
b = 1  #estimated intercept
sd = 20  #estimated variability defined by standard deviation
nsim = 400  #400 simulations
pval = numeric(nsim)  #placeholder for the second for loop output
Nvec = seq(25, 100, by = 1)  #vector for the range of sample sizes to be tested
power.N = numeric(length(Nvec))   #create placeholder for first for loop output
for (j in 1:length(Nvec)) {
  N = Nvec[j]  
  x = seq(1, 20, length = Nvec[j])  #x value length needs to match sample size (Nvec) length
  for (i in 1:nsim) {   #for this value of N, create random error 400 times
    y_det = a + b * x
    y = rnorm(N, mean = y_det, sd = sd)
    m = lm(y ~ x)
    pval[i] = coef(summary(m))["x", "Pr(>|t|)"]  #all the p values for 400 sims
  }  #cycle through all N values
  power.N[j] = sum(pval < 0.05)/nsim  #the proportion of correct p-values (i.e the power)
}
power.N
plot(Nvec, power.N)  #need about 90 - 100 samples for 80% power

Możesz także zasymulować minimalny trend, który można przetestować dla danej wielkości próbki, jak pokazano w książce

bvec = seq(-2, 2, by = 0.1)
power.b = numeric(length(bvec))
for (j in 1:length(bvec)) {
  b = bvec[j]
   for (i in 1:nsim) {
     y_det = a + b * x
     y = rnorm(N, mean = y_det, sd = sd)
     m = lm(y ~ x)
     pval[i] = coef(summary(m))["x", "Pr(>|t|)"]
     }
   power.b[j] = sum(pval < 0.05)/nsim
  }
 power.b
 plot(bvec, power.b)