Graphical validation of residuals of two models or their AICc?
Dear Andrew, Thank you for the thoughtful comments. I'm mainly trying to estimate parameters with uncertainty. The model is a multivariate multilevel model where `cal` defines the dependent variables. It does seem that the `cal` variables (there are 7 of them) don't correlate much with each other for each subject (`id`) across all subjects. As a result, I used `id = pdDiag(~cal +0)` to honor this feature. The problem I'm facing is that the `cal` variables have very different descriptive statistics (e.g., their range varies substantially). As a result, I did expect the residuals of the model to be wild. To partly mitigate this, I let the residuals be spread differently across `cal:time` combinations for each subject across all subjects, that is: `weights = varIdent(form = ~ cal:time)`, correction: it must have been: `varIdent(form = ~ 1|cal*time)` . On the issue of residuals, there is one thought that comes my mind: For a multivariate multilevel model like mine, is it even fair to look at residuals like univariate multilevel models, that is: simply `plot(fitted(m2), resid(m2, type="normalized"))` OR rather we should break the residuals down by `cal` levels across `id` levels and then look at the residuals' pattern? My suggestion in the previous paragraph is based on the following post from our archives: ( https://stat.ethz.ch/pipermail/r-sig-mixed-models/2020q2/028684.html) Thanks, Tim M
On Sun, Jan 15, 2023 at 2:32 PM Andrew Robinson <apro at unimelb.edu.au> wrote:
Neither residual plot really floats my boat. Presumably you have reasons to compare these model specifications that are unrelated to the data as we see them. Ultimately it depends on the use to which the model will be put. Are you trying to estimate parameters with uncertainty or to make a prediction of some future state or something else? I'd note that AIC is computed from LL at its maximum conditional on the model. If the model is wrong then the MLE is unknowably wrong and the AIC is unknowably wrong. So absent other context I focus on the information that is available about the model, namely the residual plots. I'd add that making time a factor is a curious decision. Cheers, Andrew -- Andrew Robinson Chief Executive Officer, CEBRA and Professor of Biosecurity, School/s of BioSciences and Mathematics & Statistics University of Melbourne, VIC 3010 Australia Tel: (+61) 0403 138 955 Email: apro at unimelb.edu.au Website: https://researchers.ms.unimelb.edu.au/~apro at unimelb/ I acknowledge the Traditional Owners of the land I inhabit, and pay my respects to their Elders. On 15 Jan 2023 at 6:45 AM +1100, Timothy MacKenzie <fswfswt at gmail.com>, wrote: Dear All, I'm deciding between two candidate models (m1 and m2). AICc shows that m2 is to be preferred over m1. However, when I plot the residuals of m1 and m2, residuals of m1 look more reasonable than those of m2. I wonder which criteria to rely on to choose my model: graphical validation of residuals of the two models or their AICc? Thank you, Tim M #-----------------------------------------# dat <- read.csv("https://raw.githubusercontent.com/fpqq/w/main/c.csv") dat <- transform(dat, time = factor(time)) library(nlme) m1 <- lme(score ~ time*cal, random = list(id = pdDiag(~cal +0)), data = dat, weights = varIdent(form = ~ cal:time), control = lmeControl(msMaxIter = 2e2)) m2 <- lme(score ~ time*cal, random = list(id = pdDiag(~cal:time +0)), data = dat, weights = varIdent(form = ~ cal:time), control = lmeControl(msMaxIter = 2e2)) library(bbmle) AICctab(m1, m2) plot(m1) plot(m2) #-----------------------------------------# [[alternative HTML version deleted]]
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models