# Test examples" In simr: Power Analysis for Generalised Linear Mixed Models by Simulation

This vignette provides examples of some of the hypothesis tests that can be specified in `simr`. The function `doTest` can be used to apply a test to an input model, which lets you check that the test works before running a power simulation.

Documentation for the test specification functions can be found in the help system at `?tests`.

```library(simr)
```
```simrOptions(progress=FALSE)
```

## Binomial GLMM with a categorical predictor

The first example comes from the help page for `glmer`. The data frame `cbpp` contains data on contagious bovine pleuropneumonia. An observation variable is added to allow for overdispersion. Note that the response is specified using `cbind` --- `simr` expects a binomial model to be in this form.

```cbpp\$obs <- 1:nrow(cbpp)
gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd) + (1|obs), data=cbpp,
family=binomial)
summary(gm1)\$coef
```

Note that `period` is a categorical variable with 4 levels, which enters the model as 3 dummy variables. To test all 3 dummy variables simultaneously, you can use a likelihood ratio test.

```doTest(gm1, fixed("period", "lr"))
```

If you were (for some reason) especially interested in the significance for the dummy variable `period2` you could use a z-test. This test uses the value `Pr(>|z|)` reported in the summary above.

```doTest(gm1, fixed("period2", "z"))
```

Suppose your model also has a continuous predictor. You can use `fixed` to choose which fixed effect to apply tests to.

```gm2 <- glmer(cbind(incidence, size - incidence) ~ period + size + (1 | herd), data=cbpp,
family=binomial)
doTest(gm2, fixed("size", "z"))
```

Once you have chosen your tests, you can run a power analysis by replacing `doTest` with `powerSim`. Don't forget to specify an appropriate effect size.

```fixef(gm2)["size"] <- 0.05
powerSim(gm2, fixed("size", "z"), nsim=50)
```

## Models with interaction or quadratic terms

As your models become more complex, it can be easier to explicitly specify your null hypothesis using the `compare` functions.

### Cake

This example uses the `cake` dataset.

```fm1 <- lmer(angle ~ recipe * temp + (1|recipe:replicate), data=cake, REML=FALSE)
```

Main effects should not be tested when they appear in an interaction term. Using the `fcompare` function, we can specify a comparison with a simpler model (without having to re-type the random effects specification).

```doTest(fm1, fcompare(~ recipe + temp))
```

This also works for polynomial terms:

```fm2 <- lmer(angle ~ recipe + poly(temp, 2) + (1|recipe:replicate), data=cake, REML=FALSE)
summary(fm2)\$coef
doTest(fm2, fcompare(~ recipe + temp))
```

### Budworms

We can do similar things with the `budworm` data in the `pbkrtest` package.

```data(budworm, package="pbkrtest")
bw1 <- glm(cbind(ndead, ntotal-ndead) ~ dose*sex, family="binomial", data=budworm)
summary(bw1)\$coef
```

Of course we don't want to retype the `cbind` boilerplate:

```doTest(bw1, compare(. ~ dose + sex))
```

Since `dose` is continous and `sex` is binary we could also use a Z-test on the single interaction term.

```doTest(bw1, fixed("dose:sexmale", "z"))
```

## Single random effects

The `random` function gives you access to tests from the `RLRsim` package. No additional arguments are needed for a single random effect.

```re1 <- lmer(Yield ~ 1|Batch, data=Dyestuff)
doTest(re1, random())
```

## Multiple random effects

Where the model has multiple random effects, `compare` can be used to test alternative specifications.

```fm1 <- lmer(Reaction ~ Days + (Days | Subject), data=sleepstudy)
```
```doTest(fm1, compare( ~ Days + (1 | Subject)))
```

The LRT is fast but only approximate. In fact, when testing random effects, the test is conservative[^1] because the null hypothesis is at a boundary of the parameter space. This means that you will underestimate power if you use the LRT. For more accurate results you can use `compare` with a parametric bootstrap test from the `pbkrtest` package. These can be quite slow, so you may want to use the LRT to exploring designs, and then double check with the parametric bootstrap.

```doTest(fm1, compare( ~ Days + (1 | Subject), "pb"))
```

Note that the shortcut `rcompare` saves you from retyping the fixed effect specification.

```doTest(fm1, rcompare( ~ (1 | Subject), "pb"))
```

## A note about errors during simulation

During a simulation study, some iterations may fail due to some sort of error. When this happens, `simr` treats that iteration as a failed (i.e. not significant) test. In the following example there are 50 simulations, with 14 successes, 34 failures, and 2 non-results. The power is calculated as 14/50, i.e. 28%:

```binFit <- glm(formula = cbind(z, 10 - z) ~ x + g, family = binomial, data = simdata)

poisSim <- glm(formula = z ~ x + g, family = poisson, data = simdata)
coef(poisSim)[1:2] <- c(1, -0.05)

powerSim(binFit, sim=poisSim, nsim=50, seed=1)
```

Rather than interrupting part-way through an analysis, `simr` traps and logs errors and warnings. You can access these logs using `\$warnings` and `\$errors`. If you didn't assign your analysis to a variable, you can recover it with the `lastResult` function.

```ps <- lastResult()
ps\$errors
```

[^1]: See, e.g., Pinheiro, J.C., Bates D.M. (2000) Mixed-Effects Models in S and S-PLUS, Springer, New York (p84).

## Try the simr package in your browser

Any scripts or data that you put into this service are public.

simr documentation built on Jan. 29, 2019, 9:04 a.m.