[R-sig-ME] gls() results for unbalanced data with corStruct correlation class change if Initialize() is called first
Ben Bolker
bbo|ker @end|ng |rom gm@||@com
Thu Jul 16 18:15:42 CEST 2020
My best guess is that you've found a bug in gls. Thanks for the
detailed test case! I did test in lme() {as you may already know, the
compound-symmetric GLS is equivalent, *if the estimated correlation is
positive*, to a random-intercepts model}, and found that the lme()
results agreed closely with the "direct call" results (standard errors
differ slightly, but that's to be expected - they're much more sensitive
to any numerical instability). I put my explorations up at:
https://gist.github.com/bbolker/501c1dc6b1ed8be5c46ef66da5a964aa
(sorry about tidyverse code, I find myself getting sucked in thiese days
...)
On 7/16/20 11:21 AM, Mychaleckyj, Josyf C (jcm6t) wrote:
> I’ve encountered a strange situation where the results of a gls model with a corStruct correlation structure are affected by whether Initialize() is explicitly called before the gls() model is executed.
>
> This seems to be restricted to unbalanced data sets. I have pasted a test case using the Orthodont data set and corCompSymm below, although I first noticed it while running a custom corStruct class so it would appear to be a more general problem than this test case as you might expect. I have not tested this in lme().
>
> I’m aware that calling Initialize() outside of gls or lme is generally discouraged but was surprised that it changed the analysis results. You could imagine calling it in an interactive modeling session to view the correlation structure without knowing that you might have changed future results. I would have expected the results to be identical.
>
> I’d also like to know - which is the “correct” result or are they considered practically equivalent, although the differences seem quite large for the latter.
>
> Just to the clear, this is not about the relative merits of balanced vs unbalanced data. Its about the veracity of the results.
>
> Thanks,
> Joe.
>
>
>
> R version 3.6.0 (Linux CentOS 7)
>
>> library(nlme)
>>
>> packageVersion('nlme')
> [1] ‘3.1.147’
>> #—————————————————————
>> # the values are fixed for all models
>> rho <- 0.3
>> fixform <- distance ~ age + factor(Sex)
>>
>> # remove the grouping and other classes to keep it simple
>> # full balanced data set
>> Obal<-as.data.frame(Orthodont)
>> class(Obal)
> [1] "data.frame"
>> # create unbalanced Orthodont data set with no singletons
>> Ounbal <-as.data.frame(Orthodont[c(3,4,5, 7,8, 11, 12, 15, 16, 19, 20:28,31,32,35,36,39,
> + 40,41,43,44,47,48,51,52,53,55,56,59,60,63,64,67,68,71,72,75,76,
> + 79,80,83:85,87,88,91,92,95,96,99:101,103:105,107,108),])
>> # sparser unbalanced data set with singletons
>> Ounbal2<-as.data.frame(Orthodont[c(1,5,9,13,17,18,21:23,30:31,44,45:47,56:58,65:71,74:76,77,81,85,89,90,
> + 93:94, 97,101:108),])
>> Subjs<-c(paste0(c("M"),sprintf("%02d",1:16)), paste0(c("F"),sprintf("%02d",1:11)))
>>
>> table(Obal$Subject)[Subjs]
> M01 M02 M03 M04 M05 M06 M07 M08 M09 M10 M11 M12 M13 M14 M15 M16 F01 F02 F03 F04
> 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
> F05 F06 F07 F08 F09 F10 F11
> 4 4 4 4 4 4 4
>> table(Ounbal$Subject)[Subjs]
> M01 M02 M03 M04 M05 M06 M07 M08 M09 M10 M11 M12 M13 M14 M15 M16 F01 F02 F03 F04
> 2 3 2 2 2 4 4 2 2 2 3 2 2 3 2 2 2 2 2 2
> F05 F06 F07 F08 F09 F10 F11
> 2 3 2 2 2 3 3
>> table(Ounbal2$Subject)[Subjs]
> M01 M02 M03 M04 M05 M06 <NA> M08 <NA> <NA> M11 M12 <NA> M14 M15 <NA>
> 1 1 1 1 2 3 2 1 3 1 2
> F01 F02 F03 F04 F05 F06 F07 F08 F09 F10 F11
> 4 3 3 1 1 1 2 2 1 4 4
>> #----
>> # test with the full balanced data set
>> # pre-calling Initialize doesn't seem to make a difference
>> csB <- corCompSymm(value = rho, form = ~ 1 | Subject)
>> csB<-Initialize(csB, data=Obal)
>>
>> print(mod<-gls(fixform, correlation=csB, data=Obal))
> Generalized least squares fit by REML
> Model: fixform
> Data: Obal
> Log-restricted-likelihood: -218.7563
>
> Coefficients:
> (Intercept) age SexFemale
> 17.7067130 0.6601852 -2.3210227
>
> Correlation Structure: Compound symmetry
> Formula: ~1 | Subject
> Parameter estimate(s):
> Rho
> 0.6144908
> Degrees of freedom: 108 total; 105 residual
> Residual standard error: 2.305696
>> # compare with direct call
>> print(mod2<-gls(fixform, correlation=corCompSymm(value = rho, form = ~ 1 | Subject), data=Obal) )
> Generalized least squares fit by REML
> Model: fixform
> Data: Obal
> Log-restricted-likelihood: -218.7563
>
> Coefficients:
> (Intercept) age SexFemale
> 17.7067130 0.6601852 -2.3210227
>
> Correlation Structure: Compound symmetry
> Formula: ~1 | Subject
> Parameter estimate(s):
> Rho
> 0.6144908
> Degrees of freedom: 108 total; 105 residual
> Residual standard error: 2.305696
>> #--------------
>> # compare Unbalanced:
>> csU <- corCompSymm(value = rho, form = ~ 1 | Subject)
>> csU<-Initialize(csU, data=Ounbal)
>>
>> print(mod3<-gls(fixform, correlation=csU, data=Ounbal) )
> Generalized least squares fit by REML
> Model: fixform
> Data: Ounbal
> Log-restricted-likelihood: -130.7703
>
> Coefficients:
> (Intercept) age SexFemale
> 18.046228 0.659606 -3.262389
>
> Correlation Structure: Compound symmetry
> Formula: ~1 | Subject
> Parameter estimate(s):
> Rho
> 0.7211814
> Degrees of freedom: 64 total; 61 residual
> Residual standard error: 2.378517
>> # compare with direct call
>> print(mod4<-gls(fixform, correlation=corCompSymm(value = rho, form = ~ 1 | Subject), data=Ounbal) )
> Generalized least squares fit by REML
> Model: fixform
> Data: Ounbal
> Log-restricted-likelihood: -127.616
>
> Coefficients:
> (Intercept) age SexFemale
> 18.330637 0.634451 -2.947638
>
> Correlation Structure: Compound symmetry
> Formula: ~1 | Subject
> Parameter estimate(s):
> Rho
> 0.7598463
> Degrees of freedom: 64 total; 61 residual
> Residual standard error: 2.34705
>> #----
>> # compare with Unbalanced 2 with singletons
>> csU2 <- corCompSymm(value = rho, form = ~ 1 | Subject)
>> csU2<-Initialize(csU2, data=Ounbal2)
>>
>> print(mod5<-gls(fixform, correlation=csU2, data=Ounbal2))
> Generalized least squares fit by REML
> Model: fixform
> Data: Ounbal2
> Log-restricted-likelihood: -74.02154
>
> Coefficients:
> (Intercept) age SexFemale
> 18.5018545 0.5265219 -1.4587968
>
> Correlation Structure: Compound symmetry
> Formula: ~1 | Subject
> Parameter estimate(s):
> Rho
> 0.8736732
> Degrees of freedom: 44 total; 41 residual
> Residual standard error: 1.94714
>> # compare with direct call
>> print(mod6<-gls(fixform, correlation=corCompSymm(value = rho, form = ~ 1 | Subject), data=Ounbal2))
> Generalized least squares fit by REML
> Model: fixform
> Data: Ounbal2
> Log-restricted-likelihood: -77.29477
>
> Coefficients:
> (Intercept) age SexFemale
> 18.4979634 0.5550817 -1.6893916
>
> Correlation Structure: Compound symmetry
> Formula: ~1 | Subject
> Parameter estimate(s):
> Rho
> 0.8516777
> Degrees of freedom: 44 total; 41 residual
> Residual standard error: 2.024012
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list