2
$\begingroup$

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
Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

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

Improve this answer
$\endgroup$
1
  • $\begingroup$ @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? $\endgroup$
    – u01arc7
    Nov 15, 2014 at 17:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.