[R-meta] AICc or variance components, which one matters more?

[Viechtbauer, Wolfgang (SP)]

Please see my comments below.

Best,
Wolfgang

>-----Original Message-----
>From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On
>Behalf Of Philippe Tadger
>Sent: Saturday, 13 November, 2021 11:03
>To: Luke Martinez; R meta
>Subject: Re: [R-meta] AICc or variance components, which one matters more?
>
>Dear Luke
>
>I will suggest to choose the model that can provide a reliable
>estimation of the variances. Usually a variance estimation equal to zero
>or a correlation estimation  equal to 1 or -1 are potential warning
>signal, your model is in the limit of the parameter space: so you can
>having an underestimation and/or a non-identifiable model. How can you
>be sure? Using the profile() function in each model. In your case model
>g1 is identifiable, but g2 is not identifiable. I'll choose g1.

This is not correct. Both variance componments are identifiable in both models, as the profile plots indicate. The estimate of the lab variance component in g2 just happens to be 0 (the peak of the likelihood profile is at 0).

>On 13/11/2021 05:24, Luke Martinez wrote:
>> Hello Colleagues,
>>
>> I've fit two candidate models (g1 and g2).
>>
>> With g1, I get two non-zero variance components for 'lab' and 'study'.
>>
>> With g2, I get a zero variance component for 'lab' and a non-zero one
>> for 'study'.
>>
>> So, from the perspective of variance components, g1 seems like a better model.

I don't know why you think a model is better if the estimates of its variance components are non-zero.

>> But when I compare the models using AICc, g2 seems like a better model.

Personally, I would rather think about the meaning of these models and use the one that makes more sense in the given context.

>> I wonder, then, which criterion should I use to choose the model (AICc
>> or the variance components)?
>>
>> Thanks,
>> Luke
>> For reproducibility, I'm showing a somewhat similar situation with a
>> small section of my data below.
>>
>> m="
>>        lab study   yi         vi es_id
>>           1     1 1.04 0.48503768 1
>>           1     1 0.96 0.51076604 2
>>           1     2 1.71 0.05767389 3
>>           2     2 1.52 0.07539841 4
>>           1     3 1.91 0.31349510 5
>>           2     4 3.01 0.67910095 6
>>           2     4 3.62 0.50670360 7
>>           9     5 0.99 0.18297170 8
>>           9     5 0.43 0.37225851 9
>>           9     6 1.68 0.39072390 10
>>           9     6 1.25 0.02879550 11"
>> dd <- read.table(text=m,h=T)
>>
>> (g1=rma.mv(yi, vi, random = ~1|lab/study, data = dd))
>>                      estim    sqrt  nlvls  fixed     factor
>> sigma^2.1  0.1245  0.3529      3     no        lab
>> sigma^2.2  0.2535  0.5035      7     no  lab/study
>>
>> (g2=rma.mv(yi, vi, random = list(~1|lab, ~1|study), data = dd))
>>                       estim    sqrt  nlvls  fixed  factor
>> sigma^2.1  0.0000  0.0000      3     no     lab
>> sigma^2.2  0.5127  0.7160      6     no   study
>>
>> fitstats(g1,g2)[5,]
>>              g1          g2
>> AICc: 30.85992 29.73897
>> (for full data, AICc of g2 is more noticeably smaller than that of g1)