So simulieren Sie eine benutzerdefinierte Leistungsanalyse eines lm-Modells (mit R)

Im Anschluss an die letzten Fragen hatten wir hier .

Ich wollte wissen, ob jemand auf R-Code gestoßen ist oder diesen für die Durchführung einer benutzerdefinierten Leistungsanalyse auf der Grundlage der Simulation für ein lineares Modell verwenden kann.

Später würde ich es natürlich gerne auf komplexere Modelle ausweiten, aber ich scheine an der richtigen Stelle zu sein, um anzufangen. Vielen Dank.

Ich bin mir nicht sicher, ob Sie für ein einfaches Regressionsmodell eine Simulation benötigen. Ein Beispiel finden Sie im Papier Portable Power . Für komplexere Modelle, insbesondere für gemischte Effekte, führt das Pamm- Paket in R Leistungsanalysen durch Simulationen durch. Siehe auch Todd Jobes Beitrag, der R-Code für die Simulation enthält.

Hier sind ein paar Quellen für Simulationscode in R. Ich bin mir nicht sicher, ob sie sich speziell mit linearen Modellen befassen, aber vielleicht bieten sie genug Beispiele, um das Wesentliche zu verstehen:

• Benjamin Bolker hat mit R ein großartiges Buch über ökologische Daten und Modelle geschrieben . Ein früher Entwurf des gesamten Buches zusammen mit dem Sweave-Code ist online verfügbar. Kapitel 5 befasst sich mit Leistungsanalyse und Simulation.

Es gibt noch ein paar Beispiele für Simulationen an folgenden Standorten:

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

Angepasst an Bolker 2009 Ökologische Modelle und Daten in R. Sie müssen die Stärke des Trends (dh der Steigung) deklarieren, den Sie testen möchten. Intuitiv erfordert ein starker Trend und eine geringe Variabilität eine kleine Stichprobengröße, ein schwacher Trend und eine große Variabilität erfordern eine große Stichprobengröße.

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

Sie können auch simulieren, was der minimale Trend ist, den Sie für eine bestimmte Stichprobengröße testen können, wie im Buch gezeigt

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)