I have a pretty involved solution that requires you specify some hyperparameters of your model. Not being an expert in the relation between snake length and age, so I picked some values that seemed vaguely possible. But, by all means feel to tweak them. This series of functions is adapted and expanded from [Gelman and Hill (2007)][1].
First, we need a function to simulate data sets with the properties you want. Here is where I "invented" some numbers (under the comment `#Specifying hyperparameters`. Feel free to adjust accordingly based on your knowledge base and expectations.
#Creating fake dataframe for simulations
fake<-function(J, K){
snake<-(1:(J*K)) #creates a snake or case id
habitat<-rep(1:J, K) #creates a balanced grouping variable
age<-rnorm(J*K, 2.5, .75) #creates a random normal age variable
#Specifying hyperparameters
mu.a.true<-12 #grand intercept
mu.g1.true<-3 #grand age slope
sigma.y.true<-1.65 #residual sd
sigma.a.true<-2.30 #intercept sd
sigma.g1.true<-2.15 #age slope sd
#habitat-level parameters
a.true<-rnorm(J, mu.a.true, sigma.a.true) #vector of random intercepts
g1.true<-rnorm(J, mu.g1.true, sigma.g1.true) #vector of random slopes
#dependent variable
y<-rnorm(J*K, a.true[habitat]+g1.true[habitat]*age, sigma.y.true)
return(data.frame(snake, y, habitat, age))
}
and if you run `print(fake(5,5))` you get the following simulated data (where y = length):
> print(fake(5,5))
snake y habitat age
1 1 21.85790 1 1.9030948
2 2 27.73040 2 2.4382916
3 3 11.60743 3 0.7640231
4 4 15.63832 4 2.9965471
5 5 26.83866 5 1.8319489
6 6 24.65482 1 3.0666732
7 7 27.42804 2 2.2352739
8 8 16.98225 3 3.5354592
9 9 15.32236 4 3.0284722
10 10 25.22085 5 2.4662031
11 11 18.39988 1 2.2474545
12 12 21.85225 2 1.2652154
13 13 14.12972 3 2.6137810
14 14 17.69987 4 1.3900062
15 15 33.13150 5 3.2806024
16 16 28.80915 1 4.6216763
17 17 35.61110 2 2.5738612
18 18 15.78296 3 2.5815051
19 19 17.76041 4 3.2101400
20 20 29.23588 5 2.5638092
21 21 21.51020 1 2.7942239
22 22 32.10491 2 2.6981291
23 23 10.05619 3 1.7800981
24 24 16.96727 4 3.0558285
25 25 34.61646 5 3.3928710
Up next? We need to simulate a bunch of data and see how many data sets "recover" a significant effect for age given the hyperparameters we input in our data simulation function. Again, those hyperparameters can be adjusted to better fit your expectations.
mixed.power<-function(J, K, n.sims=1000){
signif<-rep(NA, n.sims) #note that you can specify number of simulations - default is 1000
pb<-winProgressBar(title="Progress", min=0, max=100, width=300) #if you want to watch the progress
require(lme4)
require(arm)
for(s in 1:n.sims){
fake.data<-fake(J, K) #calls in data simulation function
lme.power<-lmer(y~1+age+(1+age|habitat), data=fake.data) #estimates mixed effect model using each simulated dataset
theta.hat<-fixef(lme.power)["age"] #saves age coefficients from each simulated dataset
theta.se<-se.fixef(lme.power)["age"] #saves standard error of age coefficients from each simulated dataset
signif[s]<-ifelse((abs(theta.hat)-qt(.975, J*K-1)*theta.se)>0, 1, 0)#assigns value of 1 to significant coefficients 0 to ns coefficients
setWinProgressBar(pb, s/n.sims*100, title=paste(round(s/n.sims*100, 0), "% done"))
}
close(pb)
power<-mean(signif, na.rm=T)#calculates proportion of significant models out of # of simulated datasets...
return(power)
}
You may get some warnings with these codes when models don't converge in your simulations. I find that happens more often when dealing with small samples like the ones you are proposing. (I tend to ignore these convergence warnings)
You can pause here and calculate power for a given number of habitats (J) with a given number of snakes per habitat (K). In that case you can run `mixed.power(5,10)` and see the following returned:
> mixed.power(5, 10)
[1] 0.804
Or perhaps more instructively you can go all out and create a power plot to get a sense of how many snakes per habitat you should consider gathering. Here is the function for that:
graph.power<-function(J, max.K){
require(ggplot2)
KK<-seq(3, max.K, by=1) #will ieterate from 3 cases to the max.K specified in the function above
Y<-rep(NA, length(KK)) #Empty vector to contain power estimates from mixed.power () function
JJ<-rep(J, length(KK)) #Repeates the number of habitats you specified to equal the length of the KK vector
pb2<-winProgressBar(title="Overall Progress", min=0, max=100, width=300) #for an overall progress bar
for(i in 1:length(KK)){
Y[i]<-mixed.power(JJ[i], KK[i]) #runs mixed power function and stores results one at a time
setWinProgressBar(pb2, i/length(KK)*100, title=paste(round(i/length(KK)*100, 0), "% Overall complete"))
}
close(pb2)
DF<-as.data.frame(cbind(KK, Y)) #creates data frame for plotting
colnames(DF)<-c("KK", "Y")
#Code below provides a basic power plot with a horizontal line at .80
g1<-ggplot(aes(x=KK, y=Y), data=DF)
g2<-g1+geom_smooth(se=F)+geom_hline(yintercept = .8, lty="dashed", col="red", lwd=2)
g3<-g2+xlab("Sample Size per Habitat")+ylab("Power")
g4<-g3+ggtitle("Power Analysis")
return(g4)
}
When you use this function, the only difference is that now the second number is the maximum number of snakes you want to consider examining in a power analysis. So in this example, let's say I wanted to look at a power curve that considered up to 15 snakes per habitat (holding the number of habitats constant at 5). Here is what that would look like: `graph.power(5,15)`... which would return the following image:
[![enter image description here][2]][2]
Now you have the ability to assess the sample size at which the power curve crosses .80 (~9 snakes per 5 habitats in this case).
This is NOT the only way to approach a power analysis in a mixed model, but I tend to like it as it gives me a great deal of control over the various hyperparameters and allows me to plot various power curves to get a sense of how changing those hyperparameters in addition to sample and cluster size change the power of my model.
A note: if you want to play around to get a sense of how the functions work or to try to identify some good starting values you may want to lower the number of simulations in the `mixed.power()` function to reduce computation time prior to running your "official" power analyses with the larger set of simulated data sets.
[1]: https://g.co/kgs/M5ELVg
[2]: https://i.sstatic.net/gLzlZ.jpg