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  fit2 <- lmer(pred ~ y+ (1|pseudo)) lgd <- c( paste("R^2 =", round(summary(fit2)$r.squared,3)), paste("Offset =", round(coef(fit2),3)), paste("Slope =", round(coef(fit2),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  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:
- Is there another way to plot the predicted lmer model?
- How would I graphically represent the fixed effects evaluation?
- 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