[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