```
library(lme4)
library(ggplot2); theme_set(theme_bw())
library(broom)
library(dotwhisker)
```

There are a lot of packages for using the built-in `simulate.merMod()`

function in power calculations: e.g.Â see the reverse-depends/imports list at the CRAN lme4 page. However, if youâ€™re an intermediate R user it shouldnâ€™t be too hard to simulate data for your particular case, with whatever experimental/observational designs, degrees of missingness, etc., you like.

Letâ€™s suppose we start with 50 locations, each observed 4 times. There is an overall temporal trend, variation in density and trend among sites due to a single covariate, as well as random variation in time and random variation among sites in slope and intercept.

```
nsites <- 50
ntimes <- 4
dd <- expand.grid(site=factor(1:nsites),time=1:ntimes)
set.seed(1001)
site_cov <- data.frame(site=dd$site,x=rnorm(nsites,mean=1,sd=0.1))
dd <- merge(dd,site_cov)
dd$dens <- 1
```

Do some preliminary model setup so we can figure out the order in which the random effects occur (`lme4`

internally reorders the random effects, in decreasing order of the number of levels of the grouping variable). Weâ€™re including `(time|site)`

to allow for random slopes among sites; `(1|time)`

to include some population-wide variation around the otherwise linear trend with time (4 time periods is probably too small to reliabily fit this random effect, but it is reasonable to include it in the simulation); `(1|site:time)`

is an observation-level random effect, which will include overdispersion in the model (if our response distribution has a scale/variance/overdispersion parameter, e.g.Â negative binomial, then we should skip this term).

```
gForm <- dens~time*x+(time|site) + (1|time) + (1|site:time)
gg <- glFormula(gForm,data=dd,family=poisson)
##
names(gg$reTrms)
```

```
## [1] "Zt" "theta" "Lind" "Gp" "lower" "Lambdat" "flist"
## [8] "cnms" "Ztlist"
```

`gg$reTrms$cnms`

```
## $`site:time`
## [1] "(Intercept)"
##
## $site
## [1] "(Intercept)" "time"
##
## $time
## [1] "(Intercept)"
```

`colnames(gg$X) ## fixed-effect parameters`

`## [1] "(Intercept)" "time" "x" "time:x"`

```
true_beta <- c(1,0.1,-1,0.2)
names(true_beta) <- colnames(gg$X)
dd$dens <- simulate(gForm[-2],
newdata=dd,
newparam=list(beta=true_beta,
theta=0.5*c(1,1,0,0.1,1)),
family=poisson)[[1]]
```

`## theta parameter vector not named: assuming same order as internal vector`

```
ggplot(dd,aes(time,dens,group=site))+geom_point()+geom_line()+
scale_y_log10()
```

```
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis
```