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:
- 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