# How to build a linear mixed-effects model in R?

I asked a question but it was a bit long and confusing so I will attempt to keep this shorter and simple (original post https://stats.stackexchange.com/questions/24971/mixed-effects-model-equations)

I am basically looking to see which factors (e.g. A, B, C such as habitat type, site, disturbance rate) most affect the response value (y). I am finally understanding the concept of models and how to compare AIC values to see which combination of factors best explain y (have greater effect on), but I am new to R so wonder if my basic coding is correct.

Firstly I have all my data in a spreadsheet saved as a .csv file, so I read this file into R. I also open the package lme4.

Then I was thinking of using the following code (although letters would be the headings fo each set of values)

m1<-lmer(y ~ A + (1|E), REML=FALSE)
m2<-lmer(y ~ B + (1|E), REML=FALSE)
m3<-lmer(y ~ A + B + (1|E), REML=FALSE)
m4<-lmer(y ~ A + C + (1|E), REML=FALSE)
m5<-lmer(y ~ B + C + (1|E), REML=FALSE)
m6<-lmer(y ~ A + B + C + (1|E), REML=FALSE)
m4<-lmer(y ~ A + D + (1|E), REML=FALSE)
m5<-lmer(y ~ B + D + (1|E), REML=FALSE)
m6<-lmer(y ~ A + B + D + (1|E), REML=FALSE)
m7<-lmer(y ~ C + (1|E), REML=FALSE)
m8<-lmer(y ~ D + (1|E), REML=FALSE)


My basic questions are:

1. Let us say that C and D are similar factors, and derived from the same data and I do not think it useful to see if they have additive effects, is it okay not to put them in the same model, i.e. mix and match as see fit, or should all combinations be incorporated?

2. I chose ML because not all sets of values are mixed and matched. Is this correct or should I use REML?

3. Some sets of values are non-continuous, e.g. brood number or habitat type (either 1 or 2). Should I be letting R know this. Someone mentioned coding each factor but I have no clue how to do this, or is the csv file enough and then perhaps code to let R know about these particular ones?

4. Lastly, is there also a way to see p-values to see if the factor has a significant effect rather than just being the best fit to explain y?

I hope that this is more simple than my first post and easier to answer all at once. I really wanted to run this model today, but don't want to have to redo it all finding that I made a simple mistake.

• What is E? The way you've currently specified your models, lmer will fit a different intercept for each level of E. What kind of variable is it, and why do you think it will be useful to fit these random intercepts for the level of E? Choosing ML instead of REML allows you to compare model fits for different fixed-effects structures (eg. by doing anova(m1, m2). All your models currently have different fixed-effects structures, and the same random effect structures, so this seems like the right choice. – Marius Mar 22 '12 at 0:35
• E is the ID of the Brood, and I have numerous non-independent data points for each brood so this needs to be accounted for. I am REALLY new to models, but have been trying to grasp them, but I literally have no time left. I only understand how they can be used in terms of looking at how different effects contribute to the measurement (y) and then how to decide which models are the best using AIC and AIC weights to realise which models best explain the values of measurement 't'. For more detail you could look at my other post. My sample size is so small that no other stats will work. – Dragonwalker Mar 22 '12 at 1:29

1. It's entirely up to you as to whether you include both factors in the same model or not. But why not try it and see if you get a significantly better fit with both in than with just one in?

2. REML works with unbalanced and incomplete designs too. I'd go with REML to reduce the bias in the variance estimates and eliminate the bias in the covariance parameters.

3. x <- as.factor(z) turns z into a factor. You can of course do DF$x <- as.factor(DF$x).

4. anova(m1, m3) will test for the significance of the terms left out of the larger model. The models have to be nested for this to work.

Edit: The comments have made me realize my answer was way too terse, so I'm adding some sample code.

The code is not doing the full model that you are doing, it's just to illustrate syntax and what happens with ANOVA:

# Construct sample data; E(y) is a function only of x1
x1 <- c("A","A","A","B","B","C","D","D","D","D")
x2 <- c("A","B","C","A","B","C","A","B","C","A")
y <- rnorm(c(0,0,0,1,1,2,3,3,3,3))  # Values for E(y) match w/ x1

# Construct data frame
df <- data.frame(list(y=y, x1=x1, x2=x2))

# Convert x1, x2 to factors
df$x1 <- as.factor(df$x1)
df$x2 <- as.factor(df$x2)

# Run regressions and perform ANOVA to evaluate effect of factor x2
m1 <- lm(y~x1, data=df)
m2 <- lm(y~x1+x2, data=df)

> anova(m1,m2)
Analysis of Variance Table

Model 1: y ~ x1
Model 2: y ~ x1 + x2
Res.Df     RSS Df Sum of Sq     F Pr(>F)
1      6 12.9004
2      4  5.3241  2    7.5763 2.846 0.1703


The "PR(>F)" column gives the p-value associated with the F-test of whether factor x2 is significant.

• Thank you jbowman. I just have two questions in regards to your answer. for 3. do I have to do that for every column of data, e.g. for the random effect, and the other effects (e.g. y or A or B), or just those that have the arbitrary numbers for ID? 4. I am not sure what you mean. What would the results actually show in regards to m1 and m3? I am sorry that I don't know this stuff. :( I am not sure either if my models are nested. I have seen a script: >pvals.fnc(test) I think test is the model that you wish the p-vals for. Thanks – Dragonwalker Mar 21 '12 at 18:25
• Also, I don't understand what DF$x means, how would this look if the name is brood rather than x? Do you type DF$brood. Does this ensure that it is still called brood? I am sorry, but coding alludes me. – Dragonwalker Mar 21 '12 at 18:34
• I don't know why my comment looks like that and I can't erase it. I basically don't understand what DS$means. Could I then use the name of the factor rather than a letter? – Dragonwalker Mar 21 '12 at 18:35 • @Dragonwalker That's pretty basic R command to address a named column in a data.frame. (I fixed your earlier comment where text between the $'s has been parsed as a $\LaTeX$ expression.) – chl Mar 21 '12 at 19:01
• So if I can just make it clear. I would do that for each column in my data file. So for example if the column of data is brood then I could type 'DF$brood<- as.factor(DF$Brood)' and this would turn brood into a factor but then it would still be called brood for when I type it into the model (or would it be called DF$brood?). For those that are not arbitrary but contain vales I would then use: 'DF$Habitat<-factor(DF\$Habitat)' – Dragonwalker Mar 21 '12 at 19:56