Solved – Graphical representation of fixed effects from lmer

data visualizationlme4-nlmemodelpredictor

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

Best 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?)

Related Question