lme4_pvals.Rmd

# Exploring pvalue options for lme4 models
Julin Maloof
Initiated Jan 21, 2014

## The problem
lme4 does not report p-values for fixed or random effects.  The package that Dan Chitwood used, languageR, is out of date, has been shown to give unreliable results, and is picky about the lme4 version that it uses.  There are several alternatives that we will explore here.

## load packages
lme4 for model fitting.  lmerTest and car for various hypothesis testing functions.
lme4 > 1.0 is required.  If you haven't updated your R for a while you will need to.

```{r}
library(lme4)
library(lmerTest)
library(car)
```


## Read in the data
Following Dan's example, we will read in the data and transform the abs_stom trait to give it a more normal distribution

```{r read_data}
stomdata <- read.delim("Modeling_example.txt")
```

```{r}
stomdata$trans_abs_stom <- sqrt(stomdata$abs_stom)
```
```{r}
head(stomdata)
summary(stomdata)
```

## create the full model
Next we fit the fully specified model.  Note that by adding the `arugment data=stomdata` to the lmer call that we can simplify the term specication (using names internal to the `stomdata` data frame)
```{r fit_models}
model1 <- lmer(trans_abs_stom ~ il + (1|tray) + (1|row) + (1|col),data=stomdata)
```

## Assessing signficance of random effects
We can use `rand()` in the lmerTest package to test the signficance of the random effects.
```{r random}
rand(model1) #depends on lmerTest package
```

## Alternative random effects testing.
If you prefer to do this by comparing models, that also can be accomplished.  lmerTest has an updated anova function.
```{r random2}
model3 <- update(model1, . ~ . - (1|row)) # remove the row term from the model.
anova(model1,model3)
```


## Assessing signficance of fixed effects
The anova function from lmerTest can also be used to test the fixed effect term(s). 
```{r fixed}
anova(model1) #default is like a SAS type 3.
anova(model1,type=1) #Alternative: SAS type 1.
anova(model1,ddf="Kenward-Roger") #Alternative way to calculate the degrees of freedom.
```

## Assessing signficance of fixed effect factor levels
What if you want to test the signficance of each factor level (each IL in this case...)?  Once lmerTest has been loaded, then the summary function works like you would want it to.
```{r factor}
summary(model1) # gives p-values for standard lme comparisions (each against the reference)
```

## Making specific comparisions among factor levels
summary() compared everything to the reference level (M82 in this case).  You can use the functions in the car package to make other comparisons (ie IL vs IL).
```{r car_Example}
linearHypothesis(model1,"ilIL_1.1.2 = ilIL_1.1") #compare IL_1.1.2 to IL_1.1.  Note that you have to use the factor names as they are listed in the summary table.  Hence "ilIL_1.1.2..."
```

to test the hypothesis that multiple ILs are equivalent:

```{r}
linearHypothesis(model1,c("ilIL_1.1.2 - ilIL_1.1","ilIL_9.3 - ilIL_9.3.1"))
```
there are many other ways to specify the comparisions.  see ```?linearHypothesis``` for more details.

## Confidence Intervals

We can calculate confidence intervals on our coefficients using functions now available in the lme4 package itself.  Two methods are shown below.

```{r}
model1.confint <- confint(model1) #uses profile method.  This is a likelihood based method.  See ?confint.merMod for details.  Takes 138 seconds on 2013 Macbook Pro
model1.confint.boot <- confint(model1,method="boot") #bootstrapped confidence intervals. Takes 130 seconds on 2013 Macbook Pro
```

#compare them

```{r}
plot(model1.confint[-1:-5,1],model1.confint.boot[-1:-5,1]) #pretty similar in this case.
```