[R] lme syntax for P&B examples
Henric Nilsson
henric.nilsson at statisticon.se
Thu Feb 9 14:25:17 CET 2006
Paul Cossens said the following on 2006-02-09 06:21:
> Hi Harold,
>
>
> Thanks for your reply. I had already looked at all the reading material
> you suggested but updated to the latest Matrix
> as recommneded then spent all day trying to figure out what is
> happening.
>
> I worked through the problems and give my workings below that others may
> find useful.
> (My notation is to use lme> to show lme commands and lmer> to show lmer
> commands.
> I worked on two sessions in parallel. My comments are preceded by double
> hashes '##' and
> questions '##??'. I haven't included the datasets.)
>
> I have a couple of comments and outstanding issues:
>
> 1. In the Pixel data set and formulas I think the formulas are printed
> incorrectly in the
> book as some use 'I(day^2)' while others use just 'day^2'. I have used
> 'I(day^2)'. I'm not sure why the I() function is used. In the fm4Pixel
> example below the answers don't match up exactly but are close.
>
> The lme example is
> fm1Pixel<-lme(pixel~day+I(day^2),data=Pixel,random = list(Dog=~day
> ,Side=~1))
> fm5Pixel <- update(fm1Pixel,pixel ~ day + I(day^2) + Side)
> which I have converted to lmer:
> fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side +(day|Dog), data = Pixel)
>
> The t-values for Side are close (sse below) but different enough to
> wonder if I am still doing something wrong?
It's wrong. Try
fm6Pixel <- lmer(pixel ~ day + I(day^2) + Side + (day|Dog) +
(1|Side:Dog), data = Pixel)
> 2. To me the specification description in the R-News article is
> confusing as it seems
> to suggest that nesting does not need to be completely specified if the
> groupings and nestings are clear in data set.
>
> Prof Bates article in R news vol 5/1 P 30 states "It happens in this
> case that the grouping factors 'id' and 'sch' are not nested but if they
> were nested there would be no change in the model specification"
>
> If the lme formula is
> fm1Oxide<-lme(Thickness~1,Oxide)
>
> I have found the formula lmer parlance should be:
> 'fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide)'
> not 'fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Wafer),data=Oxide)'
> as the article reads to me.
>
> In other words you always need to explicitly specify nesting levels.
>
> 3. I still can't figue out how to replicate the lme formula
> fm2Oxide<-update(fm1Oxide,random =~1|Lot)
> i.e
> formula(fm2Oxide)
> Thickness ~ 1
>
> If I simply drop the Lot:Wafer term as in 'fm2Oxide<-lmer(Thickness~
> (1|Lot),data=Oxide)'
> I get the error
>
> 'Error in x[[3]] : object is not subsettable'
>
> what's the solution?
fm3Oxide <- lmer(Thickness ~ 1 + (1|Lot), data = Oxide)
This seems to be a bug in `lmer', since one doesn't have to specifiy the
intercept explicitly in e.g.
lmer(Thickness ~ (1|Lot) + (1|Wafer), data = Oxide)
HTH,
Henric
>
> I'd be interested to read you article for further insights.
>
> Thanks
>
> Paul
>
>
>
> #############################################################
> #Oxide
> # P&B(2000) p167-170
>
> #NLME lme example
>
> lme>data(Oxide)
> lme>formula(Oxide)
> lme>Thickness ~ 1 | Lot/Wafer
> lme>fm1Oxide<-lme(Thickness~1,Oxide)
> lme> fm1Oxide
> Linear mixed-effects model fit by REML
> Data: Oxide
> Log-restricted-likelihood: -227.0110
> Fixed: Thickness ~ 1
> (Intercept)
> 2000.153
>
> Random effects:
> Formula: ~1 | Lot
> (Intercept)
> StdDev: 11.39768
>
> Formula: ~1 | Wafer %in% Lot
> (Intercept) Residual
> StdDev: 5.988802 3.545341
>
> Number of Observations: 72
> Number of Groups:
> Lot Wafer %in% Lot
> 8 24
>
> lme> intervals(fm1Oxide, which = "var-cov")
> Approximate 95% confidence intervals
>
> Random Effects:
> Level: Lot
> lower est. upper
> sd((Intercept)) 6.38881 11.39768 20.33355
> Level: Wafer
> lower est. upper
> sd((Intercept)) 4.063919 5.988802 8.825408
>
> Within-group standard error:
> lower est. upper
> 2.902491 3.545341 4.330572
> fm2Oxide<-update(fm1Oxide,random =~1|Lot)
> lme> fm2Oxide
> Linear mixed-effects model fit by REML
> Data: Oxide
> Log-restricted-likelihood: -245.5658
> Fixed: Thickness ~ 1
> (Intercept)
> 2000.153
>
> Random effects:
> Formula: ~1 | Lot
> (Intercept) Residual
> StdDev: 11.78447 6.282416
>
> Number of Observations: 72
> Number of Groups: 8
>
> lme>intervals(fm2Oxide, which = "var-cov")
> Approximate 95% confidence intervals
>
> Random Effects:
> Level: Lot
> lower est. upper
> sd((Intercept)) 6.864617 11.78447 20.23035
>
> Within-group standard error:
> lower est. upper
> 5.283116 6.282416 7.470733
>
> lme> coef(fm1Oxide, level=1)
> (Intercept)
> 1 1996.689
> 2 1988.931
> 3 2001.022
> 4 1995.682
> 5 2013.616
> 6 2019.561
> 7 1991.954
> 8 1993.767
> coef(fm1Oxide, level=1)
> (Intercept)
> 1 1996.689
> 2 1988.931
> 3 2001.022
> 4 1995.682
> 5 2013.616
> 6 2019.561
> 7 1991.954
> 8 1993.767
>> coef(fm1Oxide, level=2)
> (Intercept)
> 1/1 2003.235
> 1/2 1984.730
> 1/3 2001.146
> 2/1 1989.590
> 2/2 1988.097
> 2/3 1986.008
> 3/1 2002.495
> 3/2 2000.405
> 3/3 2000.405
> 4/1 1995.668
> 4/2 1998.951
> 4/3 1991.191
> 5/1 2009.184
> 5/2 2016.646
> 5/3 2018.735
> 6/1 2031.296
> 6/2 2021.745
> 6/3 2011.000
> 7/1 1990.204
> 7/2 1991.398
> 7/3 1991.995
> 8/1 1993.677
> 8/2 1995.170
> 8/3 1990.693
>
>
> ######## the is the lmer example using Matrix 0.995-5
>
>
> lmer>Oxide<-read.csv("Oxide.csv",header=TRUE);
> lmer>Oxide$Source<-as.factor(Oxide$Source)
> lmer>Oxide$Lot<-as.factor(Oxide$Lot)
> lmer>Oxide$Wafer<-as.factor(Oxide$Wafer)
> Lmer>Oxide$Site<-as.factor(Oxide$Site)
>
>
> ## Bates article in R news vol5/1 says specifying nesting explicity
> isn't necessary:
> ## P 30 "It happens in this case that the grouping factors 'id' and
> 'sch' are not
> ## nested but if they were nested there would be no change in the model
> specification"
> ## Following this one would expect that the following statement would
> automatically
> ## detect nesting
>
> lmer>fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Wafer),data=Oxide)
>
>
> lmer> fm1Oxide
> Linear mixed-effects model fit by REML
> Formula: Thickness ~ (1 | Lot) + (1 | Wafer)
> Data: Oxide
> AIC BIC logLik MLdeviance REMLdeviance
> 496.6093 503.4393 -245.3046 495.3528 490.6093
> Random effects:
> Groups Name Variance Std.Dev.
> Lot (Intercept) 138.9981 11.7897
> Wafer (Intercept) 1.4930 1.2219
> Residual 38.3490 6.1927
> # of obs: 72, groups: Lot, 8; Wafer, 3
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 2000.1528 4.2901 466.22
>
> ## The lme vs lmer std.devs for Lot are 11.39768 : 11.7987
> ## The lme vs lmer std.devs for Wafer are 5.988802 : 1.2219
> ## If my lmer specifcation is correct then the Wafer std.dev seems too
> big.
> ## also the levels aren't specifed correctly as the following ranef()
> function
> ## shows
>
> lmer> ranef(fm1Oxide)
> An object of class "lmer.ranef"
> [[1]]
> (Intercept)
> 1 -3.7058415
> 2 -12.0069264
> 3 0.9298293
> 4 -4.7839045
> 5 14.4056166
> 6 20.7661881
> 7 -8.7727375
> 8 -6.8322241
>
> [[2]]
> (Intercept)
> 1 0.9526407
> 2 -0.2750582
> 3 -0.6775825
>
> ###There is no nesting of wafers within lots
> ###Note however that the following appears to work the same as in the
> P&B text
>
> lmer> fm3Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide)
> lmer> fm3Oxide
> Linear mixed-effects model fit by REML
> Formula: Thickness ~ (1 | Lot) + (1 | Lot:Wafer)
> Data: Oxide
> AIC BIC logLik MLdeviance REMLdeviance
> 460.0221 466.8521 -227.0110 458.7378 454.0221
> Random effects:
> Groups Name Variance Std.Dev.
> Lot:Wafer (Intercept) 35.866 5.9888
> Lot (Intercept) 129.853 11.3953
> Residual 12.570 3.5454
> # of obs: 72, groups: Lot:Wafer, 24; Lot, 8
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 2000.1528 4.2309 472.75
>
> ## the following gives the coefficients for levels however the number of
> levels is
> ## in reverse order to lme
> ## Also note that the order of levels seems to have changed from the
> lmer model
> ## fm1Oxide where Lot coefficnets are in [[1]] and Wafer [[2]]
> lmer>ranef(fm3Oxide)
> An object of class "lmer.ranef"
> [[1]]
> (Intercept)
> 1:1 6.54582998
> 1:2 -11.95898390
> 1:3 4.45657680
> 2:1 0.65819690
> 2:2 -0.83412680
> 2:3 -2.92337998
> 3:1 1.47284009
> 3:2 -0.61641309
> 3:3 -0.61641309
> 4:1 -0.01366505
> 4:2 3.26944709
> 4:3 -4.49063615
> 5:1 -4.43133801
> 5:2 3.03028049
> 5:3 5.11953367
> 6:1 11.73559478
> 6:2 2.18472310
> 6:3 -8.56000754
> 7:1 -1.74970878
> 7:2 -0.55584982
> 7:3 0.04107966
> 8:1 -0.09041889
> 8:2 1.40190481
> 8:3 -3.07506629
>
> [[2]]
> (Intercept)
> 1 -3.463334
> 2 -11.221203
> 3 0.868982
> 4 -4.470850
> 5 13.462925
> 6 19.407266
> 7 -8.198657
> 8 -6.385129
>
> ##?? given that the fm3Oxide works one would expect that dropping the
> ##?? (1|Lot:Wafer) term would work however it give the following error
>
> lmer>fm2Oxide<-lmer(Thickness~ (1|Lot),data=Oxide)
> Error in x[[3]] : object is not subsettable
>
> ##?? the question is how to specify the lme model
> fm2Oxide<-update(fm1Oxide,random =~1|Lot)
> ##?? in lme syntax
>
>
> ######################################################################
>
> #Pixel
> # P&B(2000) p40-45
>
> ## the lme example is as follows
>
> lme>data(Pixel)
>
> ## firstly the book may have an error on P 42 line 1
> lme_wrong> fm1Pixel<-lme(pixel~day+day^2,data=Pixel,random = list(Dog= ~
> day ,Side= ~ 1))
> ###which gives:
> lme_wrong> intervals(fm1Pixel)
> Approximate 95% confidence intervals
>
> Fixed effects:
> lower est. upper
> (Intercept) 1071.415167 1093.21538 1115.015590
> day -1.126053 -0.14867 0.828713
> attr(,"label")
> [1] "Fixed effects:"
>
> Random Effects:
> Level: Dog
> lower est. upper
> sd((Intercept)) 17.3485293 31.4936604 57.1720308
> sd(day) 0.3586459 1.0719672 3.2040342
> cor((Intercept),day) -0.9822311 -0.7863546 0.2294876
> Level: Side
> lower est. upper
> sd((Intercept)) 8.519543 15.08995 26.72755
>
> Within-group standard error:
> lower est. upper
> 12.35975 14.53391 17.09052
>
> ## the book should probably give I(day^2) instead of day^2 on p42 line 1
> where
> ## as this gives the answer in the book
> lme> fm1Pixel<-lme(pixel~day+I(day^2),data=Pixel,random = list(Dog=~day
> ,Side=~1))
> lme> intervals(fm1Pixel)
> Approximate 95% confidence intervals
>
> Fixed effects:
> lower est. upper
> (Intercept) 1053.0968388 1073.3391382 1093.5814376
> day 4.3796925 6.1295971 7.8795016
> I(day^2) -0.4349038 -0.3673503 -0.2997967
> attr(,"label")
> [1] "Fixed effects:"
>
> Random Effects:
> Level: Dog
> lower est. upper
> sd((Intercept)) 15.9284469 28.3699038 50.5291851
> sd(day) 1.0812133 1.8437505 3.1440751
> cor((Intercept),day) -0.8944371 -0.5547222 0.1909581
> Level: Side
> lower est. upper
> sd((Intercept)) 10.41726 16.82431 27.17195
>
> Within-group standard error:
> lower est. upper
> 7.634522 8.989606 10.585209
>
> lme>VarCorr(fm1Pixel)
> Variance StdDev Corr
> Dog = pdLogChol(day)
> (Intercept) 804.851443 28.369904 (Intr)
> day 3.399416 1.843750 -0.555
> Side = pdLogChol(1)
> (Intercept) 283.057248 16.824305
> Residual 80.813009 8.989606
>
> lme> summary(fm1Pixel)
> Linear mixed-effects model fit by REML
> Data: Pixel
> AIC BIC logLik
> 841.2102 861.9712 -412.6051
>
> Random effects:
> Formula: ~day | Dog
> Structure: General positive-definite, Log-Cholesky parametrization
> StdDev Corr
> (Intercept) 28.369904 (Intr)
> day 1.843750 -0.555
>
> Formula: ~1 | Side %in% Dog
> (Intercept) Residual
> StdDev: 16.82431 8.989606
>
> Fixed effects: pixel ~ day + I(day^2)
> Value Std.Error DF t-value p-value
> (Intercept) 1073.3391 10.171686 80 105.52225 0
> day 6.1296 0.879321 80 6.97083 0
> I(day^2) -0.3674 0.033945 80 -10.82179 0
> Correlation:
> (Intr) day
> day -0.517
> I(day^2) 0.186 -0.668
>
> Standardized Within-Group Residuals:
> Min Q1 Med Q3 Max
> -2.82905723 -0.44918107 0.02554930 0.55721629 2.75196509
>
> Number of Observations: 102
> Number of Groups:
> Dog Side %in% Dog
> 10 20
>
> ## Note thate the formula in the book is on P44 sect 1.5.1 summary table
> is
> ## Fixed effects: pixel ~ day + day^2
> ## compared to above:
> ## Fixed effects: pixel ~ day + I(day^2)
> ## but the anova on page 45 is the same
>
> lme>fm2Pixel <- update(fm1Pixel,random = ~day | Dog)
> lme>anova(fm1Pixel,fm2Pixel)
> Model df AIC BIC logLik Test L.Ratio p-value
> fm1Pixel 1 8 841.2102 861.9712 -412.6051
> fm2Pixel 2 7 884.5196 902.6854 -435.2598 1 vs 2 45.3094 <.0001
>
>
> lme> fm3Pixel <- update(fm1Pixel,random = ~1 | Dog/Side)
> lme > anova(fm1Pixel,fm3Pixel)
> Model df AIC BIC logLik Test L.Ratio p-value
> fm1Pixel 1 8 841.2102 861.9712 -412.6051
> fm3Pixel 2 6 876.8390 892.4098 -432.4195 1 vs 2 39.62885 <.0001
>
>
> ##again there following seems to be a mistake in the book
>
> lme> fm4Pixel <- update(fm1Pixel,pixel ~ day + day^2 + Side)
> lme> summary(fm4Pixel)
> Linear mixed-effects model fit by REML
> Data: Pixel
> AIC BIC logLik
> 895.071 915.832 -439.5355
>
> Random effects:
> Formula: ~day | Dog
> Structure: General positive-definite, Log-Cholesky parametrization
> StdDev Corr
> (Intercept) 31.520129 (Intr)
> day 1.073342 -0.786
>
> Formula: ~1 | Side %in% Dog
> (Intercept) Residual
> StdDev: 15.01697 14.50742
>
> Fixed effects: pixel ~ day + Side
> Value Std.Error DF t-value p-value
> (Intercept) 1097.5272 11.562590 81 94.92053 0.0000
> day -0.1496 0.491218 81 -0.30451 0.7615
> SideR -8.6098 7.379984 9 -1.16664 0.2733
> Correlation:
> (Intr) day
> day -0.633
> SideR -0.319 0.000
>
> Standardized Within-Group Residuals:
> Min Q1 Med Q3 Max
> -3.73906417 -0.38367706 0.04758941 0.39690056 2.23720545
>
> Number of Observations: 102
> Number of Groups:
> Dog Side %in% Dog
> 10 20
>
> ###the book should probably give I(day^2) instead of day^2
> ### as the fixed effect coefficient names in the summary table
> ## differ from the formula given
> ## the following seems a partial correction
> lme> fm5Pixel <- update(fm1Pixel,pixel ~ day + I(day^2) + Side)
> lme>summary(fm5Pixel)
>
> Linear mixed-effects model fit by REML
> Data: Pixel
> AIC BIC logLik
> 835.8546 859.1193 -408.9273
>
> Random effects:
> Formula: ~day | Dog
> Structure: General positive-definite, Log-Cholesky parametrization
> StdDev Corr
> (Intercept) 28.463605 (Intr)
> day 1.843823 -0.553
>
> Formula: ~1 | Side %in% Dog
> (Intercept) Residual
> StdDev: 16.5072 8.983614
>
> Fixed effects: pixel ~ day + I(day^2) + Side
> Value Std.Error DF t-value p-value
> (Intercept) 1077.9484 10.862705 80 99.23388 0.0000
> day 6.1296 0.879023 80 6.97323 0.0000
> I(day^2) -0.3674 0.033923 80 -10.82914 0.0000
> SideR -9.2175 7.626768 9 -1.20858 0.2576
> Correlation:
> (Intr) day I(d^2)
> day -0.484
> I(day^2) 0.174 -0.667
> SideR -0.351 0.000 0.000
>
> Standardized Within-Group Residuals:
> Min Q1 Med Q3 Max
> -2.80982455 -0.47133415 0.02610263 0.54115378 2.77470104
>
> Number of Observations: 102
> Number of Groups:
> Dog Side %in% Dog
> 10 20
>
>
> ## Note however that the Fixed effects Values above differ from the
> values in the book
>
> ####### the is the lmer example for Pixel using Matrix 0.995-5
>
> ##the lmer example for the Pixel data is exported and reimported as a
> csv file
> lmer>Pixel<-read.csv("Pixel.csv",header=TRUE);
> lmer>Pixel$Side<-as.factor(Pixel$Side)
> lmer>Pixel$Dog<-as.factor(Pixel$Dog)
> lmer>Pixel$DS<-with(Pixel,Dog:Side)[drop=TRUE]
>
> lmer>fm1Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog)+(1|Dog:Side),
> data = Pixel)
> lmer>fm1Pixel
> Linear mixed-effects model fit by REML
> Formula: pixel ~ day + I(day^2) + (day | Dog) + (1 | Dog:Side)
> Data: Pixel
> AIC BIC logLik MLdeviance REMLdeviance
> 839.2102 857.585 -412.6051 827.3298 825.2102
> Random effects:
> Groups Name Variance Std.Dev. Corr
> Dog:Side (Intercept) 283.0567 16.8243
> Dog (Intercept) 804.8460 28.3698
> day 3.3994 1.8438 -0.555
> Residual 80.8130 8.9896
> # of obs: 102, groups: Dog:Side, 20; Dog, 10
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 1073.339126 10.171658 105.523
> day 6.129599 0.879321 6.971
> I(day^2) -0.367350 0.033945 -10.822
>
> Correlation of Fixed Effects:
> (Intr) day
> day -0.517
> I(day^2) 0.186 -0.668
>
> ##
>
> lme> VarCorr(fm1Pixel)
> $"Dog:Side"
> 1 x 1 Matrix of class "dpoMatrix"
> (Intercept)
> (Intercept) 283.0567
>
> $Dog
> 2 x 2 Matrix of class "dpoMatrix"
> (Intercept) day
> (Intercept) 804.84605 -29.015418
> day -29.01542 3.399415
>
> attr(,"sc")
> [1] 8.989606
>
> lmer> fm2Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog), data = Pixel)
> lmer> fm2Pixel
> Linear mixed-effects model fit by REML
> Formula: pixel ~ day + I(day^2) + (day | Dog)
> Data: Pixel
> AIC BIC logLik MLdeviance REMLdeviance
> 882.5196 898.2694 -435.2598 873.5964 870.5196
> Random effects:
> Groups Name Variance Std.Dev. Corr
> Dog (Intercept) 892.7720 29.8793
> day 3.0772 1.7542 -0.488
> Residual 197.5767 14.0562
> # of obs: 102, groups: Dog, 10
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 1072.92725 10.44452 102.726
> day 6.08910 1.14699 5.309
> I(day^2) -0.35677 0.05221 -6.833
>
> Correlation of Fixed Effects:
> (Intr) day
> day -0.541
> I(day^2) 0.286 -0.799
>
> lmer> anova(fm1Pixel,fm2Pixel)
>
> Data: Pixel
> Models:
> fm2Pixel: pixel ~ day + I(day^2) + (day | Dog)
> fm1Pixel: pixel ~ day + I(day^2) + (day | Dog) + (1 | Dog:Side)
> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
> fm2Pixel 6 882.52 898.27 -435.26
> fm1Pixel 7 839.21 857.59 -412.61 45.309 1 1.682e-11 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> ## Some of the statistics here are slightly different from above,
> notably the Df
> ## but I guess the result is the same
>
> lmer> fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|Dog:Side), data =
> Pixel)
> lmer> anova(fm1Pixel,fm3Pixel)
>
> Data: Pixel
> Models:
> fm3Pixel: pixel ~ day + I(day^2) + (1 | Dog:Side)
> fm1Pixel: pixel ~ day + I(day^2) + (day | Dog) + (1 | Dog:Side)
> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
> fm3Pixel 4 877.88 888.38 -434.94
> fm1Pixel 7 839.21 857.59 -412.61 44.67 3 1.087e-09 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> ## Some of the statistics are slightly different again
>
> lmer>fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side +(day|Dog), data =
> Pixel)
> lmer>fm4Pixel
>
> Linear mixed-effects model fit by REML
> Formula: pixel ~ day + I(day^2) + Side + (day | Dog)
> Data: Pixel
> AIC BIC logLik MLdeviance REMLdeviance
> 876.8204 895.1952 -431.4102 869.6765 862.8204
> Random effects:
> Groups Name Variance Std.Dev. Corr
> Dog (Intercept) 896.1278 29.9354
> day 3.0972 1.7599 -0.490
> Residual 190.9227 13.8175
> # of obs: 102, groups: Dog, 10
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 1075.649999 10.521427 102.234
> day 6.091506 1.133983 5.372
> I(day^2) -0.357334 0.051369 -6.956
> SideR -5.401961 2.736268 -1.974
>
> Correlation of Fixed Effects:
> (Intr) day I(d^2)
> day -0.535
> I(day^2) 0.279 -0.795
> SideR -0.130 0.000 0.000
>
> ##?? Fixed effects estimates are sligtly different
> ##?? As df and p-values are missing I assume that it can be concluded
> that as
> ##?? t-value is less than 1.96 that 'Side' is not sigificant.
> ##?? In the lme example for fm4Pixel the t-value for Side is -1.21
> ##?? have i specified fm4Pixel correctly?
>
> -----Original Message-----
> From: Doran, Harold [mailto:HDoran at air.org]
> Sent: Thursday, 9 February 2006 02:01
> To: Paul Cossens; r-help at stat.math.ethz.ch
> Subject: RE: [R] lme syntax for P&B examples
>
>
> Paul:
>
> It is a little difficult to understand what you are trying to translate
> since you do not show what the model would look like using lme. If you
> show lme, then it is easy to translate into lmer syntax.
>
> A few thoughts, first, use lmer in the Matrix package and not in lme4.
> Second, see the Bates article in R news at the link below for dealing
> with nesting structures. Last, a colleague and I have a paper in press
> showing how to fit models using lme which we submitted a year or so ago.
> Since lme has evolved to lmer, we created an appendix that translates
> all of our lme models to the lmer syntax so readers can see
> equivalences. I am happy to send this to you (or others) upon request.
>
> http://cran.r-project.org/doc/Rnews/Rnews_2005-1.pdf
>
> Harold
>
>
>
>
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Paul Cossens
> Sent: Wednesday, February 08, 2006 12:08 AM
> To: r-help at stat.math.ethz.ch
> Subject: [R] lme syntax for P&B examples
>
> Hi helpeRs,
>
> I've been working through some examples in Pinhiero & Bates( 2000)
> trying to understand how to translate to the new Lme4 syntax but without
> much luck. Below is what I think I should do, but either the answers
> don't come out the same or I get errors.
> In the Oxide problems I'm particularly interested in obtaining the
> levels coeficients but this options no longer seems to be available in
> lme4. How can levels infor be obtained in lme4?
>
> If someone can recreate the examples below in lme4 syntax so I can
> follow what is happening in the text I'd be grateful.
>
> Cheers
>
> Paul Cossens
>
>
> #Pixel
> # P&B(2000) p40-45
>
> Pixel<-read.csv("Pixel.csv",header=TRUE);
> Pixel$Side<-as.factor(Pixel$Side)
> Pixel$Dog<-as.factor(Pixel$Dog)
>
> (fm1Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog)+(1|Side), data =
> Pixel))
> (fm2Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog), data = Pixel))
> (fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|Dog:Side), data = Pixel))
> or should I do it this way? Pixel$DS<-with(Pixel,Dog:Side)[drop=TRUE]
> (fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|DS), data = Pixel))
>
> (fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side , data = Pixel))
>
>
> #Oxide
> # P&B(2000) p167-170
>
> Oxide<-read.csv("Oxide.csv",header=TRUE);
> Oxide$Source<-as.factor(Oxide$Source)
> Oxide$Lot<-as.factor(Oxide$Lot)
> Oxide$Wafer<-as.factor(Oxide$Wafer)
> Oxide$Site<-as.factor(Oxide$Site)
> fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide) )
> (fm2Oxide<-lmer(Thickness~ (1|Lot),data=Oxide) )
> coef(fm1Oxide)
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> 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
>
> ______________________________________________
> 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
More information about the R-help
mailing list