I am trying to evaluate the effect of a treatment (2 levels: mock or TuMV) on 4 different plant genotypes (Col-0, ATG5, GIN2, HXK1) by means of linear mixed models (
lme). I had seen in advance that one group (ATG5-TuMV) has a much higher dispersion than the rest, but without outliers.
my.1rst.attempt<-lme(TOI ~ treatment*genotype, random= ~1 | ID, data = my.data)
It is quite evident for me that a particular group has a very different variance structure:
r<-residuals(my.1rst.attempt,type="p") boxplot(r~my.data$genxtrat,xlab="TOI", ylab="Standardized residuals")
qqnorm(my.1rst.attempt, abline = c(0,1)) #QQ plot:
So, I tried to improve that model allowing different variance structures by group:
VarIdent.attempt<-update(my.1rst.attempt,weights=varIdent(form= ~ 1|treatment*genotype))
and obtained a model that seems to be much better:
A Shapiro-Wilks test on the residuals improved from extremely significant (p-value = 1.866e-08) to marginally significant (p-value = 0.044).
The new model gave me an improvement of AIC, BIC (p <.0001).
However, when assessing the observed vs fitted values, the simpler and worse model gave me a near-perfect fit:
p<-ggplot(my.data, aes(x =TOI , y = predict(my.first.attempt), colour =genxtrat)) + geom_point(size=2) p + geom_abline(intercept = 0, slope =1) + ggtitle("Predict vs observed")
With very high correlation = 0.9980763
Which is seen for the subset of data: plot(my.first.attempt, TOI ~ fitted(.)| genotype*treatment, id = 0.05, adj = -0.3 )
In contrast, the best model (VarIDent) gave me terrible graphs:
p<-ggplot(my.data, aes(x =TOI , y = predict(VarIdent.attempt), colour =genxtrat)) + geom_point(size=2) p + geom_abline(intercept = 0, slope =1) + ggtitle("Predicted vs observed")
Including the graphs for groups:
plot(VarIdent.attempt, TOI ~ fitted(.)| genotype*treatment, id = 0.05, adj = -0.3 )
Including of course a very low correlation:
Correlation = 0.5726454
How can it be that the model that is found to be as far superior by both graphical (box-plots and QQ-plots) using standardized residuals, and analytical (Shapiro test, AIC, BIC) tools, yields an (observed vs predicted) values graph that is much worse?
Why does the group ATG5-TuMV show up as a straight horizontal line in Figure 7?
What conclusion should I draw from this analysis and what should I do?
I have tried with VarPower and VarExp and not only these models are worse than VarIdent (AIC, BIC), but also the correlation between observed and predicted did not improve. Modeling VarIdent with
gls gave me identical results.