0
$\begingroup$

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

Graph of Standardized Residuals for each genotype and treatment

qqnorm(my.1rst.attempt, abline = c(0,1)) #QQ plot:

QQ plot for my first model

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:

Standardized Residuals for each genotype and treatment after VarIdent

QQ plot for my VarIdent model

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

Predicted vs observed for my.first.attempt model

cor(predict(my.first.attempt),my.data$TOI)

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 )

Predicted vs observed for my.first.attempt model for the groups

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

Predicted vs observed for my VarIdent model

Including the graphs for groups:

plot(VarIdent.attempt, TOI ~ fitted(.)| genotype*treatment, id = 0.05, adj = -0.3 )

Predicted vs observed for my VarIdent model for the groups

Including of course a very low correlation:

cor(predict(VarIdent.attempt),my.data$TOI)

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.