Como simular uma análise de energia personalizada de um modelo lm (usando R)

Seguindo as perguntas recentes que tivemos aqui .

Eu esperava saber se alguém havia se deparado ou pode compartilhar o código R para realizar uma análise de energia personalizada com base na simulação de um modelo linear?

Mais tarde, obviamente, gostaria de estendê-lo a modelos mais complexos, mas lm parece ser o lugar certo para começar. Obrigado.

Não sei se você precisa de simulação para um modelo de regressão simples. Por exemplo, consulte o documento Energia portátil . Para modelos mais complexos, especificamente efeitos mistos, o pacote pamm em R executa análises de energia por meio de simulações. Veja também o post de Todd Jobe, que tem código R para simulação.

Aqui estão algumas fontes de código de simulação em R. Não tenho certeza se algum aborda especificamente modelos lineares, mas talvez eles forneçam um exemplo suficiente para entender a essência:

• Benjamin Bolker escreveu um grande livro Ecológica de Dados e Modelos com R . Um rascunho inicial de todo o livro, juntamente com o código Sweave, está disponível online. O capítulo 5 trata da análise e simulação de energia.

Existem outros exemplos de simulação nos seguintes sites:

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

Adaptado dos modelos e dados ecológicos de Bolker 2009 em R. Você precisa declarar a força da tendência (ou seja, a inclinação) que deseja testar. Intuitivamente, uma forte tendência e baixa variabilidade exigirão um tamanho pequeno da amostra, uma tendência fraca e grande variabilidade exigirão um tamanho grande da amostra.

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

Você também pode simular qual é a tendência mínima que pode testar para um determinado tamanho de amostra, conforme mostrado no livro

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)