[R-sig-ME] Graphical validation of residuals of two models or their AICc?

[Timothy MacKenzie]

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 using 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 using unimelb.edu.au
> Website: https://researchers.ms.unimelb.edu.au/~apro@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 using 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 using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>

	[[alternative HTML version deleted]]