0
$\begingroup$

I have the following data (just an small example here). I want to know effect of time and ecotype over height. subdata <- id time ecotype height 1 17 a 73.692 1 22 a 213.010 1 25 a 343.700 1 28 a 663.030 109 32 a 1267.300 109 17 b 47.445 109 22 b 148.050 109 25 b 280.570 109 28 b 509.000 109 32 b 954.050 121 17 c 33.972 121 22 c 71.235 121 25 c 130.650 121 28 c 220.190 121 32 c 452.740 133 17 d 38.365 133 22 d 84.068 133 25 d 137.560 133 28 d 255.390 133 32 d 426.970

I have fitted the following model

m1 <- lme(height ~ ecotype*time + (time | id))

and now I want to know the Power of this analysis in order to do that I use simulation

time <- (subdata$time) ecotype <- as.factor(subdata$ecotype) id <- subdata$id height=subdata$height

m2 <- lmer(height ~ ecotype * time + (time| id), subdata) s2 <- simulate(m2) beta.hat=fixef(m2) se=sqrt(diag(vcov(m2)))

k <- c() B=1000 tstar=rep(0,B) set.seed(781) for(b in 1:B) { ystar=drop(simulate(m2)) ostar=lmer(ystar$sim_1~ecotype +(time|id)) for(i in 2:19) { k <- append(k, ((fixef(ostar)[i]-beta.hat[i])/sqrt(vcov(ostar)[i,i]))) }}`

The problem I have is I don't know how to calculate an overall p-value to say this model's power is low or high. Normally, the examples on the web always have variables (e.g ecotype) with two classes so the results of the simulation will always produce one value. However, in my case I have n number of ecotypes therefore more than one output:

fixef(ostar) (Intercept) ecotypea ecotypeb ecotypec ecotyped 78.98846 -22.39385 -37.27676 -38.66168

All I can think is to extract all the fixed effects estimates and the cv and compare them against, get a mean and compare it against a threshold (e.g. 0.05). However, I don’t know whether that’s right. Please can you give me a hand on this Please can you suggest how can I do that?

Thanks

$\endgroup$
0
$\begingroup$

I have solved my question by doing the following:

library(lmerTest)

simulateFixedEffect <- function(model,ecotype, time, id, B) { Fvalue=rep(0,B); set.seed(781); for(b in 1:B) { s1=drop(simulate(model)); ms=lmer(s1$sim_1 ~ ecotype*time + (time|id)); Pvalue[b] <- anova(ms)[[6]][3]; } return(Pvalue);
}

plotValue <- function(Fvalue, alpha = 0.05) { mean(Fvalue>alpha); hist(Fvalue,30,prob=T,xlab="",ylab="",main="") x<-seq(0,15,0.1) lines(x,df(x,1,80),col="red") grid() }

b <- calculateT(m1); Pvalue <- simulateFixedEffect(m1,ecotype, time, id, 10); plotValue(Pvalue);

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.