# Post hoc analysis for linear mixed model with nested effects

I'm having trouble in R with my Linear Mixed-Effects Model. I'm working with yeast in nectar. This is a part of my data just so you can see what is going on:

For the condition sucrose, I have 4 different samples (you can only see data for sample 4 here). For each sample I did 2 replica's (so replica is either 1 or 2). sp tells you which species it is and condition tells you whether the two yeast species were mixed together or just grew alone (single). I linked the two variables condition and sp together in treatment. host specifies the host plant of the species and cells1 is the number of yeast cells.

It is the number of cells (cells1) that I want to compare for the different treatments. So I started off by making a mixed model with nested effects.

suc <- read.csv(file=file.choose(),header = TRUE,sep = ";")
attach(suc)
## load packages 'lme4', 'lsmeans', 'pbkrtest and 'Rcpp'
fit1 <- lmer(cells1~treatment+sp+treatment:sp+
(1|cont/replica)+(1|replica/Sample)+(1|Sample/host), suc)


Next, I wanted to do a post hoc test. TukeyHSD wouldn't work. Error said something about not being able to use it with lmer. So after doing some research, I used the function lsmeans.

library("lsmeans")


When looking at the output, I get NA as outcome everywhere and I have no clue what is wrong or how to resolve this.

Does anyone know how I can fix this?

Your treatment variable represents the interaction between condition and sp, so putting treatment and sp in the model is redundant. Since the difficulties you're having are with the fixed-effect model, you can diagnose/debug more simply by working with lm() until you can have a workable model. I would try

 fit1 <- lm(cells1~condition*sp,data=suc)


(which is equivalent to a response of ~condition+sp+condition:sp), check that all the parameters are estimated, and then move on to the mixed model.

Using attach() is often a bad idea.

I don't know what cont is.

Your random effect term should probably (?) be 1|Sample/replica (i.e. replica nested within Sample). You may need to consider a random-slopes model, although it could be too difficult to fit.

• I counted the cells in 16 equal squares in a roster to get a correcter idea of how many there where. I indicate this with the variable 'cont'. – Sara S Mar 29 '15 at 21:23

I got the result I wanted by using the following commands:

suc$logcells <- log(Sucrose$cells1)
# load packages 'lmerTest', 'lsmeans', 'pbkrtest and 'Rcpp'
fit2 <- lmer(logcells~treatment + host + treatment:host +
(1|cont/replica)+(1|replica/Sample) + (1|Sample),
data=Sucrose)
summary(fit2)
anova(fit2)
ls_fit2 <- lsmeans(fit2, test.effs="treatment:host", adjust="Tukey")
ls_fit2
T_fit3 <- difflsmeans(fit2, test.effs="treatment:host", adjust="Tukey")
T_fit3


The main issues were that I first needed to transform my data logarithmically and that some commands did not work in combination with certain packages.

• Well, OK, if you think those are the right results. But you did not use the lsmeans package at all, you used the lsmeans function in lmerTest. So you probably shouldn't have attached lsmeans. – Russ Lenth Apr 2 '15 at 1:39