# Graphical representation of fixed effects from lmer

I have run a lmer model in R:

   M25<-lmer(sqrtAbund~TP1 + Temp1 + CO_21 +
TP1:CO_21 + Temp1:CO_21 +
Temp1:CO_21:TP1 + (1|pseudo), REML=FALSE, data=sqrtCyano)


Abund (continuous response) => microbial gene abundance (sqr-transformed) TP1 (factor; fixed effect) => time point - 3 levels (day 0, 7 and 28) Temp1 (factor; fixed effect) => temperature - 2 levels (12 and 16 degrees) CO_21 (factor; fixed effect) => carbon dioxide - 2 levels (380 and 750 ppmv) pseudo (factor; random effect) => this effect is based on pseudoreplication of the independent replicates (n=3).

From exploratory analysis there seems to be a seasonal effect on gene abundance i.e samples taken in spring have higher abundance than when taken in winter. From each independent replicate, 4 samples were taken. I have coded pseudo variable to explain the variation given that samples taken within an independent replicate are more similar than between independent replicates within a given treatment (i.e. CO2 + Temp, CO2:Temp)

The question I am aiming to answer with this model:

Is there an effect of CO2 and/or temperature on microbial gene abundance?

Treatments = 4, reps = 3
control   (CO2, 380; Temp, 12)
high temp (CO2, 380; Temp, 16)
high CO2  (CO2, 750; Temp, 12)
combined  (CO2, 750; Temp, 16)


sampling was carried out at 3 time points: day 0, day 7, day 28 and subsampled n=4

I finally came to the simplest model based on extracted parameter specific p-vals using the code:

coefs <- data.frame(coef(summary(M25)))
coefs$p.z <- 2*(1 - pnorm(abs(coefs$t.value)))
coefs


Based on these parameter values, I was happy that they corresponded with simple excel graphs of the data and thereby being the main drivers of the system. I carried out model validation and again, happy that this model is a good fit.

However, I did want to plot the predicted model using the following code:

    y<-(sqrtCyano$sqrtAbund) fit<-M25 pred<-predict(fit) plot(y, pred, xlim=range(c(y,pred)), ylim=range(c(y,pred)), xlab="observed", ylab="predicted") abline(0,1, lwd=2, col=8) #Line [7] fit2 <- lmer(pred ~ y+ (1|pseudo)) lgd <- c( paste("R^2 =", round(summary(fit2)$r.squared,3)),
paste("Offset =", round(coef(fit2)[1],3)),
paste("Slope =", round(coef(fit2)[2],3))
)
legend("topleft", legend=lgd)
abline(fit2, lwd=2)
legend("bottomright", legend=c("predicted ~ observed", "1:1"), col=c(1,8), lty=1, lwd=2)


but, as soon as I get to line [7] it throws out the error:

1: Some predictor variables are on very different scales: consider rescaling 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue - Rescale variables?

Is there another way to plot the predicted model?

I would like to graphically represent the model and from reading around, the packages ggplot and multcomp have been mentioned. I used the script from here but to no avail as it throws me errors because I have multiple factors to include and i'm not sure how to code it in a way that I can bind these factors together.

I am specifically refering to

as.data.frame(confint(glht(model, mcp(...= "Tukey")))\$confint)


(...) where I have to put in the model fixed effects - I have 3 single and 3 interaction terms?

What i would like to know:

1. Is there another way to plot the predicted lmer model?
2. How would I graphically represent the fixed effects evaluation?
3. Just a side question...would you recommend fitting using REML based on my description of the random factor, pseudo?

output summary from M25:

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: sqrtAbund ~ TP1 + Temp1 + CO_21 + TP1:CO_21 + Temp1:CO_21 + Temp1:CO_21:TP1 +      (1 | pseudo)
Data: sqrtCyano

AIC      BIC   logLik deviance df.resid
3207.4   3248.6  -1589.7   3179.4      126

Scaled residuals:
Min      1Q  Median      3Q     Max
-2.4108 -0.5841  0.0322  0.4319  3.6544

Random effects:
Groups   Name        Variance  Std.Dev.
pseudo   (Intercept) 187788861 13704
Residual             313451434 17705
Number of obs: 140, groups:  pseudo, 36

Fixed effects:
Estimate Std. Error t value
(Intercept)               26271       9562   2.747
TP17                      25965      13422   1.934
TP128                     28733      13422   2.141
Temp116                   15044      13422   1.121
CO_21750                  47836      13422   3.564
TP17:CO_21750            -11185      18910  -0.591
TP128:CO_21750           -89707      18982  -4.726
Temp116:CO_21750         -57713      18910  -3.052
TP17:Temp116:CO_21380    -18675      18910  -0.988
TP128:Temp116:CO_21380     4652      18982   0.245
TP17:Temp116:CO_21750     -7929      18910  -0.419
TP128:Temp116:CO_21750    64025      18910   3.386


## 1 Answer

The simply way to get predicted values is to provide a data.frame with just the minimal fixed effects you want to plot. Then, when you run predict using that as the newdata set the re.form argument to NA.

see ?predict.merMod

(BTW, I think you want 1|subject or something for your random effect. Your replications are nested within subject, not the subjects in replications?)

• @John...thanks very much for your response, I will try this and yes, reps (n=4) are nested within the independent replicates (n=3) for each treatment (n=4). I have organised the data frame so that the pseudoreps relate to the sample number taken from the treatment. Does this mean that the random factor should be coded something like: pseudo <- replicate:treatment? Nov 15, 2014 at 17:59