[R] decide between polynomial vs ordered factor model (lme)

Dimitris Rizopoulos dimitris.rizopoulos at med.kuleuven.be
Mon Jan 9 15:28:53 CET 2006


I think that these models are not nested and thus the LRT produced by 
anova.lme() will not be valid; AIC and BIC could be more relevant. In 
terms of interpretability, I'd say that a model treating 'zeitn' as a 
factor is much easier to explain than a model with 4th order 
polynomial.

I hope it helps.

Best,
Dimitris

----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/(0)16/336899
Fax: +32/(0)16/337015
Web: http://www.med.kuleuven.be/biostat/
     http://www.student.kuleuven.be/~m0390867/dimitris.htm



----- Original Message ----- 
From: "Leo Gürtler" <leog at anicca-vijja.de>
To: <r-help at stat.math.ethz.ch>
Sent: Monday, January 09, 2006 2:59 PM
Subject: [R] decide between polynomial vs ordered factor model (lme)


> Dear alltogether,
>
> two lme's, the data are available at:
>
> http://www.anicca-vijja.de/lg/hlm3_nachw.Rdata
>
> explanations of the data:
>
> nachw = post hox knowledge tests over 6 measure time points (= 
> equally
> spaced)
> zeitn = time points (n = 6)
> subgr = small learning groups (n = 28)
> gru = 4 different groups = treatment factor
>
> levels: time (=zeitn) (n=6) within subject (n=4) within smallgroups
> (=gru) (n = 28), i.e. n = 4 * 28 = 112 persons and 112 * 6 = 672 
> data points
>
> library(nlme)
> fitlme7 <- lme(nachw ~ I(zeitn-3.5) + I((zeitn-3.5)^2) +
> I((zeitn-3.5)^3) + I((zeitn-3.5)^4)*gru, random = list(subgr = ~ 1,
> subject = ~ zeitn), data = hlm3)
>
> fit5 <- lme(nachw ~ ordered(I(zeitn-3.5))*gru, random = list(subgr =
> ~ 1, subject = ~ zeitn), data = hlm3)
>
> anova( update(fit5, method="ML"), update(fitlme7, method="ML") )
>
> > anova( update(fit5, method="ML"), update(fitlme7, method="ML") )
>                                Model df      AIC      BIC    logLik 
> Test
> update(fit5, method = "ML")        1 29 2535.821 2666.619 -1238.911
> update(fitlme7, method = "ML")     2 16 2529.719 2601.883 -1248.860 
> 1 vs 2
>                                 L.Ratio p-value
> update(fit5, method = "ML")
> update(fitlme7, method = "ML") 19.89766  0.0978
> >
>
> shows that both are ~ equal, although I know about the uncertainty 
> of ML
> tests with lme(). Both models show that the ^2 and the ^4 terms are
> important parts of the model.
>
> My question is:
>
> - Is it legitim to choose a model based on these outputs according 
> to
> theoretical considerations instead of statistical tests that not 
> really
> show a superiority of one model over the other one?
>
> - Is there another criterium I've overseen to decide which model can 
> be
> clearly prefered?
>
> - The idea behind that is that in the one model (fit5) the second
> contrast of the factor (gru) is statistically significant, although 
> not
> the whole factor in the anova output.
> In the other model, this is not the case.
> Theoretically interesting is of course the significance of the 
> second
> contrast of gru, as it shows a tendency of one treatment being 
> slightly
> superior. I want to choose this model but I am not sure whether this 
> is
> proper action. Both models shows this trend, but only one model 
> clearly
> indicates that this trend bears some empirical meaning.
>
> Thanks for any suggestions,
>
> leo
>
>
> here are the outputs for each model:
>
>> fitlme7 <- lme(nachw ~ I(zeitn-3.5) + I((zeitn-3.5)^2) +
> I((zeitn-3.5)^3) + I((zeitn-3.5)^4)*gru, random = list(subgr = ~ 1,
> subject = ~ zeitn), data = hlm3)
>> plot(augPred(fitlme7), layout=c(14,8))
>> summary(fitlme7); anova(fitlme7); intervals(fitlme7)
> Linear mixed-effects model fit by REML
> Data: hlm3
>       AIC      BIC    logLik
>  2582.934 2654.834 -1275.467
>
> Random effects:
> Formula: ~1 | subgr
>        (Intercept)
> StdDev:   0.5833797
>
> Formula: ~zeitn | subject %in% subgr
> Structure: General positive-definite, Log-Cholesky parametrization
>            StdDev    Corr
> (Intercept) 0.6881908 (Intr)
> zeitn       0.1936087 -0.055
> Residual    1.3495785
>
> Fixed effects: nachw ~ I(zeitn - 3.5) + I((zeitn - 3.5)^2) + 
> I((zeitn -
> 3.5)^3) +      I((zeitn - 3.5)^4) * gru
>                            Value  Std.Error  DF   t-value p-value
> (Intercept)              4.528757 0.17749012 553 25.515542  0.0000
> I(zeitn - 3.5)           0.010602 0.08754449 553  0.121100  0.9037
> I((zeitn - 3.5)^2)       0.815693 0.09765075 553  8.353171  0.0000
> I((zeitn - 3.5)^3)       0.001336 0.01584169 553  0.084329  0.9328
> I((zeitn - 3.5)^4)      -0.089655 0.01405811 553 -6.377486  0.0000
> gru1                     0.187181 0.30805090  24  0.607630  0.5491
> gru2                     0.532665 0.30805090  24  1.729147  0.0966
> gru3                    -0.046305 0.30805090  24 -0.150317  0.8818
> I((zeitn - 3.5)^4):gru1 -0.007860 0.00600928 553 -1.307993  0.1914
> I((zeitn - 3.5)^4):gru2 -0.001259 0.00600928 553 -0.209516  0.8341
> I((zeitn - 3.5)^4):gru3 -0.000224 0.00600928 553 -0.037225  0.9703
> Correlation:
>                        (Intr) I(-3.5 I((-3.5)^2 I((-3.5)^3 
> I((z-3.5)^4)
> I(zeitn - 3.5)           0.071
> I((zeitn - 3.5)^2)      -0.465  0.000
> I((zeitn - 3.5)^3)       0.000 -0.914  0.000
> I((zeitn - 3.5)^4)       0.401  0.000 -0.977      0.000
> gru1                     0.000  0.000  0.000      0.000      0.000
> gru2                     0.000  0.000  0.000      0.000      0.000
> gru3                     0.000  0.000  0.000      0.000      0.000
> I((zeitn - 3.5)^4):gru1  0.000  0.000  0.000      0.000      0.000
> I((zeitn - 3.5)^4):gru2  0.000  0.000  0.000      0.000      0.000
> I((zeitn - 3.5)^4):gru3  0.000  0.000  0.000      0.000      0.000
>                        gru1   gru2   gru3   I((-3.5)^4):1 
> I((-3.5)^4):2
> I(zeitn - 3.5)
> I((zeitn - 3.5)^2)
> I((zeitn - 3.5)^3)
> I((zeitn - 3.5)^4)
> gru1
> gru2                     0.000
> gru3                     0.000  0.000
> I((zeitn - 3.5)^4):gru1 -0.287  0.000  0.000
> I((zeitn - 3.5)^4):gru2  0.000 -0.287  0.000  0.000
> I((zeitn - 3.5)^4):gru3  0.000  0.000 -0.287  0.000         0.000
>
> Standardized Within-Group Residuals:
>       Min         Q1        Med         Q3        Max
> -3.1326192 -0.5888543  0.0239228  0.6519002  2.1238820
>
> Number of Observations: 672
> Number of Groups:
>             subgr subject %in% subgr
>                28                112
>                       numDF denDF   F-value p-value
> (Intercept)                1   553 1426.5275  <.0001
> I(zeitn - 3.5)             1   553    0.2381  0.6258
> I((zeitn - 3.5)^2)         1   553   98.6712  <.0001
> I((zeitn - 3.5)^3)         1   553    0.0071  0.9328
> I((zeitn - 3.5)^4)         1   553   40.6723  <.0001
> gru                        3    24    1.0410  0.3924
> I((zeitn - 3.5)^4):gru     3   553    0.5854  0.6248
> Approximate 95% confidence intervals
>
> Fixed effects:
>                              lower          est.        upper
> (Intercept)              4.18011938  4.5287566579  4.877393940
> I(zeitn - 3.5)          -0.16135875  0.0106016498  0.182562052
> I((zeitn - 3.5)^2)       0.62388162  0.8156933820  1.007505144
> I((zeitn - 3.5)^3)      -0.02978133  0.0013359218  0.032453178
> I((zeitn - 3.5)^4)      -0.11726922 -0.0896553959 -0.062041570
> gru1                    -0.44860499  0.1871808283  0.822966643
> gru2                    -0.10312045  0.5326653686  1.168451183
> gru3                    -0.68209096 -0.0463051419  0.589480673
> I((zeitn - 3.5)^4):gru1 -0.01966389 -0.0078600880  0.003943709
> I((zeitn - 3.5)^4):gru2 -0.01306284 -0.0012590380  0.010544759
> I((zeitn - 3.5)^4):gru3 -0.01202749 -0.0002236923  0.011580105
> attr(,"label")
> [1] "Fixed effects:"
>
> Random Effects:
>  Level: subgr
>                    lower      est.     upper
> sd((Intercept)) 0.3459779 0.5833797 0.9836812
>  Level: subject
>                            lower        est.     upper
> sd((Intercept))         0.4388885  0.68819079 1.0791046
> sd(zeitn)               0.1320591  0.19360866 0.2838449
> cor((Intercept),zeitn) -0.4835884 -0.05541043 0.3941661
>
> Within-group standard error:
>   lower     est.    upper
> 1.267548 1.349579 1.436918
>
> #########################################################
> an the other model:
>
>> summary(fit5); anova(fit5); intervals(fit5)
> Linear mixed-effects model fit by REML
> Data: hlm3
>       AIC      BIC    logLik
>  2564.135 2693.878 -1253.067
>
> Random effects:
> Formula: ~1 | subgr
>        (Intercept)
> StdDev:   0.5833753
>
> Formula: ~zeitn | subject %in% subgr
> Structure: General positive-definite, Log-Cholesky parametrization
>            StdDev    Corr
> (Intercept) 0.6453960 (Intr)
> zeitn       0.1709843 0.13
> Residual    1.3497627
>
> Fixed effects: nachw ~ ordered(I(zeitn - 3.5)) + gru + 
> ordered(I(zeitn -
> 3.5)):gru
>                                   Value Std.Error  DF  t-value 
> p-value
> (Intercept)                     5.587313 0.1505852 540 37.10400 
> 0.0000
> ordered(I(zeitn - 3.5)).L       0.072572 0.1443422 540  0.50278 
> 0.6153
> ordered(I(zeitn - 3.5)).Q       1.266731 0.1275406 540  9.93198 
> 0.0000
> ordered(I(zeitn - 3.5)).C       0.010754 0.1275406 540  0.08432 
> 0.9328
> ordered(I(zeitn - 3.5))^4      -0.813277 0.1275406 540 -6.37662 
> 0.0000
> ordered(I(zeitn - 3.5))^5       0.070373 0.1275406 540  0.55177 
> 0.5813
> gru1                            0.056700 0.3011704  24  0.18826 
> 0.8523
> gru2                            0.679057 0.3011704  24  2.25473 
> 0.0335
> gru3                           -0.141425 0.3011704  24 -0.46958 
> 0.6429
> ordered(I(zeitn - 3.5)).L:gru1 -0.070352 0.2886844 540 -0.24370 
> 0.8076
> ordered(I(zeitn - 3.5)).Q:gru1 -0.360380 0.2550812 540 -1.41281 
> 0.1583
> ordered(I(zeitn - 3.5)).C:gru1 -0.162411 0.2550812 540 -0.63670 
> 0.5246
> ordered(I(zeitn - 3.5))^4:gru1  0.086343 0.2550812 540  0.33849 
> 0.7351
> ordered(I(zeitn - 3.5))^5:gru1 -0.017207 0.2550812 540 -0.06746 
> 0.9462
> ordered(I(zeitn - 3.5)).L:gru2  0.788896 0.2886844 540  2.73273 
> 0.0065
> ordered(I(zeitn - 3.5)).Q:gru2  0.033386 0.2550812 540  0.13089 
> 0.8959
> ordered(I(zeitn - 3.5)).C:gru2  0.089757 0.2550812 540  0.35188 
> 0.7251
> ordered(I(zeitn - 3.5))^4:gru2 -0.402616 0.2550812 540 -1.57839 
> 0.1151
> ordered(I(zeitn - 3.5))^5:gru2 -0.507855 0.2550812 540 -1.99095 
> 0.0470
> ordered(I(zeitn - 3.5)).L:gru3 -0.439200 0.2886844 540 -1.52138 
> 0.1287
> ordered(I(zeitn - 3.5)).Q:gru3  0.026105 0.2550812 540  0.10234 
> 0.9185
> ordered(I(zeitn - 3.5)).C:gru3 -0.273643 0.2550812 540 -1.07277 
> 0.2839
> ordered(I(zeitn - 3.5))^4:gru3 -0.163738 0.2550812 540 -0.64191 
> 0.5212
> ordered(I(zeitn - 3.5))^5:gru3  0.204174 0.2550812 540  0.80043 
> 0.4238
> Correlation:
>                               (Intr) or(I(-3.5)).L or(I(-3.5)).Q
> or(I(-3.5)).C or(I(-3.5))^4 or(I(-3.5))^5 gru1 gru2 gru3 
> o(I(-3.5)).L:1
> o(I(-3.5)).Q:1 o(I(-3.5)).C:1
> ordered(I(zeitn - 3.5)).L
> 0.2
>
>
> ordered(I(zeitn - 3.5)).Q      0.0
> 0.0
>
>
> ordered(I(zeitn - 3.5)).C      0.0    0.0
> 0.0
>
>
> ordered(I(zeitn - 3.5))^4      0.0    0.0           0.0
> 0.0
>
>
> ordered(I(zeitn - 3.5))^5      0.0    0.0           0.0
> 0.0
> 0.0
>
>
> gru1                           0.0    0.0           0.0
> 0.0           0.0
> 0.0
> gru2                           0.0    0.0           0.0
> 0.0           0.0           0.0
> 0.0
> gru3                           0.0    0.0           0.0
> 0.0           0.0           0.0           0.0
> 0.0
> ordered(I(zeitn - 3.5)).L:gru1 0.0    0.0           0.0
> 0.0           0.0           0.0           0.2  0.0
> 0.0
> ordered(I(zeitn - 3.5)).Q:gru1 0.0    0.0           0.0
> 0.0           0.0           0.0           0.0  0.0  0.0
> 0.0
> ordered(I(zeitn - 3.5)).C:gru1 0.0    0.0           0.0
> 0.0           0.0           0.0           0.0  0.0  0.0  0.0
> 0.0
> ordered(I(zeitn - 3.5))^4:gru1 0.0    0.0           0.0
> 0.0           0.0           0.0           0.0  0.0  0.0  0.0
> 0.0            0.0
> ordered(I(zeitn - 3.5))^5:gru1 0.0    0.0           0.0
> 0.0           0.0           0.0           0.0  0.0  0.0  0.0
> 0.0            0.0
> ordered(I(zeitn - 3.5)).L:gru2 0.0    0.0           0.0
> 0.0           0.0           0.0           0.0  0.2  0.0  0.0
> 0.0            0.0
> ordered(I(zeitn - 3.5)).Q:gru2 0.0    0.0           0.0
> 0.0           0.0           0.0           0.0  0.0  0.0  0.0
> 0.0            0.0
> ordered(I(zeitn - 3.5)).C:gru2 0.0    0.0           0.0
> 0.0           0.0           0.0           0.0  0.0  0.0  0.0
> 0.0            0.0
> ordered(I(zeitn - 3.5))^4:gru2 0.0    0.0           0.0
> 0.0           0.0           0.0           0.0  0.0  0.0  0.0
> 0.0            0.0
> ordered(I(zeitn - 3.5))^5:gru2 0.0    0.0           0.0
> 0.0           0.0           0.0           0.0  0.0  0.0  0.0
> 0.0            0.0
> ordered(I(zeitn - 3.5)).L:gru3 0.0    0.0           0.0
> 0.0           0.0           0.0           0.0  0.0  0.2  0.0
> 0.0            0.0
> ordered(I(zeitn - 3.5)).Q:gru3 0.0    0.0           0.0
> 0.0           0.0           0.0           0.0  0.0  0.0  0.0
> 0.0            0.0
> ordered(I(zeitn - 3.5)).C:gru3 0.0    0.0           0.0
> 0.0           0.0           0.0           0.0  0.0  0.0  0.0
> 0.0            0.0
> ordered(I(zeitn - 3.5))^4:gru3 0.0    0.0           0.0
> 0.0           0.0           0.0           0.0  0.0  0.0  0.0
> 0.0            0.0
> ordered(I(zeitn - 3.5))^5:gru3 0.0    0.0           0.0
> 0.0           0.0           0.0           0.0  0.0  0.0  0.0
> 0.0            0.0
>                               o(I(-3.5))^4:1 o(I(-3.5))^5:1
> o(I(-3.5)).L:2 o(I(-3.5)).Q:2 o(I(-3.5)).C:2 o(I(-3.5))^4:2
> o(I(-3.5))^5:2 o(I(-3.5)).L:3 o(I(-3.5)).Q:3
> ordered(I(zeitn -
> 3.5)).L
>
>
> ordered(I(zeitn -
> 3.5)).Q
>
>
> ordered(I(zeitn -
> 3.5)).C
>
>
> ordered(I(zeitn -
> 3.5))^4
>
>
> ordered(I(zeitn -
> 3.5))^5
>
>
> gru1
>
>
>
> gru2
>
>
>
> gru3
>
>
>
> ordered(I(zeitn -
> 3.5)).L:gru1
>
>
> ordered(I(zeitn -
> 3.5)).Q:gru1
>
>
> ordered(I(zeitn -
> 3.5)).C:gru1
>
>
> ordered(I(zeitn -
> 3.5))^4:gru1
>
>
> ordered(I(zeitn - 3.5))^5:gru1
> 0.0
>
>
> ordered(I(zeitn - 3.5)).L:gru2 0.0
> 0.0
>
>
> ordered(I(zeitn - 3.5)).Q:gru2 0.0            0.0
> 0.0
>
>
> ordered(I(zeitn - 3.5)).C:gru2 0.0            0.0
> 0.0
> 0.0
>
>
> ordered(I(zeitn - 3.5))^4:gru2 0.0            0.0
> 0.0            0.0
> 0.0
> ordered(I(zeitn - 3.5))^5:gru2 0.0            0.0
> 0.0            0.0            0.0
> 0.0
> ordered(I(zeitn - 3.5)).L:gru3 0.0            0.0
> 0.0            0.0            0.0            0.0
> 0.0
> ordered(I(zeitn - 3.5)).Q:gru3 0.0            0.0
> 0.0            0.0            0.0            0.0
> 0.0            0.0
> ordered(I(zeitn - 3.5)).C:gru3 0.0            0.0
> 0.0            0.0            0.0            0.0
> 0.0            0.0            0.0
> ordered(I(zeitn - 3.5))^4:gru3 0.0            0.0
> 0.0            0.0            0.0            0.0
> 0.0            0.0            0.0
> ordered(I(zeitn - 3.5))^5:gru3 0.0            0.0
> 0.0            0.0            0.0            0.0
> 0.0            0.0            0.0
>                               o(I(-3.5)).C:3 o(I(-3.5))^4:3
> ordered(I(zeitn - 3.5)).L
> ordered(I(zeitn - 3.5)).Q
> ordered(I(zeitn - 3.5)).C
> ordered(I(zeitn - 3.5))^4
> ordered(I(zeitn - 3.5))^5
> gru1
> gru2
> gru3
> ordered(I(zeitn - 3.5)).L:gru1
> ordered(I(zeitn - 3.5)).Q:gru1
> ordered(I(zeitn - 3.5)).C:gru1
> ordered(I(zeitn - 3.5))^4:gru1
> ordered(I(zeitn - 3.5))^5:gru1
> ordered(I(zeitn - 3.5)).L:gru2
> ordered(I(zeitn - 3.5)).Q:gru2
> ordered(I(zeitn - 3.5)).C:gru2
> ordered(I(zeitn - 3.5))^4:gru2
> ordered(I(zeitn - 3.5))^5:gru2
> ordered(I(zeitn - 3.5)).L:gru3
> ordered(I(zeitn - 3.5)).Q:gru3
> ordered(I(zeitn - 3.5)).C:gru3
> ordered(I(zeitn - 3.5))^4:gru3 0.0
> ordered(I(zeitn - 3.5))^5:gru3 0.0            0.0
>
> Standardized Within-Group Residuals:
>        Min          Q1         Med          Q3         Max
> -3.10206117 -0.62626454  0.02807962  0.64554138  2.13155536
>
> Number of Observations: 672
> Number of Groups:
>             subgr subject %in% subgr
>                28                112
>                            numDF denDF   F-value p-value
> (Intercept)                     1   540 1426.5315  <.0001
> ordered(I(zeitn - 3.5))         5   540   27.9740  <.0001
> gru                             3    24    1.0410  0.3924
> ordered(I(zeitn - 3.5)):gru    15   540    1.4115  0.1363
> Approximate 95% confidence intervals
>
> Fixed effects:
>                                    lower        est.        upper
> (Intercept)                     5.2915086  5.58731309  5.883117621
> ordered(I(zeitn - 3.5)).L      -0.2109689  0.07257212  0.356113124
> ordered(I(zeitn - 3.5)).Q       1.0161942  1.26673073  1.517267227
> ordered(I(zeitn - 3.5)).C      -0.2397825  0.01075396  0.261290456
> ordered(I(zeitn - 3.5))^4      -1.0638138 -0.81327731 -0.562740815
> ordered(I(zeitn - 3.5))^5      -0.1801634  0.07037312  0.320909612
> gru1                           -0.5648856  0.05669953  0.678284624
> gru2                            0.0574723  0.67905739  1.300642487
> gru3                           -0.7630097 -0.14142458  0.480160517
> ordered(I(zeitn - 3.5)).L:gru1 -0.6374343 -0.07035232  0.496729683
> ordered(I(zeitn - 3.5)).Q:gru1 -0.8614532 -0.36038020  0.140692783
> ordered(I(zeitn - 3.5)).C:gru1 -0.6634839 -0.16241093  0.338662057
> ordered(I(zeitn - 3.5))^4:gru1 -0.4147301  0.08634286  0.587415843
> ordered(I(zeitn - 3.5))^5:gru1 -0.5182803 -0.01720729  0.483865692
> ordered(I(zeitn - 3.5)).L:gru2  0.2218139  0.78889594  1.355977946
> ordered(I(zeitn - 3.5)).Q:gru2 -0.4676866  0.03338637  0.534459352
> ordered(I(zeitn - 3.5)).C:gru2 -0.4113159  0.08975711  0.590830099
> ordered(I(zeitn - 3.5))^4:gru2 -0.9036894 -0.40261640  0.098456584
> ordered(I(zeitn - 3.5))^5:gru2 -1.0089275 -0.50785453 -0.006781542
> ordered(I(zeitn - 3.5)).L:gru3 -1.0062815 -0.43919953  0.127882479
> ordered(I(zeitn - 3.5)).Q:gru3 -0.4749680  0.02610502  0.527178001
> ordered(I(zeitn - 3.5)).C:gru3 -0.7747163 -0.27364336  0.227429629
> ordered(I(zeitn - 3.5))^4:gru3 -0.6648114 -0.16373838  0.337334604
> ordered(I(zeitn - 3.5))^5:gru3 -0.2968991  0.20417390  0.705246883
> attr(,"label")
> [1] "Fixed effects:"
>
> Random Effects:
>  Level: subgr
>                    lower      est.     upper
> sd((Intercept)) 0.3464888 0.5833753 0.9822158
>  Level: subject
>                            lower      est.     upper
> sd((Intercept))         0.3640439 0.6453960 1.1441916
> sd(zeitn)               0.1000264 0.1709843 0.2922790
> cor((Intercept),zeitn) -0.6712236 0.1295558 0.7907922
>
> Within-group standard error:
>   lower     est.    upper
> 1.265702 1.349763 1.439406
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html
> 


Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm




More information about the R-help mailing list