Rev 8595 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
R version 4.6.0 (2026-04-24) -- "Because it was There"Copyright (C) 2026 The R Foundation for Statistical ComputingPlatform: x86_64-pc-linux-gnuR is free software and comes with ABSOLUTELY NO WARRANTY.You are welcome to redistribute it under certain conditions.Type 'license()' or 'licence()' for distribution details.R is a collaborative project with many contributors.Type 'contributors()' for more information and'citation()' on how to cite R or R packages in publications.Type 'demo()' for some demos, 'help()' for on-line help, or'help.start()' for an HTML browser interface to help.Type 'q()' to quit R.> ## run reproduction scripts from the NLME book chapters> testdir <- system.file("scripts", package = "nlme", mustWork = TRUE)> scripts <- dir(testdir, pattern = "^ch[0-9]*\\.R$")> for(f in scripts) {+ writeLines(c("", strrep("=", nchar(f)), basename(f), strrep("=", nchar(f))))+ set.seed(3)+ options(warn = 1) # chapters set digits+ source(file.path(testdir, f), echo = TRUE,+ max.deparse.length = Inf, keep.source = TRUE)+ }======ch01.R======> #-*- R -*->> library(nlme)> pdf(file = 'ch01.pdf')> options( width = 65, digits = 5 )> options( contrasts = c(unordered = "contr.helmert", ordered = "contr.poly") )> # Chapter 1 Linear Mixed-Effects Models: Basic Concepts and Examples>> # 1.1 A Simple Example of Random Effects>> RailGrouped Data: travel ~ 1 | RailRail travel1 1 552 1 533 1 544 2 265 2 376 2 327 3 788 3 919 3 8510 4 9211 4 10012 4 9613 5 4914 5 5115 5 5016 6 8017 6 8518 6 83> fm1Rail.lm <- lm( travel ~ 1, data = Rail )> fm1Rail.lmCall:lm(formula = travel ~ 1, data = Rail)Coefficients:(Intercept)66.5> fm2Rail.lm <- lm( travel ~ Rail - 1, data = Rail )> fm2Rail.lmCall:lm(formula = travel ~ Rail - 1, data = Rail)Coefficients:Rail2 Rail5 Rail1 Rail6 Rail3 Rail431.7 50.0 54.0 82.7 84.7 96.0> fm1Rail.lme <- lme(travel ~ 1, data = Rail, random = ~ 1 | Rail)> summary( fm1Rail.lme )Linear mixed-effects model fit by REMLData: RailAIC BIC logLik128.18 130.68 -61.089Random effects:Formula: ~1 | Rail(Intercept) ResidualStdDev: 24.805 4.0208Fixed effects: travel ~ 1Value Std.Error DF t-value p-value(Intercept) 66.5 10.171 12 6.5382 0Standardized Within-Group Residuals:Min Q1 Med Q3 Max-1.618827 -0.282177 0.035693 0.219558 1.614377Number of Observations: 18Number of Groups: 6> fm1Rail.lmeML <- update( fm1Rail.lme, method = "ML" )> summary( fm1Rail.lmeML )Linear mixed-effects model fit by maximum likelihoodData: RailAIC BIC logLik134.56 137.23 -64.28Random effects:Formula: ~1 | Rail(Intercept) ResidualStdDev: 22.624 4.0208Fixed effects: travel ~ 1Value Std.Error DF t-value p-value(Intercept) 66.5 9.554 12 6.9604 0Standardized Within-Group Residuals:Min Q1 Med Q3 Max-1.610981 -0.288870 0.034542 0.213728 1.622223Number of Observations: 18Number of Groups: 6> plot( fm1Rail.lme ) # produces Figure 1.4> intervals( fm1Rail.lme )Approximate 95% confidence intervalsFixed effects:lower est. upper(Intercept) 44.339 66.5 88.661Random Effects:Level: Raillower est. uppersd((Intercept)) 13.274 24.805 46.353Within-group standard error:lower est. upper2.6950 4.0208 5.9987> anova( fm1Rail.lme )numDF denDF F-value p-value(Intercept) 1 12 42.748 <.0001> # 1.2 A Randomized Block Design>> plot.design( ergoStool ) # produces Figure 1.6> contrasts( ergoStool$Type )[,1] [,2] [,3]T1 -1 -1 -1T2 1 -1 -1T3 0 2 -1T4 0 0 3> ergoStool1 <- ergoStool[ ergoStool$Subject == "1", ]> model.matrix( effort ~ Type, ergoStool1 ) # X matrix for Subject 1(Intercept) Type1 Type2 Type31 1 -1 -1 -12 1 1 -1 -13 1 0 2 -14 1 0 0 3attr(,"assign")[1] 0 1 1 1attr(,"contrasts")attr(,"contrasts")$Type[1] "contr.helmert"> fm1Stool <-+ lme(effort ~ Type, data = ergoStool, random = ~ 1 | Subject)> summary( fm1Stool )Linear mixed-effects model fit by REMLData: ergoStoolAIC BIC logLik139.49 148.28 -63.743Random effects:Formula: ~1 | Subject(Intercept) ResidualStdDev: 1.3325 1.1003Fixed effects: effort ~ TypeValue Std.Error DF t-value p-value(Intercept) 10.2500 0.48052 24 21.3309 0.0000Type1 1.9444 0.25934 24 7.4976 0.0000Type2 0.0926 0.14973 24 0.6184 0.5421Type3 -0.3426 0.10588 24 -3.2358 0.0035Correlation:(Intr) Type1 Type2Type1 0Type2 0 0Type3 0 0 0Standardized Within-Group Residuals:Min Q1 Med Q3 Max-1.802003 -0.643166 0.057831 0.700997 1.631421Number of Observations: 36Number of Groups: 9> anova( fm1Stool )numDF denDF F-value p-value(Intercept) 1 24 455.01 <.0001Type 3 24 22.36 <.0001> options( contrasts = c( factor = "contr.treatment",+ ordered = "contr.poly" ) )> contrasts( ergoStool$Type )T2 T3 T4T1 0 0 0T2 1 0 0T3 0 1 0T4 0 0 1> fm2Stool <-+ lme(effort ~ Type, data = ergoStool, random = ~ 1 | Subject)> summary( fm2Stool )Linear mixed-effects model fit by REMLData: ergoStoolAIC BIC logLik133.13 141.93 -60.565Random effects:Formula: ~1 | Subject(Intercept) ResidualStdDev: 1.3325 1.1003Fixed effects: effort ~ TypeValue Std.Error DF t-value p-value(Intercept) 8.5556 0.57601 24 14.8531 0.0000TypeT2 3.8889 0.51868 24 7.4976 0.0000TypeT3 2.2222 0.51868 24 4.2843 0.0003TypeT4 0.6667 0.51868 24 1.2853 0.2110Correlation:(Intr) TypeT2 TypeT3TypeT2 -0.45TypeT3 -0.45 0.50TypeT4 -0.45 0.50 0.50Standardized Within-Group Residuals:Min Q1 Med Q3 Max-1.802003 -0.643166 0.057831 0.700997 1.631421Number of Observations: 36Number of Groups: 9> anova( fm2Stool )numDF denDF F-value p-value(Intercept) 1 24 455.01 <.0001Type 3 24 22.36 <.0001> model.matrix( effort ~ Type - 1, ergoStool1 )TypeT1 TypeT2 TypeT3 TypeT41 1 0 0 02 0 1 0 03 0 0 1 04 0 0 0 1attr(,"assign")[1] 1 1 1 1attr(,"contrasts")attr(,"contrasts")$Type[1] "contr.treatment"> fm3Stool <-+ lme(effort ~ Type - 1, data = ergoStool, random = ~ 1 | Subject)> summary( fm3Stool )Linear mixed-effects model fit by REMLData: ergoStoolAIC BIC logLik133.13 141.93 -60.565Random effects:Formula: ~1 | Subject(Intercept) ResidualStdDev: 1.3325 1.1003Fixed effects: effort ~ Type - 1Value Std.Error DF t-value p-valueTypeT1 8.5556 0.57601 24 14.853 0TypeT2 12.4444 0.57601 24 21.605 0TypeT3 10.7778 0.57601 24 18.711 0TypeT4 9.2222 0.57601 24 16.011 0Correlation:TypeT1 TypeT2 TypeT3TypeT2 0.595TypeT3 0.595 0.595TypeT4 0.595 0.595 0.595Standardized Within-Group Residuals:Min Q1 Med Q3 Max-1.802003 -0.643166 0.057831 0.700997 1.631421Number of Observations: 36Number of Groups: 9> anova( fm3Stool )numDF denDF F-value p-valueType 4 24 130.52 <.0001> intervals( fm1Stool )Approximate 95% confidence intervalsFixed effects:lower est. upper(Intercept) 9.25825 10.250000 11.24175Type1 1.40919 1.944444 2.47970Type2 -0.21644 0.092593 0.40162Type3 -0.56111 -0.342593 -0.12408Random Effects:Level: Subjectlower est. uppersd((Intercept)) 0.74962 1.3325 2.3685Within-group standard error:lower est. upper0.82957 1.10029 1.45937> plot( fm1Stool, # produces Figure 1.8+ form = resid(., type = "p") ~ fitted(.) | Subject,+ abline = 0 )> # 1.3 Mixed-effects Models for Replicated, Blocked Designs>> with(Machines, interaction.plot( Machine, Worker, score, las = 1)) # Figure 1.10> fm1Machine <-+ lme( score ~ Machine, data = Machines, random = ~ 1 | Worker )> fm1MachineLinear mixed-effects model fit by REMLData: MachinesLog-restricted-likelihood: -143.44Fixed: score ~ Machine(Intercept) MachineB MachineC52.3556 7.9667 13.9167Random effects:Formula: ~1 | Worker(Intercept) ResidualStdDev: 5.1466 3.1616Number of Observations: 54Number of Groups: 6> fm2Machine <- update( fm1Machine, random = ~ 1 | Worker/Machine )> fm2MachineLinear mixed-effects model fit by REMLData: MachinesLog-restricted-likelihood: -107.84Fixed: score ~ Machine(Intercept) MachineB MachineC52.3556 7.9667 13.9167Random effects:Formula: ~1 | Worker(Intercept)StdDev: 4.781Formula: ~1 | Machine %in% Worker(Intercept) ResidualStdDev: 3.7295 0.96158Number of Observations: 54Number of Groups:Worker Machine %in% Worker6 18> anova( fm1Machine, fm2Machine )Model df AIC BIC logLik Test L.Ratio p-valuefm1Machine 1 5 296.88 306.54 -143.44fm2Machine 2 6 227.69 239.28 -107.84 1 vs 2 71.191 <.0001> ## delete selected rows from the Machines data> MachinesUnbal <- Machines[ -c(2,3,6,8,9,12,19,20,27,33), ]> ## check that the result is indeed unbalanced> table(MachinesUnbal$Machine, MachinesUnbal$Worker)6 2 4 1 3 5A 3 2 2 1 1 3B 3 3 3 1 2 2C 3 3 3 3 3 3> fm1MachinesU <- lme( score ~ Machine, data = MachinesUnbal,+ random = ~ 1 | Worker/Machine )> fm1MachinesULinear mixed-effects model fit by REMLData: MachinesUnbalLog-restricted-likelihood: -90.936Fixed: score ~ Machine(Intercept) MachineB MachineC52.3540 7.9624 13.9182Random effects:Formula: ~1 | Worker(Intercept)StdDev: 4.7387Formula: ~1 | Machine %in% Worker(Intercept) ResidualStdDev: 3.7728 0.9332Number of Observations: 44Number of Groups:Worker Machine %in% Worker6 18> intervals( fm1MachinesU )Approximate 95% confidence intervalsFixed effects:lower est. upper(Intercept) 47.2345 52.3540 57.474MachineB 3.0278 7.9624 12.897MachineC 8.9955 13.9182 18.841Random Effects:Level: Workerlower est. uppersd((Intercept)) 2.2162 4.7387 10.132Level: Machinelower est. uppersd((Intercept)) 2.4091 3.7728 5.9085Within-group standard error:lower est. upper0.71113 0.93320 1.22463> fm4Stool <- lme( effort ~ Type, ergoStool, ~ 1 | Subject/Type )> if (interactive()) intervals( fm4Stool )> (fm1Stool$sigma)^2[1] 1.2106> (fm4Stool$sigma)^2 + 0.79621^2[1] 0.84554> Machine1 <- Machines[ Machines$Worker == "1", ]> model.matrix( score ~ Machine, Machine1 ) # fixed-effects X_i(Intercept) MachineB MachineC1 1 0 02 1 0 03 1 0 019 1 1 020 1 1 021 1 1 037 1 0 138 1 0 139 1 0 1attr(,"assign")[1] 0 1 1attr(,"contrasts")attr(,"contrasts")$Machine[1] "contr.treatment"> model.matrix( ~ Machine - 1, Machine1 ) # random-effects Z_iMachineA MachineB MachineC1 1 0 02 1 0 03 1 0 019 0 1 020 0 1 021 0 1 037 0 0 138 0 0 139 0 0 1attr(,"assign")[1] 1 1 1attr(,"contrasts")attr(,"contrasts")$Machine[1] "contr.treatment"> fm3Machine <- update( fm1Machine, random = ~Machine - 1 |Worker)> summary( fm3Machine )Linear mixed-effects model fit by REMLData: MachinesAIC BIC logLik228.31 247.63 -104.16Random effects:Formula: ~Machine - 1 | WorkerStructure: General positive-definite, Log-Cholesky parametrizationStdDev CorrMachineA 4.07928 MachnA MachnBMachineB 8.62529 0.803MachineC 4.38948 0.623 0.771Residual 0.96158Fixed effects: score ~ MachineValue Std.Error DF t-value p-value(Intercept) 52.356 1.6807 46 31.1508 0.0000MachineB 7.967 2.4209 46 3.2909 0.0019MachineC 13.917 1.5401 46 9.0362 0.0000Correlation:(Intr) MachnBMachineB 0.463MachineC -0.374 0.301Standardized Within-Group Residuals:Min Q1 Med Q3 Max-2.393540 -0.513776 0.026908 0.472455 2.533387Number of Observations: 54Number of Groups: 6> anova( fm1Machine, fm2Machine, fm3Machine )Model df AIC BIC logLik Test L.Ratio p-valuefm1Machine 1 5 296.88 306.54 -143.44fm2Machine 2 6 227.69 239.28 -107.84 1 vs 2 71.191 <.0001fm3Machine 3 10 228.31 247.63 -104.16 2 vs 3 7.376 0.1173> # 1.4 An Analysis of Covariance Model>> names( Orthodont )[1] "distance" "age" "Subject" "Sex"> levels( Orthodont$Sex )[1] "Male" "Female"> OrthoFem <- Orthodont[ Orthodont$Sex == "Female", ]> fm1OrthF.lis <- lmList( distance ~ age, data = OrthoFem )> coef( fm1OrthF.lis )(Intercept) ageF10 13.55 0.450F09 18.10 0.275F06 17.00 0.375F01 17.25 0.375F05 19.60 0.275F07 16.95 0.550F02 14.20 0.800F08 21.45 0.175F03 14.40 0.850F04 19.65 0.475F11 18.95 0.675> intervals( fm1OrthF.lis ), , (Intercept)lower est. upperF10 10.071 13.55 17.029F09 14.621 18.10 21.579F06 13.521 17.00 20.479F01 13.771 17.25 20.729F05 16.121 19.60 23.079F07 13.471 16.95 20.429F02 10.721 14.20 17.679F08 17.971 21.45 24.929F03 10.921 14.40 17.879F04 16.171 19.65 23.129F11 15.471 18.95 22.429, , agelower est. upperF10 0.1401 0.450 0.7599F09 -0.0349 0.275 0.5849F06 0.0651 0.375 0.6849F01 0.0651 0.375 0.6849F05 -0.0349 0.275 0.5849F07 0.2401 0.550 0.8599F02 0.4901 0.800 1.1099F08 -0.1349 0.175 0.4849F03 0.5401 0.850 1.1599F04 0.1651 0.475 0.7849F11 0.3651 0.675 0.9849> plot( intervals ( fm1OrthF.lis ) ) # produces Figure 1.12> fm2OrthF.lis <- update( fm1OrthF.lis, distance ~ I( age - 11 ) )> plot( intervals( fm2OrthF.lis ) ) # produces Figure 1.13> fm1OrthF <-+ lme( distance ~ age, data = OrthoFem, random = ~ 1 | Subject )> summary( fm1OrthF )Linear mixed-effects model fit by REMLData: OrthoFemAIC BIC logLik149.22 156.17 -70.609Random effects:Formula: ~1 | Subject(Intercept) ResidualStdDev: 2.0685 0.78003Fixed effects: distance ~ ageValue Std.Error DF t-value p-value(Intercept) 17.3727 0.85874 32 20.2304 0age 0.4795 0.05259 32 9.1186 0Correlation:(Intr)age -0.674Standardized Within-Group Residuals:Min Q1 Med Q3 Max-2.27365 -0.70902 0.17282 0.41221 1.63252Number of Observations: 44Number of Groups: 11> fm1OrthFM <- update( fm1OrthF, method = "ML" )> summary( fm1OrthFM )Linear mixed-effects model fit by maximum likelihoodData: OrthoFemAIC BIC logLik146.03 153.17 -69.015Random effects:Formula: ~1 | Subject(Intercept) ResidualStdDev: 1.9699 0.76812Fixed effects: distance ~ ageValue Std.Error DF t-value p-value(Intercept) 17.3727 0.85063 32 20.4234 0age 0.4795 0.05301 32 9.0471 0Correlation:(Intr)age -0.685Standardized Within-Group Residuals:Min Q1 Med Q3 Max-2.30562 -0.71924 0.17636 0.42580 1.66894Number of Observations: 44Number of Groups: 11> fm2OrthF <- update( fm1OrthF, random = ~ age | Subject )> anova( fm1OrthF, fm2OrthF )Model df AIC BIC logLik Test L.Ratio p-valuefm1OrthF 1 4 149.22 156.17 -70.609fm2OrthF 2 6 149.43 159.85 -68.714 1 vs 2 3.7896 0.1503> random.effects( fm1OrthF )(Intercept)F10 -4.005329F09 -1.470449F06 -1.470449F01 -1.229032F05 -0.021947F07 0.340179F02 0.340179F08 0.702304F03 1.064430F04 2.150807F11 3.599309> ranef( fm1OrthFM )(Intercept)F10 -3.995835F09 -1.466964F06 -1.466964F01 -1.226119F05 -0.021895F07 0.339372F02 0.339372F08 0.700640F03 1.061907F04 2.145709F11 3.590778> coef( fm1OrthF )(Intercept) ageF10 13.367 0.47955F09 15.902 0.47955F06 15.902 0.47955F01 16.144 0.47955F05 17.351 0.47955F07 17.713 0.47955F02 17.713 0.47955F08 18.075 0.47955F03 18.437 0.47955F04 19.524 0.47955F11 20.972 0.47955> plot( compareFits(coef(fm1OrthF), coef(fm1OrthFM))) # Figure 1.15> plot( augPred(fm1OrthF), aspect = "xy", grid = TRUE ) # Figure 1.16> # 1.5 Models for Nested Classification Factors>> fm1Pixel <- lme( pixel ~ day + I(day^2), data = Pixel,+ random = list( Dog = ~ day, Side = ~ 1 ) )> intervals( fm1Pixel )Approximate 95% confidence intervalsFixed effects:lower est. upper(Intercept) 1053.0968 1073.33914 1093.5814day 4.3797 6.12960 7.8795I(day^2) -0.4349 -0.36735 -0.2998Random Effects:Level: Doglower est. uppersd((Intercept)) 15.92760 28.36990 50.53187sd(day) 1.08139 1.84375 3.14357cor((Intercept),day) -0.89465 -0.55472 0.19197Level: Sidelower est. uppersd((Intercept)) 10.417 16.824 27.173Within-group standard error:lower est. upper7.6345 8.9896 10.5852> plot( augPred( fm1Pixel ) ) # produces Figure 1.18> VarCorr( fm1Pixel )Variance StdDev CorrDog = pdLogChol(day)(Intercept) 804.8514 28.3699 (Intr)day 3.3994 1.8438 -0.555Side = pdLogChol(1)(Intercept) 283.0572 16.8243Residual 80.8130 8.9896> summary( fm1Pixel )Linear mixed-effects model fit by REMLData: PixelAIC BIC logLik841.21 861.97 -412.61Random effects:Formula: ~day | DogStructure: General positive-definite, Log-Cholesky parametrizationStdDev Corr(Intercept) 28.3699 (Intr)day 1.8438 -0.555Formula: ~1 | Side %in% Dog(Intercept) ResidualStdDev: 16.824 8.9896Fixed effects: pixel ~ day + I(day^2)Value Std.Error DF t-value p-value(Intercept) 1073.34 10.1717 80 105.522 0day 6.13 0.8793 80 6.971 0I(day^2) -0.37 0.0339 80 -10.822 0Correlation:(Intr) dayday -0.517I(day^2) 0.186 -0.668Standardized Within-Group Residuals:Min Q1 Med Q3 Max-2.829057 -0.449181 0.025549 0.557216 2.751965Number of Observations: 102Number of Groups:Dog Side %in% Dog10 20> fm2Pixel <- update( fm1Pixel, random = ~ day | Dog)> anova( fm1Pixel, fm2Pixel )Model df AIC BIC logLik Test L.Ratio p-valuefm1Pixel 1 8 841.21 861.97 -412.61fm2Pixel 2 7 884.52 902.69 -435.26 1 vs 2 45.309 <.0001> fm3Pixel <- update( fm1Pixel, random = ~ 1 | Dog/Side )> anova( fm1Pixel, fm3Pixel )Model df AIC BIC logLik Test L.Ratio p-valuefm1Pixel 1 8 841.21 861.97 -412.61fm3Pixel 2 6 876.84 892.41 -432.42 1 vs 2 39.629 <.0001> fm4Pixel <- update( fm1Pixel, pixel ~ day + I(day^2) + Side )> summary( fm4Pixel )Linear mixed-effects model fit by REMLData: PixelAIC BIC logLik835.85 859.12 -408.93Random effects:Formula: ~day | DogStructure: General positive-definite, Log-Cholesky parametrizationStdDev Corr(Intercept) 28.4636 (Intr)day 1.8438 -0.553Formula: ~1 | Side %in% Dog(Intercept) ResidualStdDev: 16.507 8.9836Fixed effects: pixel ~ day + I(day^2) + SideValue Std.Error DF t-value p-value(Intercept) 1077.95 10.8627 80 99.234 0.0000day 6.13 0.8790 80 6.973 0.0000I(day^2) -0.37 0.0339 80 -10.829 0.0000SideR -9.22 7.6268 9 -1.209 0.2576Correlation:(Intr) day I(d^2)day -0.484I(day^2) 0.174 -0.667SideR -0.351 0.000 0.000Standardized Within-Group Residuals:Min Q1 Med Q3 Max-2.809825 -0.471334 0.026103 0.541154 2.774701Number of Observations: 102Number of Groups:Dog Side %in% Dog10 20> # 1.6 A Split-Plot Experiment>> fm1Oats <- lme( yield ~ ordered(nitro) * Variety, data = Oats,+ random = ~ 1 | Block/Variety )> anova( fm1Oats )numDF denDF F-value p-value(Intercept) 1 45 245.143 <.0001ordered(nitro) 3 45 37.686 <.0001Variety 2 10 1.485 0.2724ordered(nitro):Variety 6 45 0.303 0.9322> fm2Oats <- update( fm1Oats, yield ~ ordered(nitro) + Variety )> anova( fm2Oats )numDF denDF F-value p-value(Intercept) 1 51 245.145 <.0001ordered(nitro) 3 51 41.053 <.0001Variety 2 10 1.485 0.2724> summary( fm2Oats )Linear mixed-effects model fit by REMLData: OatsAIC BIC logLik587.46 607.16 -284.73Random effects:Formula: ~1 | Block(Intercept)StdDev: 14.645Formula: ~1 | Variety %in% Block(Intercept) ResidualStdDev: 10.473 12.75Fixed effects: yield ~ ordered(nitro) + VarietyValue Std.Error DF t-value p-value(Intercept) 104.500 7.7975 51 13.4017 0.0000ordered(nitro).L 32.945 3.0052 51 10.9627 0.0000ordered(nitro).Q -5.167 3.0052 51 -1.7193 0.0916ordered(nitro).C -0.447 3.0052 51 -0.1488 0.8823VarietyMarvellous 5.292 7.0789 10 0.7475 0.4720VarietyVictory -6.875 7.0789 10 -0.9712 0.3544Correlation:(Intr) or().L or().Q or().C VrtyMrordered(nitro).L 0.000ordered(nitro).Q 0.000 0.000ordered(nitro).C 0.000 0.000 0.000VarietyMarvellous -0.454 0.000 0.000 0.000VarietyVictory -0.454 0.000 0.000 0.000 0.500Standardized Within-Group Residuals:Min Q1 Med Q3 Max-1.841341 -0.662797 -0.066943 0.638225 1.660668Number of Observations: 72Number of Groups:Block Variety %in% Block6 18> fm3Oats <- update( fm1Oats, yield ~ ordered( nitro ) )> summary( fm3Oats )Linear mixed-effects model fit by REMLData: OatsAIC BIC logLik597.61 613.14 -291.8Random effects:Formula: ~1 | Block(Intercept)StdDev: 14.506Formula: ~1 | Variety %in% Block(Intercept) ResidualStdDev: 11.039 12.75Fixed effects: yield ~ ordered(nitro)Value Std.Error DF t-value p-value(Intercept) 103.972 6.6407 51 15.6569 0.0000ordered(nitro).L 32.945 3.0052 51 10.9627 0.0000ordered(nitro).Q -5.167 3.0052 51 -1.7193 0.0916ordered(nitro).C -0.447 3.0052 51 -0.1488 0.8823Correlation:(Intr) or().L or().Qordered(nitro).L 0ordered(nitro).Q 0 0ordered(nitro).C 0 0 0Standardized Within-Group Residuals:Min Q1 Med Q3 Max-1.781556 -0.611689 0.022224 0.622007 1.681382Number of Observations: 72Number of Groups:Block Variety %in% Block6 18> fm4Oats <-+ lme( yield ~ nitro, data = Oats, random = ~ 1 | Block/Variety )> summary( fm4Oats )Linear mixed-effects model fit by REMLData: OatsAIC BIC logLik603.04 614.28 -296.52Random effects:Formula: ~1 | Block(Intercept)StdDev: 14.506Formula: ~1 | Variety %in% Block(Intercept) ResidualStdDev: 11.005 12.867Fixed effects: yield ~ nitroValue Std.Error DF t-value p-value(Intercept) 81.872 6.9453 53 11.788 0nitro 73.667 6.7815 53 10.863 0Correlation:(Intr)nitro -0.293Standardized Within-Group Residuals:Min Q1 Med Q3 Max-1.743808 -0.664752 0.017104 0.542988 1.802989Number of Observations: 72Number of Groups:Block Variety %in% Block6 18> VarCorr( fm4Oats )Variance StdDevBlock = pdLogChol(1)(Intercept) 210.42 14.506Variety = pdLogChol(1)(Intercept) 121.10 11.005Residual 165.56 12.867> intervals( fm4Oats )Approximate 95% confidence intervalsFixed effects:lower est. upper(Intercept) 67.942 81.872 95.803nitro 60.065 73.667 87.269Random Effects:Level: Blocklower est. uppersd((Intercept)) 6.6089 14.506 31.839Level: Varietylower est. uppersd((Intercept)) 6.4081 11.005 18.898Within-group standard error:lower est. upper10.637 12.867 15.565> plot(augPred(fm4Oats), aspect = 2.5, layout = c(6, 3),+ between = list(x = c(0, 0, 0.5, 0, 0))) # produces Figure 1.21> # cleanup>> summary(warnings())No warnings======ch02.R======> #-*- R -*->> library( nlme )> options( width = 65, digits = 5 )> options( contrasts = c(unordered = "contr.helmert",+ ordered = "contr.poly") )> pdf( file = 'ch02.pdf' )> # Chapter 2 Theory and Computational Methods for Linear Mixed-Effects Models>> # 2.2 Likelihood Estimation for LME Models>> Xmat <- matrix( c(1, 1, 1, 1, 8, 10, 12, 14), ncol = 2 )> Xmat[,1] [,2][1,] 1 8[2,] 1 10[3,] 1 12[4,] 1 14> Xqr <- qr( Xmat ) # creates a QR structure> qr.R( Xqr ) # returns R[,1] [,2][1,] -2 -22.0000[2,] 0 -4.4721> qr.Q( Xqr ) # returns Q-truncated[,1] [,2][1,] -0.5 0.67082[2,] -0.5 0.22361[3,] -0.5 -0.22361[4,] -0.5 -0.67082> qr.Q( Xqr, complete = TRUE ) # returns the full Q[,1] [,2] [,3] [,4][1,] -0.5 0.67082 0.023607 0.54721[2,] -0.5 0.22361 -0.439345 -0.71202[3,] -0.5 -0.22361 0.807869 -0.21760[4,] -0.5 -0.67082 -0.392131 0.38240> fm1Rail.lme <- lme( travel ~ 1, data = Rail, random = ~ 1 | Rail,+ control = list( msVerbose = TRUE ) )0: 61.048859: -1.819591: 61.048859: -1.81959> fm1Rail.lme <- lme( travel ~ 1, data = Rail, random = ~ 1 | Rail,+ control = list( msVerbose = TRUE, niterEM = 0 ))0: 67.893737: -0.4315231: 61.612483: -1.431522: 61.138913: -1.984413: 61.050114: -1.838664: 61.048866: -1.818195: 61.048859: -1.819606: 61.048859: -1.81959> fm1Machine <-+ lme( score ~ Machine, data = Machines, random = ~ 1 | Worker )> fm2Machine <- update( fm1Machine, random = ~ 1 | Worker/Machine )> anova( fm1Machine, fm2Machine )Model df AIC BIC logLik Test L.Ratio p-valuefm1Machine 1 5 300.46 310.12 -145.23fm2Machine 2 6 231.27 242.86 -109.64 1 vs 2 71.191 <.0001> OrthoFem <- Orthodont[ Orthodont$Sex == "Female", ]> fm1OrthF <- lme( distance ~ age, data = OrthoFem,+ random = ~ 1 | Subject )> fm2OrthF <- update( fm1OrthF, random = ~ age | Subject )> orthLRTsim <- simulate.lme( fm1OrthF, m2 = fm2OrthF, nsim = 1000 )> plot( orthLRTsim, df = c(1, 2) ) # produces Figure 2.3> machineLRTsim <- simulate.lme(fm1Machine, m2 = fm2Machine, nsim= 1000)> plot( machineLRTsim, df = c(0, 1), # produces Figure 2.4+ layout = c(4,1), between = list(x = c(0, 0.5, 0)) )> stoolLRTsim <-+ simulate.lme( list(fixed = effort ~ 1, data = ergoStool,+ random = ~ 1 | Subject),+ m2 = list(fixed = effort ~ Type),+ method = "ML", nsim = 1000 )> plot( stoolLRTsim, df = c(3, 4) ) # Figure 2.5> ## "partially balanced incomplete block" experiment> ## from Littell et al. 1996 (Data Set 1.5.1):> ##data( PBIB, package = 'SASmixed' ) # reproduced below> PBIB <- data.frame(+ "response" = c(+ 2.4, 2.5, 2.6, 2.0, 2.7, 2.8, 2.4, 2.7, 2.6, 2.8, 2.4, 2.4,+ 3.4, 3.1, 2.1, 2.3, 4.1, 3.3, 3.3, 2.9, 3.4, 3.2, 2.8, 3.0,+ 3.2, 2.5, 2.4, 2.6, 2.3, 2.3, 2.4, 2.7, 2.8, 2.8, 2.6, 2.5,+ 2.5, 2.7, 2.8, 2.6, 2.6, 2.6, 2.3, 2.4, 2.7, 2.7, 2.5, 2.6,+ 3.0, 3.6, 3.2, 3.2, 3.0, 2.8, 2.4, 2.5, 2.4, 2.5, 3.2, 3.1+ ),+ "Treatment" = factor(c(+ 15, 9, 1, 13, 5, 7, 8, 1, 10, 1, 14, 2, 15, 11, 2, 3,+ 6, 15, 4, 7, 12, 4, 3, 1, 12, 14, 15, 8, 6, 3, 14, 5,+ 5, 4, 2, 13, 10, 12, 13, 6, 9, 7, 10, 3, 8, 6, 2, 9,+ 5, 9, 11, 12, 7, 13, 14, 11, 10, 4, 8, 11+ )),+ "Block" = as.factor(rep(1:15, each = 4))+ )> pbibLRTsim <-+ simulate.lme(list( fixed = response ~ 1, data = PBIB,+ random = ~ 1 | Block ),+ m2 = list(fixed = response ~ Treatment, data = PBIB,+ random = ~ 1 | Block),+ method = "ML", nsim = 1000 )> plot( pbibLRTsim, df = c(14,16,18), weights = FALSE ) # Figure 2.6> summary( fm2Machine )Linear mixed-effects model fit by REMLData: MachinesAIC BIC logLik231.27 242.86 -109.64Random effects:Formula: ~1 | Worker(Intercept)StdDev: 4.781Formula: ~1 | Machine %in% Worker(Intercept) ResidualStdDev: 3.7295 0.96158Fixed effects: score ~ MachineValue Std.Error DF t-value p-value(Intercept) 59.650 2.14467 36 27.8131 0.0000Machine1 3.983 1.08849 10 3.6595 0.0044Machine2 3.311 0.62844 10 5.2688 0.0004Correlation:(Intr) Machn1Machine1 0Machine2 0 0Standardized Within-Group Residuals:Min Q1 Med Q3 Max-2.269587 -0.548466 -0.010706 0.439366 2.540058Number of Observations: 54Number of Groups:Worker Machine %in% Worker6 18> fm1PBIB <- lme(response ~ Treatment, data = PBIB, random = ~ 1 | Block)> anova( fm1PBIB )numDF denDF F-value p-value(Intercept) 1 31 1654.21 <.0001Treatment 14 31 1.53 0.1576> fm2PBIB <- update( fm1PBIB, method = "ML" )> fm3PBIB <- update( fm2PBIB, response ~ 1 )> anova( fm2PBIB, fm3PBIB )Model df AIC BIC logLik Test L.Ratio p-valuefm2PBIB 1 17 56.571 92.174 -11.285fm3PBIB 2 3 52.152 58.435 -23.076 1 vs 2 23.581 0.0514> anova( fm2Machine )numDF denDF F-value p-value(Intercept) 1 36 773.57 <.0001Machine 2 10 20.58 3e-04> ##save(orthLRTsim, machineLRTsim, pbibLRTsim, stoolLRTsim,> ## file = "sims.rda")>> summary(warnings())No warnings======ch03.R======> #-*- R -*->> # initialization>> library(nlme)> options(width = 65, digits = 5)> options(contrasts = c(unordered = "contr.helmert", ordered = "contr.poly"))> pdf(file = 'ch03.pdf')> # Chapter 3 Describing the Structure of Grouped Data>> # 3.1 The Display Formula and Its Components>> formula( Rail )travel ~ 1 | Rail> formula( ergoStool )effort ~ Type | Subject> formula( Machines )score ~ Machine | Worker> formula( Orthodont )distance ~ age | Subject> formula( Pixel )pixel ~ day | Dog/Side> formula( Oats )yield ~ nitro | Block> table( Oxboys$Subject )10 26 25 9 2 6 7 17 16 15 8 20 1 18 5 23 11 21 3 24 229 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 912 13 14 19 49 9 9 9 9> table( getGroups( Oxboys ) )10 26 25 9 2 6 7 17 16 15 8 20 1 18 5 23 11 21 3 24 229 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 912 13 14 19 49 9 9 9 9> unique( table( getGroups( Oxboys ) ) ) # a more concise result[1] 9> unique( table( getCovariate( Oxboys ), getGroups( Oxboys ) ) )[1] 1 0> length( unique( getCovariate( Oxboys ) ) )[1] 16> unique( getGroups(Pixel, level = 1) )[1] 1 2 3 4 5 6 7 8 9 10Levels: 1 10 2 3 4 5 6 7 8 9> unique( getGroups(Pixel, level = 2) )[1] 1/R 2/R 3/R 4/R 5/R 6/R 7/R 8/R 9/R 10/R 1/L 2/L[13] 3/L 4/L 5/L 6/L 7/L 8/L 9/L 10/L20 Levels: 1/R < 2/R < 3/R < 4/R < 5/R < 6/R < 7/R < ... < 10/L> Pixel.groups <- getGroups( Pixel, level = 1:2 )> class( Pixel.groups )[1] "data.frame"> names( Pixel.groups )[1] "Dog" "Side"> unique( Pixel.groups[["Side"]] )[1] R LLevels: L R> formula( PBG )deltaBP ~ dose | Rabbit> PBG.log <- update( PBG, formula = deltaBP ~ log(dose) | Rabbit )> formula(PBG.log)deltaBP ~ log(dose) | Rabbit<environment: 0x55aa561e7958>> unique( getCovariate(PBG.log) )[1] 1.8326 2.5257 3.2189 3.9120 4.6052 5.2983> unique( getCovariate(PBG) )[1] 6.25 12.50 25.00 50.00 100.00 200.00> # 3.2 Constructing groupedData Objects>> # The next line is not from the book.> # It is added to ensure that the file is available>> write.table( Oxboys, "oxboys.dat" )> Oxboys.frm <- read.table( "oxboys.dat", header = TRUE )> class( Oxboys.frm ) # check the class of the result[1] "data.frame"> dim( Oxboys.frm ) # check the dimensions[1] 234 4> Oxboys <- groupedData( height ~ age | Subject,+ data = read.table("oxboys.dat", header = TRUE),+ labels = list(x = "Centered age", y = "Height"),+ units = list(y = "(cm)") )> Oxboys # display the objectGrouped Data: height ~ age | SubjectSubject age height Occasion1 1 -1.0000 140.50 12 1 -0.7479 143.40 23 1 -0.4630 144.80 34 1 -0.1643 147.10 45 1 -0.0027 147.70 56 1 0.2466 150.20 67 1 0.5562 151.70 78 1 0.7781 153.30 89 1 0.9945 155.80 910 2 -1.0000 136.90 111 2 -0.7479 139.10 212 2 -0.4630 140.10 313 2 -0.1643 142.60 414 2 -0.0027 143.20 515 2 0.2466 144.00 616 2 0.5562 145.80 717 2 0.7781 146.80 818 2 0.9945 148.30 919 3 -1.0000 150.00 120 3 -0.7479 152.10 221 3 -0.4630 153.90 322 3 -0.1643 155.80 423 3 -0.0027 156.00 524 3 0.2466 156.90 625 3 0.5562 157.40 726 3 0.7781 159.10 827 3 0.9945 160.60 928 4 -1.0000 155.70 129 4 -0.7479 158.70 230 4 -0.4630 160.60 331 4 -0.1643 163.30 432 4 -0.0027 164.40 533 4 0.2466 167.30 634 4 0.5562 170.70 735 4 0.7781 172.00 836 4 0.9945 174.80 937 5 -1.0000 145.80 138 5 -0.7479 147.30 239 5 -0.4493 148.70 340 5 -0.1643 149.78 441 5 -0.0027 150.22 542 5 0.2466 152.50 643 5 0.5562 154.80 744 5 0.7781 156.40 845 5 0.9973 158.70 946 6 -1.0000 142.40 147 6 -0.7479 143.80 248 6 -0.4630 145.20 349 6 -0.1643 146.30 450 6 -0.0027 147.10 551 6 0.2466 148.10 652 6 0.5562 148.90 753 6 0.7781 149.10 854 6 0.9945 151.00 955 7 -1.0000 141.30 156 7 -0.7479 142.40 257 7 -0.4493 144.30 358 7 -0.1643 145.20 459 7 0.0000 146.10 560 7 0.2466 146.80 661 7 0.5562 147.90 762 7 0.7945 150.50 863 7 0.9945 151.80 964 8 -1.0000 141.70 165 8 -0.7479 143.70 266 8 -0.4630 145.10 367 8 -0.1643 147.90 468 8 -0.0027 148.10 569 8 0.2466 149.60 670 8 0.5562 150.99 771 8 0.7945 154.10 872 8 1.0055 154.90 973 9 -1.0000 132.70 174 9 -0.7479 134.10 275 9 -0.4493 135.30 376 9 -0.1643 136.60 477 9 -0.0027 137.50 578 9 0.2466 139.10 679 9 0.5562 140.90 780 9 0.7945 143.70 881 9 0.9945 144.70 982 10 -1.0000 126.20 183 10 -0.7479 128.20 284 10 -0.4630 129.00 385 10 -0.1643 129.40 486 10 -0.0027 129.59 587 10 0.2466 130.60 688 10 0.5562 132.50 789 10 0.7781 133.40 890 10 0.9945 134.20 991 11 -1.0000 142.50 192 11 -0.7479 143.80 293 11 -0.4630 145.60 394 11 -0.1643 148.30 495 11 -0.0027 149.40 596 11 0.2466 151.60 697 11 0.5562 154.80 798 11 0.7781 156.90 899 11 0.9945 159.20 9100 12 -1.0000 149.90 1101 12 -0.7479 151.70 2102 12 -0.4630 153.30 3103 12 -0.1643 156.10 4104 12 0.0000 156.70 5105 12 0.2466 157.80 6106 12 0.5562 160.70 7107 12 0.7781 162.70 8108 12 0.9945 163.80 9109 13 -1.0000 148.90 1110 13 -0.7150 149.80 2111 13 -0.4630 151.70 3112 13 -0.1643 154.40 4113 13 -0.0027 155.50 5114 13 0.2466 156.40 6115 13 0.5562 161.40 7116 13 0.7781 163.90 8117 13 0.9945 164.60 9118 14 -1.0000 151.60 1119 14 -0.7479 153.20 2120 14 -0.4630 155.20 3121 14 -0.1643 157.30 4122 14 0.0000 159.10 5123 14 0.2466 160.90 6124 14 0.5562 164.40 7125 14 0.7781 166.90 8126 14 0.9945 168.40 9127 15 -1.0000 137.50 1128 15 -0.7479 139.30 2129 15 -0.4630 140.90 3130 15 -0.1643 142.70 4131 15 -0.0027 144.20 5132 15 0.2466 145.70 6133 15 0.5562 147.09 7134 15 0.7781 150.20 8135 15 0.9945 152.30 9136 16 -1.0000 142.80 1137 16 -0.7479 144.90 2138 16 -0.4630 145.00 3139 16 -0.1643 146.70 4140 16 -0.0027 147.20 5141 16 0.2466 148.90 6142 16 0.5562 150.10 7143 16 0.7781 151.00 8144 16 0.9945 152.20 9145 17 -1.0000 134.90 1146 17 -0.7479 137.40 2147 17 -0.4630 138.20 3148 17 -0.1643 140.20 4149 17 -0.0027 143.60 5150 17 0.2466 144.20 6151 17 0.5562 147.90 7152 17 0.7781 150.30 8153 17 0.9945 151.80 9154 18 -1.0000 145.50 1155 18 -0.7479 146.20 2156 18 -0.4630 148.20 3157 18 -0.1643 150.30 4158 18 -0.0027 152.00 5159 18 0.2466 152.30 6160 18 0.5562 154.30 7161 18 0.7781 156.20 8162 18 0.9945 156.80 9163 19 -1.0000 156.90 1164 19 -0.7479 157.90 2165 19 -0.4630 160.30 3166 19 -0.1643 161.90 4167 19 0.0000 163.80 5168 19 0.2466 165.50 6169 19 0.5562 169.90 7170 19 0.7781 172.40 8171 19 0.9945 174.40 9172 20 -1.0000 146.50 1173 20 -0.7479 148.40 2174 20 -0.4630 149.30 3175 20 -0.1643 151.20 4176 20 -0.0027 152.10 5177 20 0.2466 152.40 6178 20 0.5562 153.90 7179 20 0.7781 154.90 8180 20 0.9945 155.40 9181 21 -1.0000 143.90 1182 21 -0.7479 145.10 2183 21 -0.4630 147.00 3184 21 -0.1643 149.20 4185 21 -0.0027 149.80 5186 21 0.2466 151.50 6187 21 0.5562 153.17 7188 21 0.7781 156.90 8189 21 0.9945 159.60 9190 22 -1.0000 147.40 1191 22 -0.7479 148.80 2192 22 -0.4630 150.10 3193 22 -0.1643 152.50 4194 22 -0.0027 154.70 5195 22 0.2466 156.00 6196 22 0.5562 158.40 7197 22 0.7781 161.50 8198 22 0.9945 163.30 9199 23 -1.0000 144.50 1200 23 -0.7479 146.00 2201 23 -0.4630 147.40 3202 23 -0.1643 149.20 4203 23 -0.0027 150.80 5204 23 0.2466 152.50 6205 23 0.5562 155.00 7206 23 0.7781 156.80 8207 23 0.9945 158.80 9208 24 -1.0000 147.80 1209 24 -0.7479 148.20 2210 24 -0.4630 150.20 3211 24 -0.1643 151.00 4212 24 -0.0027 152.20 5213 24 0.2466 153.60 6214 24 0.5562 155.80 7215 24 0.7781 159.20 8216 24 0.9945 161.60 9217 25 -1.0000 135.50 1218 25 -0.7479 136.60 2219 25 -0.4630 137.30 3220 25 -0.1643 138.20 4221 25 -0.0027 139.00 5222 25 0.2466 139.50 6223 25 0.5562 141.00 7224 25 0.7808 142.70 8225 25 0.9945 143.90 9226 26 -1.0000 132.20 1227 26 -0.7479 134.30 2228 26 -0.4630 135.10 3229 26 -0.1643 136.70 4230 26 -0.0027 138.40 5231 26 0.2466 138.90 6232 26 0.5562 141.80 7233 26 0.7781 142.60 8234 26 1.0055 143.10 9> unique( getGroups( Oxboys ) )[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20[21] 21 22 23 24 25 2626 Levels: 10 < 26 < 25 < 9 < 2 < 6 < 7 < 17 < 16 < 15 < ... < 4> plot( BodyWeight, outer = ~ Diet, aspect = 3 ) # Figure 3.3> plot( BodyWeight, outer = TRUE, aspect = 3 )> plot( Soybean, outer = ~ Year * Variety ) # Figure 6.10> plot( Soybean, outer = ~ Variety * Year )> gsummary( BodyWeight, invar = TRUE )Rat Diet2 2 13 3 14 4 11 1 18 8 15 5 16 6 17 7 111 11 29 9 210 10 212 12 213 13 315 15 314 14 316 16 3> plot( PBG, inner = ~ Treatment, scales = list(x = list(log = 2)))> ergoStool.mat <- asTable( ergoStool )> ergoStool.matT1 T2 T3 T48 7 11 8 75 8 11 8 74 7 11 10 99 9 13 10 86 9 11 11 103 7 14 13 97 8 12 12 111 12 15 12 102 10 14 13 12> ergoStool.new <- balancedGrouped( effort ~ Type | Subject,+ data = ergoStool.mat )> ergoStool.newGrouped Data: effort ~ Type | SubjectType Subject effort1 T1 8 72 T2 8 113 T3 8 84 T4 8 75 T1 5 86 T2 5 117 T3 5 88 T4 5 79 T1 4 710 T2 4 1111 T3 4 1012 T4 4 913 T1 9 914 T2 9 1315 T3 9 1016 T4 9 817 T1 6 918 T2 6 1119 T3 6 1120 T4 6 1021 T1 3 722 T2 3 1423 T3 3 1324 T4 3 925 T1 7 826 T2 7 1227 T3 7 1228 T4 7 1129 T1 1 1230 T2 1 1531 T3 1 1232 T4 1 1033 T1 2 1034 T2 2 1435 T3 2 1336 T4 2 12> # 3.3 Controlling Trellis Graphics Presentations of Grouped Data>> plot(CO2, layout=c(6,2), between=list(x=c(0,0,0.5,0,0)))> plot( Spruce, layout = c(7, 4, 3),+ skip = c(rep(FALSE, 27), TRUE, rep(FALSE, 27), TRUE,+ rep(FALSE, 12), rep(TRUE, 2), rep(FALSE,13)) )> plot( Spruce, layout = c(9, 3, 3),+ skip = c(rep(FALSE, 66), TRUE, TRUE, rep(FALSE, 13)) )> unique( getCovariate(DNase) )[1] 0.048828 0.195312 0.390625 0.781250 1.562500 3.125000[7] 6.250000 12.500000> log( unique(getCovariate(DNase)), 2 )[1] -4.35614 -2.35614 -1.35614 -0.35614 0.64386 1.64386[7] 2.64386 3.64386> plot( DNase, layout=c(6,2), scales = list(x=list(log=2)) )> plot(Pixel, layout = c(4,5),+ between = list(x = c(0, 0.5, 0), y = 0.5))> plot( Pixel, displayLevel = 1 )> plot( Wafer, display = 1, collapse = 1 )> plot( Wafer, display = 1, collapse = 1,+ FUN = function(x) sqrt(var(x)), layout = c(10,1) )> # 3.4 Summaries>> sapply( ergoStool, data.class )effort Type Subject"numeric" "factor" "ordered"> gsummary( Theoph, inv = TRUE )Subject Wt Dose6 6 80.0 4.007 7 64.6 4.958 8 70.5 4.5311 11 65.0 4.923 3 70.5 4.532 2 72.4 4.404 4 72.7 4.409 9 86.4 3.1012 12 60.5 5.3010 10 58.2 5.501 1 79.6 4.025 5 54.6 5.86> gsummary( Theoph, omit = TRUE, inv = TRUE )Wt Dose6 80.0 4.007 64.6 4.958 70.5 4.5311 65.0 4.923 70.5 4.532 72.4 4.404 72.7 4.409 86.4 3.1012 60.5 5.3010 58.2 5.501 79.6 4.025 54.6 5.86> is.null(gsummary(Theoph, inv = TRUE, omit = TRUE)) # invariants present[1] FALSE> is.null(gsummary(Oxboys, inv = TRUE, omit = TRUE)) # no invariants[1] TRUE> gsummary( Theoph )Subject Wt Dose Time conc6 6 80.0 4.00 5.8882 3.52557 7 64.6 4.95 5.8655 3.91098 8 70.5 4.53 5.8900 4.271811 11 65.0 4.92 5.8718 4.51093 3 70.5 4.53 5.9073 5.08642 2 72.4 4.40 5.8691 4.82364 4 72.7 4.40 5.9400 4.94009 9 86.4 3.10 5.8682 4.893612 12 60.5 5.30 5.8764 5.410010 10 58.2 5.50 5.9155 5.93091 1 79.6 4.02 5.9500 6.43915 5 54.6 5.86 5.8936 5.7827> gsummary( Theoph, FUN = max, omit = TRUE )Wt Dose Time conc6 80.0 4.00 23.85 6.447 64.6 4.95 24.22 7.098 70.5 4.53 24.12 7.5611 65.0 4.92 24.08 8.003 70.5 4.53 24.17 8.202 72.4 4.40 24.30 8.334 72.7 4.40 24.65 8.609 86.4 3.10 24.43 9.0312 60.5 5.30 24.15 9.7510 58.2 5.50 23.70 10.211 79.6 4.02 24.37 10.505 54.6 5.86 24.35 11.40> Quin.sum <- gsummary( Quinidine, omit = TRUE, FUN = mean )> dim( Quin.sum )[1] 136 13> Quin.sum[1:10, ]time conc dose interval Age Height Weight Race109 30.2633 NA NA NA 70 67 58.000 Caucasian70 0.7500 NA NA NA 68 69 75.000 Caucasian23 52.0262 NA NA NA 75 72 108.000 Caucasian92 8.8571 NA NA NA 68 72 65.000 Caucasian111 18.1638 NA NA NA 68 66 56.000 Latin5 24.3750 NA NA NA 62 71 66.000 Caucasian18 196.8438 NA NA NA 87 69 85.375 Caucasian24 31.2500 NA NA NA 55 69 89.000 Latin2 12.2000 NA NA NA 58 69 85.000 Latin88 4.7900 NA NA NA 85 72 77.000 CaucasianSmoke Ethanol Heart Creatinine glyco109 no none No/Mild >= 50 0.4600070 no former No/Mild >= 50 1.1500023 yes none No/Mild >= 50 0.8387592 yes former No/Mild >= 50 1.27000111 yes former No/Mild >= 50 1.230005 yes none Severe >= 50 1.3900018 no none No/Mild < 50 1.2600024 no former No/Mild >= 50 0.570002 no current Moderate >= 50 0.8200088 no none Moderate >= 50 0.61000> Quinidine[Quinidine[["Subject"]] == 3, 1:8]Grouped Data: conc ~ time | SubjectSubject time conc dose interval Age Height Weight17 3 0.00 NA 201 NA 67 69 6918 3 8.00 NA 201 NA 67 69 6919 3 16.00 NA 201 NA 67 69 6920 3 24.00 NA 201 NA 67 69 6921 3 32.00 NA 201 NA 67 69 6922 3 41.25 2.4 NA NA 67 69 6923 3 104.00 NA 201 8 67 69 6924 3 113.00 2.3 NA NA 67 69 6925 3 3865.00 NA 201 6 67 69 6226 3 3873.00 NA 201 NA 67 69 6227 3 3881.00 NA 201 NA 67 69 6228 3 3889.00 NA 201 NA 67 69 6229 3 3897.00 NA 201 NA 67 69 6230 3 3900.00 NA NA NA 67 69 6231 3 3905.00 NA 201 NA 67 69 6232 3 3909.00 4.7 NA NA 67 69 6233 3 4073.00 NA 201 8 67 69 62> Quin.sum1 <- gsummary( Quinidine, omit = TRUE )> Quin.sum1[1:10, 1:7]time conc dose interval Age Height Weight109 30.2633 0.50000 100.00 NaN 70 67 58.00070 0.7500 0.60000 201.00 8 68 69 75.00023 52.0262 0.56667 166.00 6 75 72 108.00092 8.8571 0.70000 83.00 NaN 68 72 65.000111 18.1638 0.90000 249.00 NaN 68 66 56.0005 24.3750 0.70000 301.00 NaN 62 71 66.00018 196.8438 0.93333 201.00 6 87 69 85.37524 31.2500 1.10000 187.88 NaN 55 69 89.0002 12.2000 1.20000 166.00 NaN 58 69 85.00088 4.7900 1.20000 201.00 8 85 72 77.000> summary( Quin.sum1 )time conc dose intervalMin. : 0.065 Min. :0.50 Min. : 83 Min. : 5.001st Qu.: 19.285 1st Qu.:1.70 1st Qu.:198 1st Qu.: 6.00Median : 47.169 Median :2.24 Median :201 Median : 6.00Mean : 251.486 Mean :2.36 Mean :224 Mean : 6.993rd Qu.: 171.148 3rd Qu.:2.92 3rd Qu.:249 3rd Qu.: 8.00Max. :5364.890 Max. :5.77 Max. :498 Max. :12.00NAs :29Age Height Weight RaceMin. :42.0 Min. :60.0 Min. : 41.0 Caucasian:911st Qu.:61.0 1st Qu.:67.0 1st Qu.: 67.8 Latin :35Median :66.0 Median :70.0 Median : 79.2 Black :10Mean :66.9 Mean :69.6 Mean : 79.23rd Qu.:73.0 3rd Qu.:72.0 3rd Qu.: 88.2Max. :92.0 Max. :79.0 Max. :119.0Smoke Ethanol Heart Creatinine glycono :94 none :90 No/Mild :55 < 50 : 36 Min. :0.390yes:42 current:16 Moderate:40 >= 50:100 1st Qu.:0.885former :30 Severe :41 Median :1.174Mean :1.2123rd Qu.:1.453Max. :2.995> summary( Quinidine )Subject time conc dose223 : 47 Min. : 0 Min. :0.40 Min. : 83110 : 41 1st Qu.: 16 1st Qu.:1.60 1st Qu.:16681 : 40 Median : 60 Median :2.30 Median :201136 : 32 Mean : 373 Mean :2.45 Mean :2257 : 31 3rd Qu.: 241 3rd Qu.:3.00 3rd Qu.:24976 : 28 Max. :8096 Max. :9.40 Max. :603(Other):1252 NAs :1110 NAs :443interval Age Height WeightMin. : 4.00 Min. :42.0 Min. :60.0 Min. : 41.01st Qu.: 6.00 1st Qu.:60.0 1st Qu.:67.0 1st Qu.: 69.5Median : 6.00 Median :66.0 Median :69.0 Median : 78.0Mean : 7.11 Mean :66.7 Mean :69.2 Mean : 79.73rd Qu.: 8.00 3rd Qu.:74.0 3rd Qu.:72.0 3rd Qu.: 89.0Max. :12.00 Max. :92.0 Max. :79.0 Max. :119.0NAs :1222Race Smoke Ethanol HeartCaucasian:968 no :1024 none :991 No/Mild :598Latin :384 yes: 447 current:191 Moderate:375Black :119 former :289 Severe :498Creatinine glyco< 50 : 418 Min. :0.39>= 50:1053 1st Qu.:0.93Median :1.23Mean :1.283rd Qu.:1.54Max. :3.16> sum( ifelse(is.na(Quinidine[["conc"]]), 0, 1) )[1] 361> sum( !is.na(Quinidine[["conc"]]) )[1] 361> sum( !is.na(Quinidine[["dose"]]) )[1] 1028> gapply( Quinidine, "conc", function(x) sum(!is.na(x)) )109 70 23 92 111 5 18 24 2 88 91 117 120 13 89 271 1 3 1 1 2 3 1 1 1 1 3 2 1 3 153 122 129 132 16 106 15 22 57 77 115 121 123 11 48 1261 1 2 3 1 1 1 1 3 1 4 1 1 2 2 2223 19 38 42 52 56 63 83 104 118 137 17 29 34 46 736 1 1 2 1 1 4 1 2 2 1 1 1 1 3 287 103 138 45 44 97 36 37 72 100 8 71 6 14 26 752 1 2 3 7 2 2 3 1 3 1 5 1 3 1 320 96 99 134 12 49 67 85 112 127 55 68 124 1 35 472 3 2 1 1 3 3 1 3 3 6 3 1 2 2 579 95 114 135 105 116 62 65 107 130 66 139 33 80 125 1103 3 2 2 1 3 4 7 4 3 1 3 3 2 1 11128 136 21 43 90 102 40 84 98 30 82 93 108 119 32 1332 11 2 1 1 2 2 6 2 1 3 4 1 3 1 27 9 76 94 58 113 50 39 78 25 61 3 64 60 59 106 2 6 5 1 2 3 2 10 2 2 3 4 4 3 669 4 81 54 41 74 28 512 6 11 4 3 3 4 6> table( gapply(Quinidine, "conc", function(x) sum(!is.na(x))) )1 2 3 4 5 6 7 10 1146 33 31 9 3 8 2 1 3> changeRecords <- gapply( Quinidine, FUN = function(frm)+ any(is.na(frm[["conc"]]) & is.na(frm[["dose"]])) )> changeRecords109 70 23 92 111 5 18 24 2 88FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE91 117 120 13 89 27 53 122 129 132FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE16 106 15 22 57 77 115 121 123 11FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE48 126 223 19 38 42 52 56 63 83TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE104 118 137 17 29 34 46 73 87 103FALSE TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE138 45 44 97 36 37 72 100 8 71FALSE TRUE TRUE TRUE FALSE TRUE FALSE FALSE FALSE TRUE6 14 26 75 20 96 99 134 12 49FALSE TRUE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE67 85 112 127 55 68 124 1 35 47FALSE FALSE TRUE TRUE TRUE FALSE FALSE FALSE FALSE TRUE79 95 114 135 105 116 62 65 107 130TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE FALSE FALSE66 139 33 80 125 110 128 136 21 43FALSE TRUE TRUE TRUE FALSE TRUE FALSE TRUE FALSE FALSE90 102 40 84 98 30 82 93 108 119FALSE FALSE TRUE TRUE TRUE FALSE TRUE FALSE FALSE TRUE32 133 7 9 76 94 58 113 50 39FALSE TRUE TRUE FALSE TRUE TRUE FALSE FALSE TRUE FALSE78 25 61 3 64 60 59 10 69 4FALSE FALSE TRUE TRUE TRUE FALSE FALSE TRUE FALSE TRUE81 54 41 74 28 51TRUE TRUE TRUE FALSE TRUE FALSE> sort( as.numeric( names(changeRecords)[changeRecords] ) )[1] 3 4 7 10 14 18 28 33 37 40 41 44 45 46 47[16] 48 50 54 55 61 62 63 64 65 71 75 76 77 79 80[31] 81 82 84 94 95 96 97 98 110 112 114 118 119 127 132[46] 133 135 136 139 223> Quinidine[29:31,]Grouped Data: conc ~ time | SubjectSubject time conc dose interval Age Height Weight Race29 3 3897 NA 201 NA 67 69 62 Caucasian30 3 3900 NA NA NA 67 69 62 Caucasian31 3 3905 NA 201 NA 67 69 62 CaucasianSmoke Ethanol Heart Creatinine glyco29 yes former Moderate < 50 1.7130 yes former Moderate < 50 1.7131 yes former Moderate < 50 1.71> Quinidine[Quinidine[["Subject"]] == 4, ]Grouped Data: conc ~ time | SubjectSubject time conc dose interval Age Height Weight Race45 4 0.00 NA 332 NA 88 66 103 Black46 4 7.00 NA 332 NA 88 66 103 Black47 4 13.00 NA 332 NA 88 66 103 Black48 4 19.00 NA 332 NA 88 66 103 Black49 4 21.50 3.1 NA NA 88 66 103 Black50 4 85.00 NA 249 6 88 66 103 Black51 4 91.00 5.8 NA NA 88 66 103 Black52 4 91.08 NA 249 NA 88 66 103 Black53 4 97.00 NA 249 NA 88 66 103 Black54 4 103.00 NA 249 NA 88 66 103 Black55 4 105.00 NA NA NA 88 66 92 Black56 4 109.00 NA 249 NA 88 66 92 Black57 4 115.00 NA 249 NA 88 66 92 Black58 4 145.00 NA 166 NA 88 66 92 Black59 4 151.00 NA 166 NA 88 66 92 Black60 4 156.00 3.1 NA NA 88 66 92 Black61 4 157.00 NA 166 NA 88 66 92 Black62 4 163.00 NA 166 NA 88 66 92 Black63 4 169.00 NA 166 NA 88 66 92 Black64 4 174.75 NA 201 NA 88 66 92 Black65 4 177.00 NA NA NA 88 66 92 Black66 4 181.50 3.1 NA NA 88 66 92 Black67 4 245.00 NA 201 8 88 66 92 Black68 4 249.00 NA NA NA 88 66 86 Black69 4 252.50 3.2 NA NA 88 66 86 Black70 4 317.00 NA 201 8 88 66 86 Black71 4 326.00 1.9 NA NA 88 66 86 BlackSmoke Ethanol Heart Creatinine glyco45 yes none Severe >= 50 1.4846 yes none Severe >= 50 1.4847 yes none Severe >= 50 1.4848 yes none Severe >= 50 1.4849 yes none Severe >= 50 1.4850 yes none Severe >= 50 1.6151 yes none Severe >= 50 1.6152 yes none Severe >= 50 1.6153 yes none Severe >= 50 1.6154 yes none Severe >= 50 1.6155 yes none Severe >= 50 1.6156 yes none Severe >= 50 1.6157 yes none Severe >= 50 1.6158 yes none Severe >= 50 1.8859 yes none Severe >= 50 1.8860 yes none Severe >= 50 1.8861 yes none Severe >= 50 1.8862 yes none Severe >= 50 1.8863 yes none Severe >= 50 1.8864 yes none Severe >= 50 1.8865 yes none Severe >= 50 1.6866 yes none Severe >= 50 1.6867 yes none Severe >= 50 1.8768 yes none Severe >= 50 1.8769 yes none Severe >= 50 1.8770 yes none Severe >= 50 1.8371 yes none Severe >= 50 1.83> # cleanup>> summary(warnings())No warnings======ch04.R======> #-*- R -*->> # initialization>> library(nlme)> library(lattice)> options(width = 65,+ ## reduce platform dependence in printed output when testing+ digits = if(nzchar(Sys.getenv("R_TESTS"))) 3 else 5)> options(contrasts = c(unordered = "contr.helmert", ordered = "contr.poly"))> pdf(file = 'ch04.pdf')> # Chapter 4 Fitting Linear Mixed-Effects Models>> # 4.1 Fitting Linear Models in S with lm and lmList>> fm1Orth.lm <- lm(distance ~ age, Orthodont)> fm1Orth.lmCall:lm(formula = distance ~ age, data = Orthodont)Coefficients:(Intercept) age16.76 0.66> par(mfrow=c(2,2))> plot(fm1Orth.lm) # Figure 4.1> fm2Orth.lm <- update(fm1Orth.lm, formula = distance ~ Sex*age)> summary(fm2Orth.lm)Call:lm(formula = distance ~ Sex + age + Sex:age, data = Orthodont)Residuals:Min 1Q Median 3Q Max-5.616 -1.322 -0.168 1.330 5.247Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 16.8567 1.1094 15.19 < 2e-16 ***Sex1 0.5161 1.1094 0.47 0.64age 0.6320 0.0988 6.39 4.7e-09 ***Sex1:age -0.1524 0.0988 -1.54 0.13---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 2.26 on 104 degrees of freedomMultiple R-squared: 0.423, Adjusted R-squared: 0.406F-statistic: 25.4 on 3 and 104 DF, p-value: 2.11e-12> fm3Orth.lm <- update(fm2Orth.lm, formula = . ~ . - Sex)> summary(fm3Orth.lm)Call:lm(formula = distance ~ age + Sex:age, data = Orthodont)Residuals:Min 1Q Median 3Q Max-5.742 -1.242 -0.189 1.268 5.267Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 16.7611 1.0861 15.43 < 2e-16 ***age 0.6403 0.0968 6.61 1.6e-09 ***age:Sex1 -0.1074 0.0196 -5.47 3.0e-07 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 2.25 on 105 degrees of freedomMultiple R-squared: 0.422, Adjusted R-squared: 0.411F-statistic: 38.3 on 2 and 105 DF, p-value: 3.31e-13> bwplot(getGroups(Orthodont)~resid(fm2Orth.lm)) # Figure 4.2> fm1Orth.lis <- lmList(distance ~ age | Subject, Orthodont)> getGroupsFormula(Orthodont)~Subject> fm1Orth.lis <- lmList(distance ~ age, Orthodont)> formula(Orthodont)distance ~ age | Subject> fm1Orth.lis <- lmList(Orthodont)> fm1Orth.lisCall:Model: distance ~ age | SubjectData: OrthodontCoefficients:(Intercept) ageM16 17.0 0.550M05 13.7 0.850M02 14.9 0.775M11 20.1 0.325M07 15.0 0.800M08 19.8 0.375M03 16.0 0.750M12 13.2 1.000M13 2.8 1.950M14 19.1 0.525M09 14.4 0.975M15 13.5 1.125M06 19.0 0.675M04 24.7 0.175M01 17.3 0.950M10 21.2 0.750F10 13.5 0.450F09 18.1 0.275F06 17.0 0.375F01 17.2 0.375F05 19.6 0.275F07 16.9 0.550F02 14.2 0.800F08 21.4 0.175F03 14.4 0.850F04 19.7 0.475F11 19.0 0.675Degrees of freedom: 108 total; 54 residualResidual standard error: 1.31> summary(fm1Orth.lis)Call:Model: distance ~ age | SubjectData: OrthodontCoefficients:(Intercept)Estimate Std. Error t value Pr(>|t|)M16 17.0 3.29 5.155 3.70e-06M05 13.7 3.29 4.151 1.18e-04M02 14.9 3.29 4.516 3.46e-05M11 20.1 3.29 6.098 1.19e-07M07 15.0 3.29 4.547 3.12e-05M08 19.8 3.29 6.006 1.67e-07M03 16.0 3.29 4.866 1.03e-05M12 13.2 3.29 4.030 1.76e-04M13 2.8 3.29 0.852 3.98e-01M14 19.1 3.29 5.809 3.45e-07M09 14.4 3.29 4.379 5.51e-05M15 13.5 3.29 4.106 1.37e-04M06 19.0 3.29 5.763 4.08e-07M04 24.7 3.29 7.512 6.08e-10M01 17.3 3.29 5.261 2.52e-06M10 21.2 3.29 6.463 3.07e-08F10 13.5 3.29 4.121 1.31e-04F09 18.1 3.29 5.505 1.05e-06F06 17.0 3.29 5.170 3.50e-06F01 17.2 3.29 5.246 2.67e-06F05 19.6 3.29 5.961 1.97e-07F07 16.9 3.29 5.155 3.70e-06F02 14.2 3.29 4.319 6.76e-05F08 21.4 3.29 6.523 2.44e-08F03 14.4 3.29 4.379 5.51e-05F04 19.7 3.29 5.976 1.86e-07F11 19.0 3.29 5.763 4.08e-07ageEstimate Std. Error t value Pr(>|t|)M16 0.550 0.293 1.878 6.58e-02M05 0.850 0.293 2.902 5.36e-03M02 0.775 0.293 2.646 1.07e-02M11 0.325 0.293 1.109 2.72e-01M07 0.800 0.293 2.731 8.51e-03M08 0.375 0.293 1.280 2.06e-01M03 0.750 0.293 2.560 1.33e-02M12 1.000 0.293 3.414 1.22e-03M13 1.950 0.293 6.657 1.49e-08M14 0.525 0.293 1.792 7.87e-02M09 0.975 0.293 3.328 1.58e-03M15 1.125 0.293 3.840 3.25e-04M06 0.675 0.293 2.304 2.51e-02M04 0.175 0.293 0.597 5.53e-01M01 0.950 0.293 3.243 2.03e-03M10 0.750 0.293 2.560 1.33e-02F10 0.450 0.293 1.536 1.30e-01F09 0.275 0.293 0.939 3.52e-01F06 0.375 0.293 1.280 2.06e-01F01 0.375 0.293 1.280 2.06e-01F05 0.275 0.293 0.939 3.52e-01F07 0.550 0.293 1.878 6.58e-02F02 0.800 0.293 2.731 8.51e-03F08 0.175 0.293 0.597 5.53e-01F03 0.850 0.293 2.902 5.36e-03F04 0.475 0.293 1.622 1.11e-01F11 0.675 0.293 2.304 2.51e-02Residual standard error: 1.31 on 54 degrees of freedom> pairs(fm1Orth.lis, id = 0.01, adj = -0.5) # Figure 4.3> fm2Orth.lis <- update(fm1Orth.lis, distance ~ I(age-11))> intervals(fm2Orth.lis), , (Intercept)lower est. upperM16 21.7 23.0 24.3M05 21.7 23.0 24.3M02 22.1 23.4 24.7M11 22.3 23.6 24.9M07 22.4 23.8 25.1M08 22.6 23.9 25.2M03 22.9 24.2 25.6M12 22.9 24.2 25.6M13 22.9 24.2 25.6M14 23.6 24.9 26.2M09 23.8 25.1 26.4M15 24.6 25.9 27.2M06 25.1 26.4 27.7M04 25.3 26.6 27.9M01 26.4 27.8 29.1M10 28.2 29.5 30.8F10 17.2 18.5 19.8F09 19.8 21.1 22.4F06 19.8 21.1 22.4F01 20.1 21.4 22.7F05 21.3 22.6 23.9F07 21.7 23.0 24.3F02 21.7 23.0 24.3F08 22.1 23.4 24.7F03 22.4 23.8 25.1F04 23.6 24.9 26.2F11 25.1 26.4 27.7, , I(age - 11)lower est. upperM16 -0.0373 0.550 1.137M05 0.2627 0.850 1.437M02 0.1877 0.775 1.362M11 -0.2623 0.325 0.912M07 0.2127 0.800 1.387M08 -0.2123 0.375 0.962M03 0.1627 0.750 1.337M12 0.4127 1.000 1.587M13 1.3627 1.950 2.537M14 -0.0623 0.525 1.112M09 0.3877 0.975 1.562M15 0.5377 1.125 1.712M06 0.0877 0.675 1.262M04 -0.4123 0.175 0.762M01 0.3627 0.950 1.537M10 0.1627 0.750 1.337F10 -0.1373 0.450 1.037F09 -0.3123 0.275 0.862F06 -0.2123 0.375 0.962F01 -0.2123 0.375 0.962F05 -0.3123 0.275 0.862F07 -0.0373 0.550 1.137F02 0.2127 0.800 1.387F08 -0.4123 0.175 0.762F03 0.2627 0.850 1.437F04 -0.1123 0.475 1.062F11 0.0877 0.675 1.262> plot(intervals(fm2Orth.lis)) # Figure 4.5> IGFGrouped Data: conc ~ age | LotLot age conc1 1 7 4.902 1 7 5.683 1 8 5.324 1 8 5.505 1 13 4.946 1 13 5.197 1 14 5.188 1 14 5.679 1 15 5.0210 1 15 5.8811 1 22 5.1212 1 23 5.2413 1 24 5.8814 1 27 5.4015 1 28 5.5916 1 28 5.7717 1 30 5.5718 1 34 5.8619 1 34 5.8720 1 35 4.6521 1 35 5.3422 1 36 4.9323 1 36 5.3324 1 36 4.9925 1 41 3.3826 1 42 5.4427 1 42 5.2428 1 43 5.3929 2 3 5.3430 2 3 5.2731 2 3 5.4832 2 6 5.1533 2 11 4.2334 2 11 5.7735 2 11 5.0636 2 12 5.3337 2 12 5.7838 2 13 5.0139 2 13 4.8540 2 13 4.9441 2 18 5.1442 2 24 5.4343 2 24 5.6644 2 25 5.6245 2 25 5.5346 2 26 6.2047 2 27 5.3048 2 27 4.0949 2 32 5.7850 2 32 5.6651 2 34 5.0752 2 38 5.4553 2 40 4.7654 2 42 4.8155 2 45 4.9256 2 46 4.3257 2 47 3.3058 3 1 5.8859 3 2 5.9160 3 5 0.8661 3 6 5.4062 3 7 4.9463 3 8 5.4264 3 13 5.4065 3 15 5.6866 3 15 5.7167 3 21 9.5568 3 21 5.9469 3 21 6.1770 3 22 5.3471 3 22 8.1472 3 27 5.5173 3 28 5.3174 3 28 4.8175 3 28 5.2676 3 29 4.7277 3 30 5.0878 3 30 3.9979 3 33 4.8780 3 34 4.9281 3 34 6.1382 3 35 6.3083 3 36 5.9784 3 37 5.9885 3 41 6.6886 3 42 5.3387 3 43 6.0888 3 44 4.7689 3 47 5.3190 3 47 6.6691 3 48 5.5292 3 49 5.4893 3 50 5.1094 4 5 5.1295 4 5 5.0896 4 5 4.6397 4 5 5.3898 4 7 5.7899 4 9 9.34100 4 11 5.58101 4 11 5.19102 4 12 5.25103 4 12 5.44104 4 14 5.31105 4 14 4.71106 4 14 5.67107 4 14 4.65108 4 14 5.05109 4 15 4.23110 4 19 5.02111 4 19 4.98112 4 20 5.08113 4 20 4.84114 4 22 4.84115 4 22 5.53116 4 25 5.85117 4 25 5.32118 4 26 5.47119 5 1 5.49120 5 2 5.43121 5 6 5.02122 5 6 5.29123 5 7 6.25124 5 9 4.63125 5 10 5.18126 5 15 5.17127 5 15 4.98128 5 15 5.38129 5 15 3.76130 5 17 5.63131 5 21 6.12132 5 22 4.00133 5 23 6.53134 5 24 4.67135 5 24 5.55136 5 24 5.62137 5 29 4.58138 5 30 5.41139 5 35 4.84140 5 37 4.83141 5 37 5.36142 5 37 4.81143 5 37 5.35144 5 42 5.46145 5 43 5.09146 5 44 4.78147 5 44 4.44148 5 45 4.67149 5 48 4.98150 6 2 4.56151 6 3 5.83152 6 3 5.27153 6 4 4.90154 7 1 4.94155 7 2 4.78156 7 3 5.42157 7 4 5.42158 7 5 5.38159 7 7 5.55160 7 10 5.81161 7 10 5.62162 7 11 6.08163 7 15 4.80164 7 16 5.32165 7 17 4.95166 7 17 5.44167 7 18 5.48168 7 21 5.26169 7 22 5.21170 7 23 4.65171 7 24 4.62172 7 24 5.15173 7 26 4.71174 7 27 5.02175 7 29 5.38176 7 31 5.34177 7 31 5.10178 7 32 5.69179 7 36 5.00180 7 37 5.02181 7 38 9.74182 7 38 9.60183 7 39 5.58184 7 42 4.94185 7 43 4.66186 7 43 5.23187 7 45 5.62188 7 45 5.53189 7 45 5.45190 7 45 4.63191 7 47 5.01192 7 50 5.43193 8 1 6.17194 8 1 5.57195 8 2 4.82196 8 3 5.84197 8 6 5.55198 8 9 5.17199 8 9 6.50200 8 9 5.36201 9 4 5.47202 9 4 5.57203 9 5 5.36204 9 7 4.93205 9 8 5.49206 9 11 3.25207 9 13 5.53208 9 13 4.91209 9 13 5.74210 9 14 4.95211 9 15 5.07212 9 19 5.54213 9 20 5.29214 9 21 4.59215 9 25 5.66216 9 26 4.69217 9 26 5.18218 9 27 5.19219 9 27 5.35220 9 29 5.28221 9 29 5.50222 9 29 5.00223 9 30 5.47224 9 33 5.55225 9 34 5.75226 9 35 5.41227 9 35 5.65228 9 35 5.25229 9 36 5.81230 9 40 4.71231 9 41 4.95232 10 4 6.00233 10 5 5.74234 10 6 5.68235 10 6 5.83236 10 11 5.30237 10 13 5.63> plot(IGF) # Figure 4.6> fm1IGF.lis <- lmList(IGF)> coef(fm1IGF.lis)(Intercept) age9 5.10 0.005736 4.63 0.170001 5.49 -0.0077910 6.05 -0.047332 5.48 -0.014438 5.59 0.006065 5.37 -0.009514 5.58 -0.016663 5.28 0.010087 5.21 0.00931> plot(intervals(fm1IGF.lis)) # Figure 4.7> fm1IGF.lm <- lm(conc ~ age, data = IGF)> summary(fm1IGF.lm)Call:lm(formula = conc ~ age, data = IGF)Residuals:Min 1Q Median 3Q Max-4.488 -0.374 -0.009 0.258 4.414Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 5.351059 0.103734 51.58 <2e-16 ***age -0.000669 0.003943 -0.17 0.87---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 0.833 on 235 degrees of freedomMultiple R-squared: 0.000123, Adjusted R-squared: -0.00413F-statistic: 0.0288 on 1 and 235 DF, p-value: 0.865> # 4.2 Fitting Linear Mixed-Effects Models with lme>> fm1Orth.lme <- lme(distance ~ I(age-11), data = Orthodont,+ random = ~ I(age-11) | Subject)> fm1Orth.lme <- lme(distance ~ I(age-11), data = Orthodont)> fm1Orth.lme <- lme(fm2Orth.lis)> fm1Orth.lmeLinear mixed-effects model fit by REMLData: OrthodontLog-restricted-likelihood: -221Fixed: distance ~ I(age - 11)(Intercept) I(age - 11)24.02 0.66Random effects:Formula: ~I(age - 11) | SubjectStructure: General positive-definite, Log-Cholesky parametrizationStdDev Corr(Intercept) 2.134 (Intr)I(age - 11) 0.226 0.503Residual 1.310Number of Observations: 108Number of Groups: 27> fm2Orth.lme <- update(fm1Orth.lme, distance~Sex*I(age-11))> summary(fm2Orth.lme)Linear mixed-effects model fit by REMLData: OrthodontAIC BIC logLik451 473 -218Random effects:Formula: ~I(age - 11) | SubjectStructure: General positive-definite, Log-Cholesky parametrizationStdDev Corr(Intercept) 1.83 (Intr)I(age - 11) 0.18 0.206Residual 1.31Fixed effects: distance ~ Sex + I(age - 11) + Sex:I(age - 11)Value Std.Error DF t-value p-value(Intercept) 23.81 0.381 79 62.5 0.0000Sex1 -1.16 0.381 25 -3.0 0.0054I(age - 11) 0.63 0.067 79 9.4 0.0000Sex1:I(age - 11) -0.15 0.067 79 -2.3 0.0264Correlation:(Intr) Sex1 I(-11)Sex1 0.185I(age - 11) 0.102 0.019Sex1:I(age - 11) 0.019 0.102 0.185Standardized Within-Group Residuals:Min Q1 Med Q3 Max-3.1681 -0.3859 0.0071 0.4452 3.8495Number of Observations: 108Number of Groups: 27> fitted(fm2Orth.lme, level = 0:1)fixed Subject1 22.6 24.82 24.2 26.63 25.8 28.34 27.3 30.05 22.6 21.36 24.2 22.87 25.8 24.38 27.3 25.89 22.6 22.010 24.2 23.611 25.8 25.112 27.3 26.613 22.6 24.514 24.2 25.815 25.8 27.016 27.3 28.317 22.6 20.918 24.2 22.519 25.8 24.020 27.3 25.621 22.6 23.922 24.2 25.423 25.8 27.024 27.3 28.525 22.6 21.626 24.2 23.127 25.8 24.728 27.3 26.229 22.6 22.030 24.2 23.331 25.8 24.632 27.3 26.033 22.6 22.634 24.2 24.335 25.8 26.036 27.3 27.637 22.6 26.538 24.2 28.139 25.8 29.840 27.3 31.541 22.6 21.842 24.2 23.143 25.8 24.444 27.3 25.745 22.6 21.846 24.2 23.547 25.8 25.248 27.3 26.849 22.6 21.250 24.2 23.351 25.8 25.552 27.3 27.753 22.6 22.754 24.2 24.255 25.8 25.656 27.3 27.057 22.6 23.158 24.2 24.959 25.8 26.760 27.3 28.561 22.6 21.162 24.2 22.563 25.8 23.964 27.3 25.365 21.2 20.266 22.2 21.167 23.1 21.968 24.1 22.869 21.2 21.370 22.2 22.471 23.1 23.672 24.1 24.773 21.2 21.974 22.2 23.175 23.1 24.276 24.1 25.477 21.2 23.178 22.2 24.179 23.1 25.180 24.1 26.181 21.2 21.382 22.2 22.283 23.1 23.084 24.1 23.985 21.2 20.086 22.2 20.987 23.1 21.788 24.1 22.689 21.2 21.590 22.2 22.591 23.1 23.592 24.1 24.593 21.2 22.094 22.2 22.995 23.1 23.796 24.1 24.597 21.2 20.198 22.2 20.999 23.1 21.7100 24.1 22.5101 21.2 17.7102 22.2 18.6103 23.1 19.4104 24.1 20.2105 21.2 24.2106 22.2 25.4107 23.1 26.5108 24.1 27.7> resid(fm2Orth.lme, level = 1)M01 M01 M01 M01 M02 M02 M021.15428 -1.57649 0.69275 0.96198 0.22522 -0.29641 -1.31803M02 M03 M03 M03 M03 M04 M040.66034 0.96689 -1.06449 -1.09588 0.87274 1.03549 1.74867M04 M04 M05 M05 M05 M05 M06-0.53814 -1.32495 -0.90249 1.04571 -1.50610 0.44210 0.61473M06 M06 M06 M07 M07 M07 M070.06728 0.01983 -0.02762 0.42649 -1.11840 -0.16330 0.29181M08 M08 M08 M08 M09 M09 M092.00813 -1.81291 -0.13395 -0.45500 0.39248 -3.78229 5.04295M09 M10 M10 M10 M10 M11 M11-1.63182 1.02728 -0.14284 1.18705 0.01693 1.18276 -0.10495M11 M11 M12 M12 M12 M12 M13-0.89267 -0.68038 -0.34919 -0.01420 -1.17920 1.15579 -4.15031M13 M13 M13 M14 M14 M14 M141.17692 0.50416 1.83139 -0.22716 1.34520 -0.08244 -1.01008M15 M15 M15 M15 M16 M16 M16-0.13140 -0.40616 -0.68091 1.54433 0.87681 -1.01465 -0.40610M16 F01 F01 F01 F01 F02 F02-0.29756 0.79027 -1.07931 -0.44889 0.18152 -0.27124 -0.91092F02 F02 F03 F03 F03 F03 F040.44940 0.80973 -1.36869 0.94509 0.25887 0.57265 0.40409F04 F04 F04 F05 F05 F05 F050.38858 -0.12694 0.35754 0.15965 0.81049 -0.53868 -0.38784F06 F06 F06 F06 F07 F07 F070.00168 0.13870 -0.72427 -0.08725 0.04484 0.03879 -0.46727F07 F08 F08 F08 F08 F09 F090.52667 0.95185 0.13632 -0.17921 -0.49475 -0.07189 0.11859F09 F09 F10 F10 F10 F10 F110.30906 -1.00047 -1.22334 0.44296 -0.39073 -0.72443 0.28277F11 F11 F11-0.37929 1.45866 0.29661attr(,"label")[1] "Residuals (mm)"> resid(fm2Orth.lme, level = 1, type = "pearson")M01 M01 M01 M01 M02 M02 M020.88111 -1.20339 0.52880 0.73431 0.17192 -0.22626 -1.00610M02 M03 M03 M03 M03 M04 M040.50406 0.73806 -0.81257 -0.83652 0.66619 0.79042 1.33482M04 M04 M05 M05 M05 M05 M06-0.41078 -1.01139 -0.68890 0.79822 -1.14966 0.33747 0.46925M06 M06 M06 M07 M07 M07 M070.05136 0.01514 -0.02108 0.32556 -0.85372 -0.12465 0.22275M08 M08 M08 M08 M09 M09 M091.53288 -1.38386 -0.10225 -0.34732 0.29959 -2.88715 3.84946M09 M10 M10 M10 M10 M11 M11-1.24562 0.78416 -0.10903 0.90612 0.01293 0.90284 -0.08012M11 M11 M12 M12 M12 M12 M13-0.68140 -0.51936 -0.26655 -0.01084 -0.90013 0.88226 -3.16808M13 M13 M13 M14 M14 M14 M140.89839 0.38484 1.39796 -0.17340 1.02684 -0.06293 -0.77103M15 M15 M15 M15 M16 M16 M16-0.10030 -0.31003 -0.51977 1.17884 0.66930 -0.77452 -0.30999M16 F01 F01 F01 F01 F02 F02-0.22714 0.60324 -0.82388 -0.34266 0.13856 -0.20705 -0.69534F02 F02 F03 F03 F03 F03 F040.34305 0.61809 -1.04477 0.72142 0.19761 0.43712 0.30846F04 F04 F04 F05 F05 F05 F050.29661 -0.09690 0.27293 0.12187 0.61867 -0.41119 -0.29605F06 F06 F06 F06 F07 F07 F070.00128 0.10588 -0.55286 -0.06660 0.03423 0.02961 -0.35668F07 F08 F08 F08 F08 F09 F090.40203 0.72658 0.10406 -0.13680 -0.37766 -0.05488 0.09052F09 F09 F10 F10 F10 F10 F110.23592 -0.76369 -0.93382 0.33813 -0.29826 -0.55298 0.21585F11 F11 F11-0.28952 1.11345 0.22641attr(,"label")[1] "Standardized residuals"> newOrth <- data.frame(Subject = rep(c("M11","F03"), c(3, 3)),+ Sex = rep(c("Male", "Female"), c(3, 3)),+ age = rep(16:18, 2))> predict(fm2Orth.lme, newdata = newOrth)M11 M11 M11 F03 F03 F0327.0 27.6 28.3 26.6 27.2 27.8attr(,"label")[1] "Predicted values (mm)"> predict(fm2Orth.lme, newdata = newOrth, level = 0:1)Subject predict.fixed predict.Subject1 M11 28.9 27.02 M11 29.7 27.63 M11 30.5 28.34 F03 25.0 26.65 F03 25.5 27.26 F03 26.0 27.8> fm2Orth.lmeM <- update(fm2Orth.lme, method = "ML")> summary(fm2Orth.lmeM)Linear mixed-effects model fit by maximum likelihoodData: OrthodontAIC BIC logLik444 465 -214Random effects:Formula: ~I(age - 11) | SubjectStructure: General positive-definite, Log-Cholesky parametrizationStdDev Corr(Intercept) 1.752 (Intr)I(age - 11) 0.154 0.234Residual 1.310Fixed effects: distance ~ Sex + I(age - 11) + Sex:I(age - 11)Value Std.Error DF t-value p-value(Intercept) 23.81 0.373 79 63.8 0.0000Sex1 -1.16 0.373 25 -3.1 0.0046I(age - 11) 0.63 0.066 79 9.6 0.0000Sex1:I(age - 11) -0.15 0.066 79 -2.3 0.0237Correlation:(Intr) Sex1 I(-11)Sex1 0.185I(age - 11) 0.102 0.019Sex1:I(age - 11) 0.019 0.102 0.185Standardized Within-Group Residuals:Min Q1 Med Q3 Max-3.3360 -0.4154 0.0104 0.4917 3.8582Number of Observations: 108Number of Groups: 27> compOrth <-+ compareFits(coef(fm2Orth.lis), coef(fm1Orth.lme))> compOrth, , (Intercept)coef(fm2Orth.lis) coef(fm1Orth.lme)M16 23.0 23.1M05 23.0 23.1M02 23.4 23.5M11 23.6 23.6M07 23.8 23.8M08 23.9 23.8M03 24.2 24.2M12 24.2 24.3M13 24.2 24.4M14 24.9 24.8M09 25.1 25.1M15 25.9 25.8M06 26.4 26.2M04 26.6 26.3M01 27.8 27.4M10 29.5 29.0F10 18.5 19.0F09 21.1 21.3F06 21.1 21.4F01 21.4 21.6F05 22.6 22.7F07 23.0 23.1F02 23.0 23.1F08 23.4 23.4F03 23.8 23.8F04 24.9 24.8F11 26.4 26.2, , I(age - 11)coef(fm2Orth.lis) coef(fm1Orth.lme)M16 0.550 0.591M05 0.850 0.686M02 0.775 0.675M11 0.325 0.541M07 0.800 0.695M08 0.375 0.565M03 0.750 0.696M12 1.000 0.775M13 1.950 1.074M14 0.525 0.646M09 0.975 0.796M15 1.125 0.868M06 0.675 0.743M04 0.175 0.594M01 0.950 0.876M10 0.750 0.871F10 0.450 0.410F09 0.275 0.442F06 0.375 0.474F01 0.375 0.482F05 0.275 0.492F07 0.550 0.591F02 0.800 0.670F08 0.175 0.486F03 0.850 0.711F04 0.475 0.630F11 0.675 0.743> plot(compOrth, mark = fixef(fm1Orth.lme)) # Figure 4.8> ## Figure 4.9> plot(comparePred(fm2Orth.lis, fm1Orth.lme, length.out = 2),+ layout = c(8,4), between = list(y = c(0, 0.5, 0)))> plot(compareFits(ranef(fm2Orth.lme), ranef(fm2Orth.lmeM)),+ mark = c(0, 0))> fm4Orth.lm <- lm(distance ~ Sex * I(age-11), Orthodont)> summary(fm4Orth.lm)Call:lm(formula = distance ~ Sex * I(age - 11), data = Orthodont)Residuals:Min 1Q Median 3Q Max-5.616 -1.322 -0.168 1.330 5.247Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 23.8082 0.2210 107.73 < 2e-16 ***Sex1 -1.1605 0.2210 -5.25 8.1e-07 ***I(age - 11) 0.6320 0.0988 6.39 4.7e-09 ***Sex1:I(age - 11) -0.1524 0.0988 -1.54 0.13---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 2.26 on 104 degrees of freedomMultiple R-squared: 0.423, Adjusted R-squared: 0.406F-statistic: 25.4 on 3 and 104 DF, p-value: 2.11e-12> anova(fm2Orth.lme, fm4Orth.lm)Model df AIC BIC logLik Test L.Ratio p-valuefm2Orth.lme 1 8 451 473 -218fm4Orth.lm 2 5 496 510 -243 1 vs 2 51 <.0001> #fm1IGF.lme <- lme(fm1IGF.lis)> #fm1IGF.lme> #intervals(fm1IGF.lme)> #summary(fm1IGF.lme)> pd1 <- pdDiag(~ age)> pd1Uninitialized positive definite matrix structure of class pdDiag.> formula(pd1)~age> #fm2IGF.lme <- update(fm1IGF.lme, random = pdDiag(~age))> (fm2IGF.lme <- lme(conc ~ age, IGF,+ random = pdDiag(~age)))Linear mixed-effects model fit by REMLData: IGFLog-restricted-likelihood: -297Fixed: conc ~ age(Intercept) age5.36904 -0.00193Random effects:Formula: ~age | LotStructure: Diagonal(Intercept) age ResidualStdDev: 3.62e-05 0.00537 0.822Number of Observations: 237Number of Groups: 10> #anova(fm1IGF.lme, fm2IGF.lme)> anova(fm2IGF.lme)numDF denDF F-value p-value(Intercept) 1 226 6439 <.0001age 1 226 0 0.673> #update(fm1IGF.lme, random = list(Lot = pdDiag(~ age)))> pd2 <- pdDiag(value = diag(2), form = ~ age)> pd2Positive definite matrix structure of class pdDiag representing[,1] [,2][1,] 1 0[2,] 0 1> formula(pd2)~age> lme(conc ~ age, IGF, pdDiag(diag(2), ~age))Linear mixed-effects model fit by REMLData: IGFLog-restricted-likelihood: -297Fixed: conc ~ age(Intercept) age5.36904 -0.00193Random effects:Formula: ~age | LotStructure: Diagonal(Intercept) age ResidualStdDev: 3.12e-05 0.00537 0.822Number of Observations: 237Number of Groups: 10> fm4OatsB <- lme(yield ~ nitro, data = Oats,+ random =list(Block = pdCompSymm(~ Variety - 1)))> summary(fm4OatsB)Linear mixed-effects model fit by REMLData: OatsAIC BIC logLik603 614 -297Random effects:Formula: ~Variety - 1 | BlockStructure: Compound SymmetryStdDev CorrVarietyGolden Rain 18.2VarietyMarvellous 18.2 0.635VarietyVictory 18.2 0.635 0.635Residual 12.9Fixed effects: yield ~ nitroValue Std.Error DF t-value p-value(Intercept) 81.9 6.95 65 11.8 0nitro 73.7 6.78 65 10.9 0Correlation:(Intr)nitro -0.293Standardized Within-Group Residuals:Min Q1 Med Q3 Max-1.7438 -0.6648 0.0171 0.5430 1.8030Number of Observations: 72Number of Groups: 6> corMatrix(fm4OatsB$modelStruct$reStruct$Block)[1,2][1] 0.635> fm4OatsC <- lme(yield ~ nitro, data = Oats,+ random=list(Block=pdBlocked(list(pdIdent(~ 1),+ pdIdent(~ Variety-1)))))> summary(fm4OatsC)Linear mixed-effects model fit by REMLData: OatsAIC BIC logLik603 614 -297Random effects:Composite Structure: BlockedBlock 1: (Intercept)Formula: ~1 | Block(Intercept)StdDev: 14.5Block 2: VarietyGolden Rain, VarietyMarvellous, VarietyVictoryFormula: ~Variety - 1 | BlockStructure: Multiple of an IdentityVarietyGolden Rain VarietyMarvellous VarietyVictoryStdDev: 11 11 11ResidualStdDev: 12.9Fixed effects: yield ~ nitroValue Std.Error DF t-value p-value(Intercept) 81.9 6.95 65 11.8 0nitro 73.7 6.78 65 10.9 0Correlation:(Intr)nitro -0.293Standardized Within-Group Residuals:Min Q1 Med Q3 Max-1.7438 -0.6648 0.0171 0.5430 1.8030Number of Observations: 72Number of Groups: 6> ## establishing the desired parameterization for contrasts> options(contrasts = c("contr.treatment", "contr.poly"))> fm1Assay <- lme(logDens ~ sample * dilut, Assay,+ random = pdBlocked(list(pdIdent(~ 1), pdIdent(~ sample - 1),+ pdIdent(~ dilut - 1))))> fm1AssayLinear mixed-effects model fit by REMLData: AssayLog-restricted-likelihood: 38.5Fixed: logDens ~ sample * dilut(Intercept) sampleb samplec sampled-0.18279 0.08075 0.13398 0.20770samplee samplef dilut2 dilut3-0.02367 0.07357 0.20443 0.40586dilut4 dilut5 sampleb:dilut2 samplec:dilut20.57319 0.72064 0.00894 -0.00850sampled:dilut2 samplee:dilut2 samplef:dilut2 sampleb:dilut30.00108 -0.04192 0.01935 -0.02507samplec:dilut3 sampled:dilut3 samplee:dilut3 samplef:dilut30.01865 0.00399 -0.02771 0.05432sampleb:dilut4 samplec:dilut4 sampled:dilut4 samplee:dilut40.06079 0.00526 -0.01649 0.04980samplef:dilut4 sampleb:dilut5 samplec:dilut5 sampled:dilut50.06337 -0.04576 -0.07260 -0.17776samplee:dilut5 samplef:dilut50.01361 0.00402Random effects:Composite Structure: BlockedBlock 1: (Intercept)Formula: ~1 | Block(Intercept)StdDev: 0.00981Block 2: samplea, sampleb, samplec, sampled, samplee, samplefFormula: ~sample - 1 | BlockStructure: Multiple of an Identitysamplea sampleb samplec sampled samplee samplefStdDev: 0.0253 0.0253 0.0253 0.0253 0.0253 0.0253Block 3: dilut1, dilut2, dilut3, dilut4, dilut5Formula: ~dilut - 1 | BlockStructure: Multiple of an Identitydilut1 dilut2 dilut3 dilut4 dilut5 ResidualStdDev: 0.00913 0.00913 0.00913 0.00913 0.00913 0.0416Number of Observations: 60Number of Groups: 2> anova(fm1Assay)numDF denDF F-value p-value(Intercept) 1 29 538 <.0001sample 5 29 11 <.0001dilut 4 29 421 <.0001sample:dilut 20 29 2 0.119> formula(Oxide)Thickness ~ 1 | Lot/Wafer> fm1Oxide <- lme(Thickness ~ 1, Oxide)> fm1OxideLinear mixed-effects model fit by REMLData: OxideLog-restricted-likelihood: -227Fixed: Thickness ~ 1(Intercept)2000Random effects:Formula: ~1 | Lot(Intercept)StdDev: 11.4Formula: ~1 | Wafer %in% Lot(Intercept) ResidualStdDev: 5.99 3.55Number of Observations: 72Number of Groups:Lot Wafer %in% Lot8 24> intervals(fm1Oxide, which = "var-cov")Approximate 95% confidence intervalsRandom Effects:Level: Lotlower est. uppersd((Intercept)) 6.39 11.4 20.3Level: Waferlower est. uppersd((Intercept)) 4.06 5.99 8.82Within-group standard error:lower est. upper2.90 3.55 4.33> fm2Oxide <- update(fm1Oxide, random = ~ 1 | Lot)> anova(fm1Oxide, fm2Oxide)Model df AIC BIC logLik Test L.Ratio p-valuefm1Oxide 1 4 462 471 -227fm2Oxide 2 3 497 504 -246 1 vs 2 37.1 <.0001> coef(fm1Oxide, level = 1)(Intercept)1 19972 19893 20014 19965 20146 20207 19928 1994> coef(fm1Oxide, level = 2)(Intercept)1/1 20031/2 19851/3 20012/1 19902/2 19882/3 19863/1 20023/2 20003/3 20004/1 19964/2 19994/3 19915/1 20095/2 20175/3 20196/1 20316/2 20226/3 20117/1 19907/2 19917/3 19928/1 19948/2 19958/3 1991> ranef(fm1Oxide, level = 1:2)Level: Lot(Intercept)1 -3.4632 -11.2223 0.8694 -4.4715 13.4636 19.4087 -8.1998 -6.385Level: Wafer %in% Lot(Intercept)1/1 6.54601/2 -11.95891/3 4.45672/1 0.65862/2 -0.83372/3 -2.92303/1 1.47283/2 -0.61643/3 -0.61644/1 -0.01354/2 3.26964/3 -4.49055/1 -4.43185/2 3.02985/3 5.11916/1 11.73506/2 2.18416/3 -8.56077/1 -1.74947/2 -0.55567/3 0.04148/1 -0.09028/2 1.40218/3 -3.0749> fm1Wafer <- lme(current ~ voltage + I(voltage^2), data = Wafer,+ random = list(Wafer = pdDiag(~voltage + I(voltage^2)),+ Site = pdDiag(~voltage + I(voltage^2))))> ## IGNORE_RDIFF_BEGIN> summary(fm1Wafer)Linear mixed-effects model fit by REMLData: WaferAIC BIC logLik-282 -242 151Random effects:Formula: ~voltage + I(voltage^2) | WaferStructure: Diagonal(Intercept) voltage I(voltage^2)StdDev: 2.81e-05 0.187 0.025Formula: ~voltage + I(voltage^2) | Site %in% WaferStructure: Diagonal(Intercept) voltage I(voltage^2) ResidualStdDev: 8.17e-06 0.136 2.45e-08 0.115Fixed effects: current ~ voltage + I(voltage^2)Value Std.Error DF t-value p-value(Intercept) -4.46 0.0513 318 -87.0 0voltage 5.90 0.0927 318 63.7 0I(voltage^2) 1.17 0.0230 318 51.0 0Correlation:(Intr) voltagvoltage -0.735I(voltage^2) 0.884 -0.698Standardized Within-Group Residuals:Min Q1 Med Q3 Max-1.8966 -0.5354 0.0249 0.7985 1.7777Number of Observations: 400Number of Groups:Wafer Site %in% Wafer10 80> ## IGNORE_RDIFF_END> fitted(fm1Wafer, level = 0)1 1 1 1 1 1 1 1 1 11.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.451 1 1 1 1 1 1 1 1 11.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.451 1 1 1 1 1 1 1 1 11.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.451 1 1 1 1 1 1 1 1 11.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.452 2 2 2 2 2 2 2 2 21.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.452 2 2 2 2 2 2 2 2 21.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.452 2 2 2 2 2 2 2 2 21.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.452 2 2 2 2 2 2 2 2 21.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.453 3 3 3 3 3 3 3 3 31.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.453 3 3 3 3 3 3 3 3 31.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.453 3 3 3 3 3 3 3 3 31.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.453 3 3 3 3 3 3 3 3 31.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.454 4 4 4 4 4 4 4 4 41.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.454 4 4 4 4 4 4 4 4 41.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.454 4 4 4 4 4 4 4 4 41.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.454 4 4 4 4 4 4 4 4 41.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.455 5 5 5 5 5 5 5 5 51.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.455 5 5 5 5 5 5 5 5 51.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.455 5 5 5 5 5 5 5 5 51.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.455 5 5 5 5 5 5 5 5 51.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.456 6 6 6 6 6 6 6 6 61.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.456 6 6 6 6 6 6 6 6 61.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.456 6 6 6 6 6 6 6 6 61.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.456 6 6 6 6 6 6 6 6 61.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.457 7 7 7 7 7 7 7 7 71.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.457 7 7 7 7 7 7 7 7 71.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.457 7 7 7 7 7 7 7 7 71.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.457 7 7 7 7 7 7 7 7 71.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.458 8 8 8 8 8 8 8 8 81.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.458 8 8 8 8 8 8 8 8 81.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.458 8 8 8 8 8 8 8 8 81.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.458 8 8 8 8 8 8 8 8 81.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.459 9 9 9 9 9 9 9 9 91.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.459 9 9 9 9 9 9 9 9 91.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.459 9 9 9 9 9 9 9 9 91.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.459 9 9 9 9 9 9 9 9 91.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.4510 10 10 10 10 10 10 10 10 101.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.4510 10 10 10 10 10 10 10 10 101.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.4510 10 10 10 10 10 10 10 10 101.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.4510 10 10 10 10 10 10 10 10 101.01 4.31 7.98 12.03 16.45 1.01 4.31 7.98 12.03 16.45attr(,"label")[1] "Fitted values (mA)"> resid(fm1Wafer, level = 1:2)Wafer Site1 0.061492 0.0680622 -0.189869 -0.1800133 -0.015086 -0.0019444 0.103762 0.1201895 -0.053726 -0.0340146 0.192612 0.0737367 0.044131 -0.1341838 0.275914 0.0381639 0.431762 0.13457310 0.306274 -0.05035311 0.084612 0.06017712 -0.150069 -0.18672213 0.045314 -0.00355614 0.185762 0.12467515 0.054274 -0.01903116 0.042212 0.07367117 -0.237069 -0.18988018 -0.077086 -0.01416719 0.035762 0.11441020 -0.121726 -0.02734821 0.092692 0.07669622 -0.149069 -0.17306223 0.033314 0.00132424 0.159762 0.11977425 0.014274 -0.03371226 -0.057768 0.11184127 -0.453669 -0.19925628 -0.379286 -0.04006829 -0.339838 0.08418530 -0.553726 -0.04489931 0.047012 0.08804832 -0.238069 -0.17651433 -0.090886 -0.00881234 0.007762 0.11035435 -0.165726 -0.04261636 0.079392 0.08484137 -0.172269 -0.16409538 -0.006486 0.00441339 0.101762 0.11538540 -0.063726 -0.04737841 0.038702 0.06547642 -0.209573 -0.16941243 -0.048411 0.00513744 0.036789 0.10372445 -0.117373 -0.03705246 0.266102 0.15185247 0.114027 -0.05734948 0.302789 0.07428849 0.390789 0.10516350 0.226627 -0.11612551 0.299502 0.20513552 0.129627 -0.01192553 0.280989 0.09225454 0.324789 0.08887055 0.120627 -0.16247756 -0.032838 0.07244957 -0.343973 -0.18604358 -0.225411 -0.01483859 -0.171211 0.09200560 -0.353373 -0.03751461 0.262902 0.16078662 0.095827 -0.05734863 0.274189 0.06995664 0.356789 0.10149765 0.188627 -0.11772466 0.000342 0.08786767 -0.298573 -0.16728668 -0.178211 -0.00316169 -0.127211 0.09160170 -0.315373 -0.05279971 0.100502 0.12728572 -0.153973 -0.11380073 -0.025611 0.02795474 0.022789 0.08974575 -0.169373 -0.08902776 0.032102 0.09718677 -0.244373 -0.14674878 -0.120611 0.00955679 -0.071211 0.09149880 -0.261373 -0.06612381 -0.004099 0.05271782 -0.278076 -0.19285283 -0.127696 -0.01406484 -0.029418 0.11262285 -0.197444 -0.02699586 0.052321 0.08924987 -0.208276 -0.15288488 -0.067296 0.00656189 0.014582 0.10690290 -0.171444 -0.06065991 0.118641 0.06278292 -0.064476 -0.14826493 0.134904 0.02318794 0.266582 0.12693595 0.120556 -0.04701996 -0.041079 0.05126597 -0.346676 -0.20816098 -0.212496 -0.02780799 -0.121418 0.109442100 -0.297444 -0.020411101 0.128041 0.079868102 -0.066076 -0.138336103 0.121104 0.024758104 0.240582 0.120149105 0.088556 -0.055963106 -0.091839 0.070304107 -0.452476 -0.209261108 -0.361896 -0.037608109 -0.311418 0.093941110 -0.519444 -0.033013111 0.286041 0.146703112 0.154524 -0.054483113 0.353704 0.075028114 0.468582 0.120237115 0.298556 -0.119458116 0.253641 0.183845117 0.066124 -0.038569118 0.211104 0.071514119 0.274582 0.100094120 0.062556 -0.146829121 0.113168 0.059522122 -0.082704 -0.163173123 0.123749 0.016457124 0.262907 0.128791125 0.124569 -0.036370126 0.199348 0.075597127 0.057096 -0.128531128 0.288549 0.041047129 0.444907 0.135529130 0.316569 -0.054685131 0.010568 0.105606132 -0.309104 -0.166546133 -0.198251 -0.008174134 -0.139093 0.098502135 -0.349431 -0.064316136 0.000368 0.076116137 -0.314704 -0.201082138 -0.178051 -0.026555139 -0.083093 0.106277140 -0.251431 -0.024187141 0.016268 0.116152142 -0.315904 -0.166078143 -0.212251 -0.012483144 -0.155093 0.094617145 -0.363431 -0.063779146 0.004348 0.054446147 -0.286504 -0.211357148 -0.125651 -0.025456149 -0.009093 0.116151150 -0.161431 -0.011138151 0.096848 0.080552152 -0.138304 -0.162749153 0.039349 0.006756154 0.158907 0.118165155 0.006569 -0.042321156 0.118788 0.080347157 -0.096904 -0.154565158 0.090949 0.014067159 0.218907 0.122804160 0.068569 -0.046754161 -0.029651 0.042434162 -0.299821 -0.191694163 -0.157165 -0.012996164 -0.067684 0.112527165 -0.246778 -0.030524166 0.116949 0.129128167 -0.114421 -0.096153168 0.013235 0.037592169 0.072316 0.102762170 -0.146778 -0.110243171 0.197149 0.101805172 0.049179 -0.093837173 0.245235 0.054547174 0.362316 0.123955175 0.195222 -0.090810176 0.048749 0.058063177 -0.177021 -0.163051178 -0.010365 0.008261179 0.094316 0.117599180 -0.072778 -0.044839181 0.214149 0.102694182 0.073779 -0.093404183 0.277635 0.054723184 0.402316 0.123676185 0.249222 -0.085146186 -0.092031 0.056118187 -0.426021 -0.203798188 -0.326965 -0.030668189 -0.271684 0.098687190 -0.478778 -0.034333191 0.187949 0.129497192 0.004979 -0.082699193 0.168635 0.051730194 0.256316 0.110185195 0.069222 -0.106135196 0.095349 0.120157197 -0.145621 -0.108410198 -0.019765 0.029849199 0.040316 0.102333200 -0.174778 -0.100357201 0.115311 0.075105202 -0.097595 -0.157904203 0.094646 0.014234204 0.223635 0.123120205 0.077572 -0.043047206 0.121051 0.100980207 -0.110195 -0.140302208 0.058846 0.018704209 0.165635 0.115457210 -0.004428 -0.064642211 0.079591 0.081229212 -0.172795 -0.170338213 -0.001954 0.001323214 0.113635 0.117730215 -0.046428 -0.041514216 0.007011 0.076714217 -0.304595 -0.200040218 -0.163954 -0.024547219 -0.066365 0.107893220 -0.234428 -0.025319221 0.066991 0.085934222 -0.202595 -0.174180223 -0.042754 -0.004867224 0.065635 0.112994225 -0.096428 -0.039598226 -0.020549 0.093776227 -0.371395 -0.199907228 -0.261754 -0.033103229 -0.188365 0.097449230 -0.376428 -0.033452231 0.124251 0.092105232 -0.097195 -0.145414233 0.081646 0.017355234 0.199635 0.119271235 0.039572 -0.056865236 0.104871 0.083043237 -0.123995 -0.156738238 0.055046 0.011389239 0.173635 0.119064240 0.017572 -0.047914241 0.227356 0.097058242 0.136724 -0.058724243 0.348539 0.087942244 0.457002 0.131256245 0.268913 -0.121982246 -0.049644 -0.001886247 -0.250476 -0.178840248 -0.082661 0.012853249 0.007002 0.126395250 -0.185087 -0.041815251 0.491556 0.164445252 0.535924 0.045257253 0.814739 0.160517254 0.963002 0.145224255 0.798913 -0.182420256 0.035556 -0.000644257 -0.106476 -0.160777258 0.103139 0.030738259 0.229002 0.138501260 0.066913 -0.041688261 0.084356 0.047445262 -0.064476 -0.119844263 0.122139 0.048316264 0.219002 0.126723265 0.030913 -0.079821266 -0.102844 -0.008156267 -0.348676 -0.206645268 -0.197861 -0.008486269 -0.112998 0.123721270 -0.311087 -0.027023271 -0.104044 0.032614272 -0.381676 -0.176689273 -0.275861 -0.002545274 -0.234998 0.106647275 -0.471087 -0.061112276 -0.127844 0.030339277 -0.422076 -0.184802278 -0.325461 -0.009095279 -0.292998 0.102459280 -0.531087 -0.056538281 0.272748 0.047840282 0.262060 -0.075302283 0.546385 0.096569284 0.718724 0.156454285 0.586276 -0.088447286 0.249948 0.062457287 0.206660 -0.074577288 0.464785 0.089803289 0.616724 0.147996290 0.466276 -0.096197291 -0.032652 0.011344292 -0.243540 -0.177547293 -0.075615 0.012376294 0.014724 0.124713295 -0.175724 -0.043737296 -0.108452 -0.012938297 -0.355740 -0.212470298 -0.201215 -0.010187299 -0.113276 0.125508300 -0.309724 -0.023182301 -0.096052 0.018185302 -0.362540 -0.191185303 -0.234015 -0.005541304 -0.171276 0.114316305 -0.387724 -0.045013306 -0.123652 -0.009999307 -0.389340 -0.218861308 -0.245215 -0.017909309 -0.163276 0.120856310 -0.359724 -0.018765311 -0.108452 0.037755312 -0.402940 -0.183630313 -0.300215 -0.007802314 -0.259276 0.106241315 -0.497724 -0.059104316 0.285348 0.119276317 0.217660 -0.031449318 0.435785 0.103641319 0.548724 0.133543320 0.356276 -0.141940321 0.066249 0.061990322 -0.106937 -0.113325323 0.058058 0.049541324 0.133235 0.122588325 -0.084807 -0.097583326 0.058049 0.013390327 -0.080137 -0.147125328 0.128858 0.039540329 0.251235 0.139588330 0.077193 -0.056784331 0.004649 0.041019332 -0.199737 -0.145182333 -0.044542 0.028198334 0.029235 0.120160335 -0.182807 -0.073697336 0.088449 -0.002738337 -0.008937 -0.145718338 0.230458 0.048083339 0.377235 0.149266340 0.225193 -0.048369341 0.017249 0.032907342 -0.167737 -0.144250343 0.000858 0.032174344 0.085235 0.124380345 -0.116807 -0.069833346 -0.084751 -0.026585347 -0.303337 -0.216088348 -0.124142 -0.007810349 -0.012765 0.132649350 -0.184807 -0.010310351 -0.104351 0.042779352 -0.394337 -0.173642353 -0.296142 -0.001881354 -0.264765 0.103060355 -0.508807 -0.067416356 0.278649 0.082220357 0.247263 -0.047380358 0.498058 0.105201359 0.635235 0.144163360 0.469193 -0.120093361 -0.107404 -0.012451362 -0.354401 -0.211972363 -0.199807 -0.009901364 -0.112820 0.124562365 -0.307642 -0.022784366 0.231396 0.061385367 0.179399 -0.075617368 0.428793 0.088772369 0.571180 0.146153370 0.410358 -0.099674371 0.162596 0.050096372 0.066399 -0.102351373 0.292993 0.067994374 0.421180 0.139931375 0.252358 -0.085141376 -0.022804 -0.020045377 -0.198801 -0.194662378 0.005393 0.010912379 0.131180 0.138078380 -0.027642 -0.019364381 0.015396 0.011892382 -0.160401 -0.165658383 0.030993 0.023985384 0.141180 0.132419385 -0.035642 -0.046155386 -0.105404 0.014675387 -0.373401 -0.193283388 -0.246407 -0.006248389 -0.184820 0.115377390 -0.405642 -0.045405391 0.164196 0.113448392 0.015399 -0.060723393 0.177393 0.075897394 0.243180 0.116310395 0.016358 -0.135886396 -0.008404 0.037007397 -0.216401 -0.148285398 -0.065007 0.025815399 0.005180 0.118707400 -0.207642 -0.071409> newWafer <-+ data.frame(Wafer = rep(1, 4), voltage = c(1, 1.5, 3, 3.5))> predict(fm1Wafer, newWafer, level = 0:1)Wafer predict.fixed predict.Wafer1 1 2.61 2.402 1 7.03 6.723 1 23.78 23.234 1 30.54 29.92> newWafer2 <- data.frame(Wafer = rep(1, 4), Site = rep(3, 4),+ voltage = c(1, 1.5, 3, 3.5))> predict(fm1Wafer, newWafer2, level = 0:2)Wafer Site predict.fixed predict.Wafer predict.Site1 1 1/3 2.61 2.40 2.432 1 1/3 7.03 6.72 6.773 1 1/3 23.78 23.23 23.324 1 1/3 30.54 29.92 30.03> # 4.3 Examining a Fitted Model>> plot(fm2Orth.lme, Subject~resid(.), abline = 0)> plot(fm2Orth.lme, resid(., type = "p") ~ fitted(.) | Sex,+ id = 0.05, adj = -0.3)> fm3Orth.lme <-+ update(fm2Orth.lme, weights = varIdent(form = ~ 1 | Sex))> fm3Orth.lmeLinear mixed-effects model fit by REMLData: OrthodontLog-restricted-likelihood: -206Fixed: distance ~ Sex + I(age - 11) + Sex:I(age - 11)(Intercept) SexFemale24.969 -2.321I(age - 11) SexFemale:I(age - 11)0.784 -0.305Random effects:Formula: ~I(age - 11) | SubjectStructure: General positive-definite, Log-Cholesky parametrizationStdDev Corr(Intercept) 1.855 (Intr)I(age - 11) 0.157 0.394Residual 1.630Variance function:Structure: Different standard deviations per stratumFormula: ~1 | SexParameter estimates:Male Female1.000 0.409Number of Observations: 108Number of Groups: 27> plot(fm3Orth.lme, distance ~ fitted(.),+ id = 0.05, adj = -0.3)> anova(fm2Orth.lme, fm3Orth.lme)Model df AIC BIC logLik Test L.Ratio p-valuefm2Orth.lme 1 8 451 473 -218fm3Orth.lme 2 9 430 453 -206 1 vs 2 23.8 <.0001> qqnorm(fm3Orth.lme, ~resid(.) | Sex)> plot(fm2IGF.lme, resid(., type = "p") ~ fitted(.) | Lot,+ layout = c(5,2))> qqnorm(fm2IGF.lme, ~ resid(.), id = 0.05, adj = -0.75)> plot(fm1Oxide)> qqnorm(fm1Oxide)> plot(fm1Wafer, resid(.) ~ voltage | Wafer)> plot(fm1Wafer, resid(.) ~ voltage | Wafer,+ panel = function(x, y, ...) {+ panel.grid()+ panel.xyplot(x, y)+ panel.loess(x, y, lty = 2)+ panel.abline(0, 0)+ })> with(Wafer,+ coef(lm(resid(fm1Wafer) ~ cos(4.19*voltage)+sin(4.19*voltage)-1)))cos(4.19 * voltage) sin(4.19 * voltage)-0.0519 0.1304> nls(resid(fm1Wafer) ~ b3*cos(w*voltage) + b4*sin(w*voltage), Wafer,+ start = list(b3 = -0.0519, b4 = 0.1304, w = 4.19))Nonlinear regression modelmodel: resid(fm1Wafer) ~ b3 * cos(w * voltage) + b4 * sin(w * voltage)data: Waferb3 b4 w-0.1117 0.0777 4.5679residual sum-of-squares: 0.729Number of iterations to convergence: 6Achieved convergence tolerance: 1.12e-06> fm2Wafer <- update(fm1Wafer,+ . ~ . + cos(4.5679*voltage) + sin(4.5679*voltage),+ random = list(Wafer=pdDiag(~voltage+I(voltage^2)),+ Site=pdDiag(~voltage+I(voltage^2))))> summary(fm2Wafer)Linear mixed-effects model fit by REMLData: WaferAIC BIC logLik-1233 -1185 628Random effects:Formula: ~voltage + I(voltage^2) | WaferStructure: Diagonal(Intercept) voltage I(voltage^2)StdDev: 0.129 0.349 0.0491Formula: ~voltage + I(voltage^2) | Site %in% WaferStructure: Diagonal(Intercept) voltage I(voltage^2) ResidualStdDev: 0.0397 0.234 0.0475 0.0113Fixed effects: current ~ voltage + I(voltage^2) + cos(4.5679 * voltage) + sin(4.5679 * voltage)Value Std.Error DF t-value p-value(Intercept) -4.26 0.0422 316 -100.8 0voltage 5.62 0.1142 316 49.2 0I(voltage^2) 1.26 0.0170 316 74.2 0cos(4.5679 * voltage) -0.10 0.0011 316 -85.0 0sin(4.5679 * voltage) 0.10 0.0015 316 69.4 0Correlation:(Intr) voltag I(v^2) c(4.*vvoltage -0.029I(voltage^2) 0.060 -0.031cos(4.5679 * voltage) 0.162 -0.082 0.172sin(4.5679 * voltage) 0.200 -0.101 0.212 0.567Standardized Within-Group Residuals:Min Q1 Med Q3 Max-2.4272 -0.4032 0.0253 0.3936 2.8427Number of Observations: 400Number of Groups:Wafer Site %in% Wafer10 80> ## IGNORE_RDIFF_BEGIN> intervals(fm2Wafer)Approximate 95% confidence intervalsFixed effects:lower est. upper(Intercept) -4.3385 -4.2554 -4.1723voltage 5.3977 5.6224 5.8470I(voltage^2) 1.2251 1.2585 1.2919cos(4.5679 * voltage) -0.0978 -0.0956 -0.0933sin(4.5679 * voltage) 0.1014 0.1043 0.1073Random Effects:Level: Waferlower est. uppersd((Intercept)) 0.0802 0.1289 0.207sd(voltage) 0.2135 0.3487 0.569sd(I(voltage^2)) 0.0290 0.0491 0.083Level: Sitelower est. uppersd((Intercept)) 0.0220 0.0397 0.0717sd(voltage) 0.1909 0.2344 0.2878sd(I(voltage^2)) 0.0383 0.0475 0.0590Within-group standard error:lower est. upper0.00927 0.01133 0.01383> ## IGNORE_RDIFF_END> qqnorm(fm2Wafer)> qqnorm(fm2Orth.lme, ~ranef(.), id = 0.10, cex = 0.7)> pairs(fm2Orth.lme, ~ranef(.) | Sex,+ id = ~ Subject == "M13", adj = -0.3)> fm2IGF.lmeLinear mixed-effects model fit by REMLData: IGFLog-restricted-likelihood: -297Fixed: conc ~ age(Intercept) age5.36904 -0.00193Random effects:Formula: ~age | LotStructure: Diagonal(Intercept) age ResidualStdDev: 3.62e-05 0.00537 0.822Number of Observations: 237Number of Groups: 10> c(0.00031074, 0.0053722)/abs(fixef(fm2IGF.lme))(Intercept) age5.79e-05 2.78e+00> fm3IGF.lme <- update(fm2IGF.lme, random = ~ age - 1)> anova(fm2IGF.lme, fm3IGF.lme)Model df AIC BIC logLik Test L.Ratio p-valuefm2IGF.lme 1 5 605 622 -297fm3IGF.lme 2 4 603 617 -297 1 vs 2 1.47e-07 1> qqnorm(fm1Oxide, ~ranef(., level = 1), id=0.10)> qqnorm(fm1Oxide, ~ranef(., level = 2), id=0.10)> #fm3Wafer <- update(fm2Wafer,> # random = list(Wafer = ~voltage+I(voltage^2),> # Site = pdDiag(~voltage+I(voltage^2))),> # control = list(msVerbose = TRUE, msMaxIter = 200)> # )> #fm3Wafer> #anova(fm2Wafer, fm3Wafer)> #fm4Wafer <- update(fm2Wafer,> # random = list(Wafer = ~ voltage + I(voltage^2),> # Site = pdBlocked(list(~1,> # ~voltage+I(voltage^2) - 1))),> # control = list(msVerbose = TRUE,> # msMaxIter = 200))> #fm4Wafer> #anova(fm3Wafer, fm4Wafer)> #qqnorm(fm4Wafer, ~ranef(., level = 2), id = 0.05,> # cex = 0.7, layout = c(3, 1))>> # The next line is not in the book but is needed to get fm1Machine>> fm1Machine <-+ lme(score ~ Machine, data = Machines, random = ~ 1 | Worker)> (fm3Machine <- update(fm1Machine, random = ~Machine-1|Worker))Linear mixed-effects model fit by REMLData: MachinesLog-restricted-likelihood: -104Fixed: score ~ Machine(Intercept) MachineB MachineC52.36 7.97 13.92Random effects:Formula: ~Machine - 1 | WorkerStructure: General positive-definite, Log-Cholesky parametrizationStdDev CorrMachineA 4.079 MachnA MachnBMachineB 8.625 0.803MachineC 4.389 0.623 0.771Residual 0.962Number of Observations: 54Number of Groups: 6> # cleanup>> summary(warnings())No warnings======ch05.R======> #-*- R -*->> # initialization>> library(nlme)> options(width = 65,+ ## reduce platform dependence in printed output when testing+ digits = if(nzchar(Sys.getenv("R_TESTS"))) 3 else 5)> options(contrasts = c(unordered = "contr.helmert", ordered = "contr.poly"))> pdf(file = "ch05.pdf")> # Chapter 5 Extending the Basic Linear Mixed-Effects Models>> # 5.1 General Formulation of the Extended Model>> vf1Fixed <- varFixed(~ age)> vf1Fixed <- Initialize(vf1Fixed, data = Orthodont)> varWeights(vf1Fixed)[1] 0.354 0.316 0.289 0.267 0.354 0.316 0.289 0.267 0.354 0.316[11] 0.289 0.267 0.354 0.316 0.289 0.267 0.354 0.316 0.289 0.267[21] 0.354 0.316 0.289 0.267 0.354 0.316 0.289 0.267 0.354 0.316[31] 0.289 0.267 0.354 0.316 0.289 0.267 0.354 0.316 0.289 0.267[41] 0.354 0.316 0.289 0.267 0.354 0.316 0.289 0.267 0.354 0.316[51] 0.289 0.267 0.354 0.316 0.289 0.267 0.354 0.316 0.289 0.267[61] 0.354 0.316 0.289 0.267 0.354 0.316 0.289 0.267 0.354 0.316[71] 0.289 0.267 0.354 0.316 0.289 0.267 0.354 0.316 0.289 0.267[81] 0.354 0.316 0.289 0.267 0.354 0.316 0.289 0.267 0.354 0.316[91] 0.289 0.267 0.354 0.316 0.289 0.267 0.354 0.316 0.289 0.267[101] 0.354 0.316 0.289 0.267 0.354 0.316 0.289 0.267> vf1Ident <- varIdent(c(Female = 0.5), ~ 1 | Sex)> vf1Ident <- Initialize(vf1Ident, Orthodont)> varWeights(vf1Ident)Male Male Male Male Male Male Male Male Male1 1 1 1 1 1 1 1 1Male Male Male Male Male Male Male Male Male1 1 1 1 1 1 1 1 1Male Male Male Male Male Male Male Male Male1 1 1 1 1 1 1 1 1Male Male Male Male Male Male Male Male Male1 1 1 1 1 1 1 1 1Male Male Male Male Male Male Male Male Male1 1 1 1 1 1 1 1 1Male Male Male Male Male Male Male Male Male1 1 1 1 1 1 1 1 1Male Male Male Male Male Male Male Male Male1 1 1 1 1 1 1 1 1Male Female Female Female Female Female Female Female Female1 2 2 2 2 2 2 2 2Female Female Female Female Female Female Female Female Female2 2 2 2 2 2 2 2 2Female Female Female Female Female Female Female Female Female2 2 2 2 2 2 2 2 2Female Female Female Female Female Female Female Female Female2 2 2 2 2 2 2 2 2Female Female Female Female Female Female Female Female Female2 2 2 2 2 2 2 2 2> vf2Ident <- varIdent(form = ~ 1 | Sex, fixed = c(Female = 0.5))> vf2Ident <- Initialize(vf2Ident, Orthodont)> varWeights(vf2Ident)Male Male Male Male Male Male Male Male Male1 1 1 1 1 1 1 1 1Male Male Male Male Male Male Male Male Male1 1 1 1 1 1 1 1 1Male Male Male Male Male Male Male Male Male1 1 1 1 1 1 1 1 1Male Male Male Male Male Male Male Male Male1 1 1 1 1 1 1 1 1Male Male Male Male Male Male Male Male Male1 1 1 1 1 1 1 1 1Male Male Male Male Male Male Male Male Male1 1 1 1 1 1 1 1 1Male Male Male Male Male Male Male Male Male1 1 1 1 1 1 1 1 1Male Female Female Female Female Female Female Female Female1 2 2 2 2 2 2 2 2Female Female Female Female Female Female Female Female Female2 2 2 2 2 2 2 2 2Female Female Female Female Female Female Female Female Female2 2 2 2 2 2 2 2 2Female Female Female Female Female Female Female Female Female2 2 2 2 2 2 2 2 2Female Female Female Female Female Female Female Female Female2 2 2 2 2 2 2 2 2> vf3Ident <- varIdent(form = ~ 1 | Sex * age)> vf3Ident <- Initialize(vf3Ident, Orthodont)> varWeights(vf3Ident)Male*8 Male*10 Male*12 Male*14 Male*8 Male*101 1 1 1 1 1Male*12 Male*14 Male*8 Male*10 Male*12 Male*141 1 1 1 1 1Male*8 Male*10 Male*12 Male*14 Male*8 Male*101 1 1 1 1 1Male*12 Male*14 Male*8 Male*10 Male*12 Male*141 1 1 1 1 1Male*8 Male*10 Male*12 Male*14 Male*8 Male*101 1 1 1 1 1Male*12 Male*14 Male*8 Male*10 Male*12 Male*141 1 1 1 1 1Male*8 Male*10 Male*12 Male*14 Male*8 Male*101 1 1 1 1 1Male*12 Male*14 Male*8 Male*10 Male*12 Male*141 1 1 1 1 1Male*8 Male*10 Male*12 Male*14 Male*8 Male*101 1 1 1 1 1Male*12 Male*14 Male*8 Male*10 Male*12 Male*141 1 1 1 1 1Male*8 Male*10 Male*12 Male*14 Female*8 Female*101 1 1 1 1 1Female*12 Female*14 Female*8 Female*10 Female*12 Female*141 1 1 1 1 1Female*8 Female*10 Female*12 Female*14 Female*8 Female*101 1 1 1 1 1Female*12 Female*14 Female*8 Female*10 Female*12 Female*141 1 1 1 1 1Female*8 Female*10 Female*12 Female*14 Female*8 Female*101 1 1 1 1 1Female*12 Female*14 Female*8 Female*10 Female*12 Female*141 1 1 1 1 1Female*8 Female*10 Female*12 Female*14 Female*8 Female*101 1 1 1 1 1Female*12 Female*14 Female*8 Female*10 Female*12 Female*141 1 1 1 1 1> vf1Power <- varPower(1)> formula(vf1Power)~fitted(.)<environment: 0x55aa58bbc708>> vf2Power <- varPower(fixed = 0.5)> vf3Power <- varPower(form = ~ fitted(.) | Sex,+ fixed = list(Male = 0.5, Female = 0))> vf1Exp <- varExp(form = ~ age | Sex, fixed = c(Female = 0))> vf1ConstPower <- varConstPower(power = 0.5,+ fixed = list(const = 1))> vf1Comb <- varComb(varIdent(c(Female = 0.5), ~ 1 | Sex),+ varExp(1, ~ age))> vf1Comb <- Initialize(vf1Comb, Orthodont)> varWeights(vf1Comb)[1] 3.35e-04 4.54e-05 6.14e-06 8.32e-07 3.35e-04 4.54e-05[7] 6.14e-06 8.32e-07 3.35e-04 4.54e-05 6.14e-06 8.32e-07[13] 3.35e-04 4.54e-05 6.14e-06 8.32e-07 3.35e-04 4.54e-05[19] 6.14e-06 8.32e-07 3.35e-04 4.54e-05 6.14e-06 8.32e-07[25] 3.35e-04 4.54e-05 6.14e-06 8.32e-07 3.35e-04 4.54e-05[31] 6.14e-06 8.32e-07 3.35e-04 4.54e-05 6.14e-06 8.32e-07[37] 3.35e-04 4.54e-05 6.14e-06 8.32e-07 3.35e-04 4.54e-05[43] 6.14e-06 8.32e-07 3.35e-04 4.54e-05 6.14e-06 8.32e-07[49] 3.35e-04 4.54e-05 6.14e-06 8.32e-07 3.35e-04 4.54e-05[55] 6.14e-06 8.32e-07 3.35e-04 4.54e-05 6.14e-06 8.32e-07[61] 3.35e-04 4.54e-05 6.14e-06 8.32e-07 6.71e-04 9.08e-05[67] 1.23e-05 1.66e-06 6.71e-04 9.08e-05 1.23e-05 1.66e-06[73] 6.71e-04 9.08e-05 1.23e-05 1.66e-06 6.71e-04 9.08e-05[79] 1.23e-05 1.66e-06 6.71e-04 9.08e-05 1.23e-05 1.66e-06[85] 6.71e-04 9.08e-05 1.23e-05 1.66e-06 6.71e-04 9.08e-05[91] 1.23e-05 1.66e-06 6.71e-04 9.08e-05 1.23e-05 1.66e-06[97] 6.71e-04 9.08e-05 1.23e-05 1.66e-06 6.71e-04 9.08e-05[103] 1.23e-05 1.66e-06 6.71e-04 9.08e-05 1.23e-05 1.66e-06> fm1Dial.lme <-+ lme(rate ~(pressure + I(pressure^2) + I(pressure^3) + I(pressure^4))*QB,+ Dialyzer, ~ pressure + I(pressure^2))> fm1Dial.lmeLinear mixed-effects model fit by REMLData: DialyzerLog-restricted-likelihood: -326Fixed: rate ~ (pressure + I(pressure^2) + I(pressure^3) + I(pressure^4)) * QB(Intercept) pressure I(pressure^2)-16.5980 88.6733 -42.7320I(pressure^3) I(pressure^4) QB19.2165 -0.7756 -0.6317pressure:QB1 I(pressure^2):QB1 I(pressure^3):QB10.3104 1.5742 0.0509I(pressure^4):QB1-0.0860Random effects:Formula: ~pressure + I(pressure^2) | SubjectStructure: General positive-definite, Log-Cholesky parametrizationStdDev Corr(Intercept) 1.50 (Intr) pressrpressure 4.91 -0.507I(pressure^2) 1.47 0.311 -0.944Residual 1.82Number of Observations: 140Number of Groups: 20> plot(fm1Dial.lme, resid(.) ~ pressure, abline = 0)> fm2Dial.lme <- update(fm1Dial.lme,+ weights = varPower(form = ~ pressure))> fm2Dial.lmeLinear mixed-effects model fit by REMLData: DialyzerLog-restricted-likelihood: -310Fixed: rate ~ (pressure + I(pressure^2) + I(pressure^3) + I(pressure^4)) * QB(Intercept) pressure I(pressure^2)-17.680 93.711 -49.187I(pressure^3) I(pressure^4) QB112.245 -1.243 -0.921pressure:QB1 I(pressure^2):QB1 I(pressure^3):QB11.353 0.480 0.491I(pressure^4):QB1-0.146Random effects:Formula: ~pressure + I(pressure^2) | SubjectStructure: General positive-definite, Log-Cholesky parametrizationStdDev Corr(Intercept) 1.86 (Intr) pressrpressure 5.33 -0.522I(pressure^2) 1.65 0.362 -0.954Residual 1.26Variance function:Structure: Power of variance covariateFormula: ~pressureParameter estimates:power0.749Number of Observations: 140Number of Groups: 20> anova(fm1Dial.lme, fm2Dial.lme)Model df AIC BIC logLik Test L.Ratio p-valuefm1Dial.lme 1 17 687 736 -326fm2Dial.lme 2 18 655 707 -310 1 vs 2 33.8 <.0001> plot(fm2Dial.lme, resid(., type = "p") ~ pressure,+ abline = 0)> ## IGNORE_RDIFF_BEGIN> intervals(fm2Dial.lme)Approximate 95% confidence intervalsFixed effects:lower est. upper(Intercept) -19.148 -17.680 -16.212pressure 87.231 93.711 100.192I(pressure^2) -57.616 -49.187 -40.757I(pressure^3) 7.967 12.245 16.523I(pressure^4) -1.953 -1.243 -0.533QB1 -2.478 -0.921 0.636pressure:QB1 -5.127 1.353 7.833I(pressure^2):QB1 -7.949 0.480 8.910I(pressure^3):QB1 -3.787 0.491 4.769I(pressure^4):QB1 -0.856 -0.146 0.564Random Effects:Level: Subjectlower est. uppersd((Intercept)) 1.256 1.857 2.7466sd(pressure) 3.623 5.328 7.8363sd(I(pressure^2)) 1.091 1.648 2.4909cor((Intercept),pressure) -0.803 -0.522 -0.0525cor((Intercept),I(pressure^2)) -0.166 0.362 0.7292cor(pressure,I(pressure^2)) -0.985 -0.954 -0.8624Variance function:lower est. upperpower 0.508 0.749 0.991Within-group standard error:lower est. upper1.06 1.26 1.50> ## IGNORE_RDIFF_END> plot(fm2Dial.lme, resid(.) ~ pressure|QB, abline = 0)> fm3Dial.lme <- update(fm2Dial.lme,+ weights=varPower(form = ~ pressure | QB))> fm3Dial.lmeLinear mixed-effects model fit by REMLData: DialyzerLog-restricted-likelihood: -309Fixed: rate ~ (pressure + I(pressure^2) + I(pressure^3) + I(pressure^4)) * QB(Intercept) pressure I(pressure^2)-17.695 93.759 -49.231I(pressure^3) I(pressure^4) QB112.260 -1.244 -1.017pressure:QB1 I(pressure^2):QB1 I(pressure^3):QB11.840 -0.194 0.827I(pressure^4):QB1-0.200Random effects:Formula: ~pressure + I(pressure^2) | SubjectStructure: General positive-definite, Log-Cholesky parametrizationStdDev Corr(Intercept) 1.82 (Intr) pressrpressure 5.24 -0.502I(pressure^2) 1.64 0.338 -0.951Residual 1.26Variance function:Structure: Power of variance covariate, different strataFormula: ~pressure | QBParameter estimates:200 3000.648 0.838Number of Observations: 140Number of Groups: 20> anova(fm2Dial.lme, fm3Dial.lme)Model df AIC BIC logLik Test L.Ratio p-valuefm2Dial.lme 1 18 655 707 -310fm3Dial.lme 2 19 656 711 -309 1 vs 2 0.711 0.399> fm4Dial.lme <- update(fm2Dial.lme,+ weights = varConstPower(form = ~ pressure))> anova(fm2Dial.lme, fm4Dial.lme)Model df AIC BIC logLik Test L.Ratio p-valuefm2Dial.lme 1 18 655 707 -310fm4Dial.lme 2 19 657 711 -309 1 vs 2 0.159 0.69> plot(augPred(fm2Dial.lme), grid = TRUE)> anova(fm2Dial.lme)numDF denDF F-value p-value(Intercept) 1 112 553 <.0001pressure 1 112 2329 <.0001I(pressure^2) 1 112 1175 <.0001I(pressure^3) 1 112 360 <.0001I(pressure^4) 1 112 12 0.0006QB 1 18 5 0.0414pressure:QB 1 112 80 <.0001I(pressure^2):QB 1 112 1 0.2476I(pressure^3):QB 1 112 2 0.1370I(pressure^4):QB 1 112 0 0.6839> anova(fm2Dial.lme, Terms = 8:10)F-test for: I(pressure^2):QB, I(pressure^3):QB, I(pressure^4):QBnumDF denDF F-value p-value1 3 112 1.25 0.294> options(contrasts = c("contr.treatment", "contr.poly"))> fm1BW.lme <- lme(weight ~ Time * Diet, BodyWeight,+ random = ~ Time)> fm1BW.lmeLinear mixed-effects model fit by REMLData: BodyWeightLog-restricted-likelihood: -576Fixed: weight ~ Time * Diet(Intercept) Time Diet2 Diet3 Time:Diet2251.652 0.360 200.665 252.072 0.606Time:Diet30.298Random effects:Formula: ~Time | RatStructure: General positive-definite, Log-Cholesky parametrizationStdDev Corr(Intercept) 36.939 (Intr)Time 0.248 -0.149Residual 4.444Number of Observations: 176Number of Groups: 16> fm2BW.lme <- update(fm1BW.lme, weights = varPower())> fm2BW.lmeLinear mixed-effects model fit by REMLData: BodyWeightLog-restricted-likelihood: -571Fixed: weight ~ Time * Diet(Intercept) Time Diet2 Diet3 Time:Diet2251.602 0.361 200.777 252.170 0.602Time:Diet30.295Random effects:Formula: ~Time | RatStructure: General positive-definite, Log-Cholesky parametrizationStdDev Corr(Intercept) 36.898 (Intr)Time 0.244 -0.145Residual 0.175Variance function:Structure: Power of variance covariateFormula: ~fitted(.)Parameter estimates:power0.543Number of Observations: 176Number of Groups: 16> anova(fm1BW.lme, fm2BW.lme)Model df AIC BIC logLik Test L.Ratio p-valuefm1BW.lme 1 10 1172 1203 -576fm2BW.lme 2 11 1164 1198 -571 1 vs 2 9.8 0.0017> summary(fm2BW.lme)Linear mixed-effects model fit by REMLData: BodyWeightAIC BIC logLik1164 1198 -571Random effects:Formula: ~Time | RatStructure: General positive-definite, Log-Cholesky parametrizationStdDev Corr(Intercept) 36.898 (Intr)Time 0.244 -0.145Residual 0.175Variance function:Structure: Power of variance covariateFormula: ~fitted(.)Parameter estimates:power0.543Fixed effects: weight ~ Time * DietValue Std.Error DF t-value p-value(Intercept) 251.6 13.07 157 19.25 0.0000Time 0.4 0.09 157 4.09 0.0001Diet2 200.8 22.66 13 8.86 0.0000Diet3 252.2 22.66 13 11.13 0.0000Time:Diet2 0.6 0.16 157 3.87 0.0002Time:Diet3 0.3 0.16 157 1.89 0.0601Correlation:(Intr) Time Diet2 Diet3 Tm:Dt2Time -0.152Diet2 -0.577 0.088Diet3 -0.577 0.088 0.333Time:Diet2 0.087 -0.569 -0.157 -0.050Time:Diet3 0.086 -0.567 -0.050 -0.158 0.322Standardized Within-Group Residuals:Min Q1 Med Q3 Max-2.9374 -0.4439 0.0799 0.5808 2.2649Number of Observations: 176Number of Groups: 16> anova(fm2BW.lme, L = c("Time:Diet2" = 1, "Time:Diet3" = -1))F-test for linear combination(s)Time:Diet2 Time:Diet31 -1numDF denDF F-value p-value1 1 157 2.86 0.0926> cs1CompSymm <- corCompSymm(value = 0.3, form = ~ 1 | Subject)> cs2CompSymm <- corCompSymm(value = 0.3, form = ~ age | Subject)> cs1CompSymm <- Initialize(cs1CompSymm, data = Orthodont)> corMatrix(cs1CompSymm)$M01[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M02[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M03[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M04[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M05[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M06[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M07[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M08[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M09[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M10[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M11[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M12[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M13[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M14[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M15[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$M16[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$F01[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$F02[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$F03[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$F04[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$F05[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$F06[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$F07[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$F08[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$F09[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$F10[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0$F11[,1] [,2] [,3] [,4][1,] 1.0 0.3 0.3 0.3[2,] 0.3 1.0 0.3 0.3[3,] 0.3 0.3 1.0 0.3[4,] 0.3 0.3 0.3 1.0> cs1Symm <- corSymm(value = c(0.2, 0.1, -0.1, 0, 0.2, 0),+ form = ~ 1 | Subject)> cs1Symm <- Initialize(cs1Symm, data = Orthodont)> corMatrix(cs1Symm)$M01[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M02[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M03[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M04[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M05[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M06[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M07[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M08[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M09[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M10[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M11[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M12[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M13[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M14[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M15[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$M16[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$F01[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$F02[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$F03[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$F04[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$F05[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$F06[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$F07[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$F08[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$F09[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$F10[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00$F11[,1] [,2] [,3] [,4][1,] 1.0 2.00e-01 1.00e-01 -1.00e-01[2,] 0.2 1.00e+00 9.02e-17 2.00e-01[3,] 0.1 9.02e-17 1.00e+00 -1.04e-16[4,] -0.1 2.00e-01 -1.04e-16 1.00e+00> cs1AR1 <- corAR1(0.8, form = ~ 1 | Subject)> cs1AR1 <- Initialize(cs1AR1, data = Orthodont)> corMatrix(cs1AR1)$M01[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M02[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M03[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M04[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M05[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M06[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M07[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M08[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M09[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M10[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M11[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M12[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M13[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M14[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M15[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$M16[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$F01[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$F02[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$F03[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$F04[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$F05[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$F06[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$F07[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$F08[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$F09[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$F10[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000$F11[,1] [,2] [,3] [,4][1,] 1.000 0.80 0.64 0.512[2,] 0.800 1.00 0.80 0.640[3,] 0.640 0.80 1.00 0.800[4,] 0.512 0.64 0.80 1.000> cs1ARMA <- corARMA(0.4, form = ~ 1 | Subject, q = 1)> cs1ARMA <- Initialize(cs1ARMA, data = Orthodont)> corMatrix(cs1ARMA)$M01[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M02[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M03[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M04[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M05[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M06[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M07[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M08[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M09[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M10[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M11[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M12[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M13[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M14[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M15[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$M16[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$F01[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$F02[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$F03[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$F04[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$F05[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$F06[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$F07[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$F08[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$F09[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$F10[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000$F11[,1] [,2] [,3] [,4][1,] 1.000 0.345 0.000 0.000[2,] 0.345 1.000 0.345 0.000[3,] 0.000 0.345 1.000 0.345[4,] 0.000 0.000 0.345 1.000> cs2ARMA <- corARMA(c(0.8, 0.4), form = ~ 1 | Subject, p=1, q=1)> cs2ARMA <- Initialize(cs2ARMA, data = Orthodont)> corMatrix(cs2ARMA)$M01[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M02[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M03[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M04[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M05[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M06[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M07[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M08[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M09[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M10[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M11[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M12[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M13[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M14[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M15[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$M16[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$F01[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$F02[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$F03[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$F04[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$F05[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$F06[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$F07[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$F08[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$F09[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$F10[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000$F11[,1] [,2] [,3] [,4][1,] 1.000 0.880 0.704 0.563[2,] 0.880 1.000 0.880 0.704[3,] 0.704 0.880 1.000 0.880[4,] 0.563 0.704 0.880 1.000> spatDat <- data.frame(x = (0:4)/4, y = (0:4)/4)> cs1Exp <- corExp(1, form = ~ x + y)> cs1Exp <- Initialize(cs1Exp, spatDat)> corMatrix(cs1Exp)[,1] [,2] [,3] [,4] [,5][1,] 1.000 0.702 0.493 0.346 0.243[2,] 0.702 1.000 0.702 0.493 0.346[3,] 0.493 0.702 1.000 0.702 0.493[4,] 0.346 0.493 0.702 1.000 0.702[5,] 0.243 0.346 0.493 0.702 1.000> cs2Exp <- corExp(1, form = ~ x + y, metric = "man")> cs2Exp <- Initialize(cs2Exp, spatDat)> corMatrix(cs2Exp)[,1] [,2] [,3] [,4] [,5][1,] 1.000 0.607 0.368 0.223 0.135[2,] 0.607 1.000 0.607 0.368 0.223[3,] 0.368 0.607 1.000 0.607 0.368[4,] 0.223 0.368 0.607 1.000 0.607[5,] 0.135 0.223 0.368 0.607 1.000> cs3Exp <- corExp(c(1, 0.2), form = ~ x + y, nugget = TRUE)> cs3Exp <- Initialize(cs3Exp, spatDat)> corMatrix(cs3Exp)[,1] [,2] [,3] [,4] [,5][1,] 1.000 0.562 0.394 0.277 0.194[2,] 0.562 1.000 0.562 0.394 0.277[3,] 0.394 0.562 1.000 0.562 0.394[4,] 0.277 0.394 0.562 1.000 0.562[5,] 0.194 0.277 0.394 0.562 1.000> fm1Ovar.lme <- lme(follicles ~ sin(2*pi*Time) + cos(2*pi*Time),+ data = Ovary, random = pdDiag(~sin(2*pi*Time)))> fm1Ovar.lmeLinear mixed-effects model fit by REMLData: OvaryLog-restricted-likelihood: -813Fixed: follicles ~ sin(2 * pi * Time) + cos(2 * pi * Time)(Intercept) sin(2 * pi * Time) cos(2 * pi * Time)12.182 -3.299 -0.862Random effects:Formula: ~sin(2 * pi * Time) | MareStructure: Diagonal(Intercept) sin(2 * pi * Time) ResidualStdDev: 3.05 2.08 3.11Number of Observations: 308Number of Groups: 11> ACF(fm1Ovar.lme)lag ACF1 0 1.00002 1 0.37953 2 0.17974 3 0.03575 4 0.05986 5 0.00217 6 0.06438 7 0.07169 8 0.048610 9 0.027811 10 -0.034312 11 -0.077213 12 -0.161114 13 -0.196015 14 -0.2893> plot(ACF(fm1Ovar.lme, maxLag = 10), alpha = 0.01)> fm2Ovar.lme <- update(fm1Ovar.lme, correlation = corAR1())> anova(fm1Ovar.lme, fm2Ovar.lme)Model df AIC BIC logLik Test L.Ratio p-valuefm1Ovar.lme 1 6 1638 1660 -813fm2Ovar.lme 2 7 1563 1589 -775 1 vs 2 76.6 <.0001> if (interactive()) intervals(fm2Ovar.lme)> fm3Ovar.lme <- update(fm1Ovar.lme, correlation = corARMA(q = 2))> fm3Ovar.lmeLinear mixed-effects model fit by REMLData: OvaryLog-restricted-likelihood: -778Fixed: follicles ~ sin(2 * pi * Time) + cos(2 * pi * Time)(Intercept) sin(2 * pi * Time) cos(2 * pi * Time)12.194 -3.115 -0.869Random effects:Formula: ~sin(2 * pi * Time) | MareStructure: Diagonal(Intercept) sin(2 * pi * Time) ResidualStdDev: 2.97 1.67 3.24Correlation Structure: ARMA(0,2)Formula: ~1 | MareParameter estimate(s):Theta1 Theta20.475 0.257Number of Observations: 308Number of Groups: 11> anova(fm2Ovar.lme, fm3Ovar.lme, test = F)Model df AIC BIC logLikfm2Ovar.lme 1 7 1563 1589 -775fm3Ovar.lme 2 8 1571 1601 -778> fm4Ovar.lme <- update(fm1Ovar.lme,+ correlation = corCAR1(form = ~Time))> anova(fm2Ovar.lme, fm4Ovar.lme, test = F)Model df AIC BIC logLikfm2Ovar.lme 1 7 1563 1589 -775fm4Ovar.lme 2 7 1566 1592 -776> (fm5Ovar.lme <- update(fm1Ovar.lme,+ corr = corARMA(p = 1, q = 1)))Linear mixed-effects model fit by REMLData: OvaryLog-restricted-likelihood: -772Fixed: follicles ~ sin(2 * pi * Time) + cos(2 * pi * Time)(Intercept) sin(2 * pi * Time) cos(2 * pi * Time)12.125 -2.920 -0.849Random effects:Formula: ~sin(2 * pi * Time) | MareStructure: Diagonal(Intercept) sin(2 * pi * Time) ResidualStdDev: 2.61 1 3.73Correlation Structure: ARMA(1,1)Formula: ~1 | MareParameter estimate(s):Phi1 Theta10.787 -0.279Number of Observations: 308Number of Groups: 11> anova(fm2Ovar.lme, fm5Ovar.lme)Model df AIC BIC logLik Test L.Ratio p-valuefm2Ovar.lme 1 7 1563 1589 -775fm5Ovar.lme 2 8 1560 1590 -772 1 vs 2 5.55 0.0184> plot(ACF(fm5Ovar.lme, maxLag = 10, resType = "n"), alpha = 0.01)> Variogram(fm2BW.lme, form = ~ Time)variog dist n.pairs1 0.345 1 162 0.993 6 163 0.762 7 1444 0.685 8 165 0.682 13 166 0.951 14 1287 0.900 15 168 1.694 20 169 1.125 21 11210 1.088 22 1611 0.897 28 9612 0.932 29 1613 0.851 35 8014 0.755 36 1615 1.082 42 6416 1.567 43 1617 0.644 49 4818 0.674 56 3219 0.587 63 16> plot(Variogram(fm2BW.lme, form = ~ Time, maxDist = 42))> fm3BW.lme <- update(fm2BW.lme,+ correlation = corExp(form = ~ Time))> ## IGNORE_RDIFF_BEGIN> intervals(fm3BW.lme)Approximate 95% confidence intervalsFixed effects:lower est. upper(Intercept) 2.26e+02 251.487 277.336Time 1.93e-01 0.363 0.532Diet2 1.52e+02 200.786 249.841Diet3 2.04e+02 252.590 301.667Time:Diet2 3.22e-01 0.624 0.926Time:Diet3 2.63e-03 0.307 0.610Random Effects:Level: Ratlower est. uppersd((Intercept)) 25.023 36.919 54.471sd(Time) 0.147 0.233 0.368cor((Intercept),Time) -0.637 -0.147 0.428Correlation structure:lower est. upperrange 2.46 4.89 9.7Variance function:lower est. upperpower 0.244 0.594 0.944Within-group standard error:lower est. upper0.0181 0.1384 1.0593> ## IGNORE_RDIFF_END> anova(fm2BW.lme, fm3BW.lme)Model df AIC BIC logLik Test L.Ratio p-valuefm2BW.lme 1 11 1164 1198 -571fm3BW.lme 2 12 1145 1183 -561 1 vs 2 20.8 <.0001> fm4BW.lme <-+ update(fm3BW.lme, correlation = corExp(form = ~ Time,+ nugget = TRUE))> anova(fm3BW.lme, fm4BW.lme)Model df AIC BIC logLik Test L.Ratio p-valuefm3BW.lme 1 12 1145 1183 -561fm4BW.lme 2 13 1138 1178 -556 1 vs 2 9.5 0.0021> plot(Variogram(fm3BW.lme, form = ~ Time, maxDist = 42))> plot(Variogram(fm3BW.lme, form = ~ Time, maxDist = 42,+ resType = "n", robust = TRUE))> fm5BW.lme <- update(fm3BW.lme, correlation = corRatio(form = ~ Time))> fm6BW.lme <- update(fm3BW.lme, correlation = corSpher(form = ~ Time))> fm7BW.lme <- update(fm3BW.lme, correlation = corLin(form = ~ Time))> fm8BW.lme <- update(fm3BW.lme, correlation = corGaus(form = ~ Time))> anova(fm3BW.lme, fm5BW.lme, fm6BW.lme, fm7BW.lme, fm8BW.lme)Model df AIC BIC logLikfm3BW.lme 1 12 1145 1183 -561fm5BW.lme 2 12 1149 1186 -562fm6BW.lme 3 12 1151 1188 -563fm7BW.lme 4 12 1151 1188 -563fm8BW.lme 5 12 1151 1188 -563> fm1Orth.gls <- gls(distance ~ Sex * I(age - 11), Orthodont,+ correlation = corSymm(form = ~ 1 | Subject),+ weights = varIdent(form = ~ 1 | age))> fm1Orth.glsGeneralized least squares fit by REMLModel: distance ~ Sex * I(age - 11)Data: OrthodontLog-restricted-likelihood: -212Coefficients:(Intercept) SexFemale24.937 -2.272I(age - 11) SexFemale:I(age - 11)0.827 -0.350Correlation Structure: GeneralFormula: ~1 | SubjectParameter estimate(s):Correlation:1 2 32 0.5683 0.659 0.5814 0.522 0.725 0.740Variance function:Structure: Different standard deviations per stratumFormula: ~1 | ageParameter estimates:8 10 12 141.000 0.879 1.074 0.959Degrees of freedom: 108 total; 104 residualResidual standard error: 2.33> ## IGNORE_RDIFF_BEGIN> intervals(fm1Orth.gls)Approximate 95% confidence intervalsCoefficients:lower est. upper(Intercept) 23.999 24.937 25.875SexFemale -3.741 -2.272 -0.803I(age - 11) 0.664 0.827 0.990SexFemale:I(age - 11) -0.606 -0.350 -0.095Correlation structure:lower est. uppercor(1,2) 0.253 0.568 0.774cor(1,3) 0.385 0.659 0.826cor(1,4) 0.184 0.522 0.749cor(2,3) 0.272 0.581 0.781cor(2,4) 0.481 0.725 0.865cor(3,4) 0.512 0.740 0.870Variance function:lower est. upper10 0.633 0.879 1.2212 0.801 1.074 1.4414 0.686 0.959 1.34Residual standard error:lower est. upper1.77 2.33 3.07> ## IGNORE_RDIFF_END> fm2Orth.gls <-+ update(fm1Orth.gls, corr = corCompSymm(form = ~ 1 | Subject))> anova(fm1Orth.gls, fm2Orth.gls)Model df AIC BIC logLik Test L.Ratio p-valuefm1Orth.gls 1 14 453 490 -212fm2Orth.gls 2 9 450 474 -216 1 vs 2 7.43 0.191> intervals(fm2Orth.gls)Approximate 95% confidence intervalsCoefficients:lower est. upper(Intercept) 23.930 24.868 25.8071SexFemale -3.668 -2.197 -0.7266I(age - 11) 0.642 0.794 0.9470SexFemale:I(age - 11) -0.555 -0.316 -0.0763Correlation structure:lower est. upperRho 0.446 0.635 0.778Variance function:lower est. upper10 0.638 0.862 1.1712 0.771 1.034 1.3914 0.683 0.920 1.24Residual standard error:lower est. upper1.81 2.39 3.15> fm3Orth.gls <- update(fm2Orth.gls, weights = NULL)> anova(fm2Orth.gls, fm3Orth.gls)Model df AIC BIC logLik Test L.Ratio p-valuefm2Orth.gls 1 9 450 474 -216fm3Orth.gls 2 6 446 462 -217 1 vs 2 1.78 0.618> plot(fm3Orth.gls, resid(., type = "n") ~ age | Sex)> fm4Orth.gls <- update(fm3Orth.gls,+ weights = varIdent(form = ~ 1 | Sex))> anova(fm3Orth.gls, fm4Orth.gls)Model df AIC BIC logLik Test L.Ratio p-valuefm3Orth.gls 1 6 446 462 -217fm4Orth.gls 2 7 436 455 -211 1 vs 2 11.6 7e-04> qqnorm(fm4Orth.gls, ~resid(., type = "n"))> # not in book but needed for the following command> fm3Orth.lme <-+ lme(distance~Sex*I(age-11), data = Orthodont,+ random = ~ I(age-11) | Subject,+ weights = varIdent(form = ~ 1 | Sex))> anova(fm3Orth.lme, fm4Orth.gls, test = FALSE)Model df AIC BIC logLikfm3Orth.lme 1 9 430 453 -206fm4Orth.gls 2 7 436 455 -211> fm1Dial.gls <-+ gls(rate ~(pressure + I(pressure^2) + I(pressure^3) + I(pressure^4))*QB,+ Dialyzer)> plot(fm1Dial.gls, resid(.) ~ pressure, abline = 0)> fm2Dial.gls <- update(fm1Dial.gls,+ weights = varPower(form = ~ pressure))> anova(fm1Dial.gls, fm2Dial.gls)Model df AIC BIC logLik Test L.Ratio p-valuefm1Dial.gls 1 11 761 793 -370fm2Dial.gls 2 12 738 773 -357 1 vs 2 24.9 <.0001> ACF(fm2Dial.gls, form = ~ 1 | Subject)lag ACF1 0 1.00002 1 0.77093 2 0.63234 3 0.40835 4 0.20076 5 0.07317 6 0.0778> plot(ACF(fm2Dial.gls, form = ~ 1 | Subject), alpha = 0.01)> (fm3Dial.gls <- update(fm2Dial.gls,+ corr = corAR1(0.771, form = ~ 1 | Subject)))Generalized least squares fit by REMLModel: rate ~ (pressure + I(pressure^2) + I(pressure^3) + I(pressure^4)) * QBData: DialyzerLog-restricted-likelihood: -308Coefficients:(Intercept) pressure I(pressure^2)-16.818 92.334 -49.265I(pressure^3) I(pressure^4) QB30011.400 -1.020 -1.594pressure:QB300 I(pressure^2):QB300 I(pressure^3):QB3001.705 2.127 0.480I(pressure^4):QB300-0.221Correlation Structure: AR(1)Formula: ~1 | SubjectParameter estimate(s):Phi0.753Variance function:Structure: Power of variance covariateFormula: ~pressureParameter estimates:power0.518Degrees of freedom: 140 total; 130 residualResidual standard error: 3.05> intervals(fm3Dial.gls)Approximate 95% confidence intervalsCoefficients:lower est. upper(Intercept) -18.90 -16.818 -14.7401pressure 81.91 92.334 102.7541I(pressure^2) -63.10 -49.265 -35.4263I(pressure^3) 4.56 11.400 18.2345I(pressure^4) -2.12 -1.020 0.0856QB300 -4.76 -1.594 1.5681pressure:QB300 -13.64 1.705 17.0518I(pressure^2):QB300 -17.95 2.127 22.2020I(pressure^3):QB300 -9.35 0.480 10.3097I(pressure^4):QB300 -1.80 -0.221 1.3608Correlation structure:lower est. upperPhi 0.628 0.753 0.839Variance function:lower est. upperpower 0.381 0.518 0.656Residual standard error:lower est. upper2.50 3.05 3.71> anova(fm2Dial.gls, fm3Dial.gls)Model df AIC BIC logLik Test L.Ratio p-valuefm2Dial.gls 1 12 738 773 -357fm3Dial.gls 2 13 643 680 -308 1 vs 2 97.5 <.0001> anova(fm3Dial.gls, fm2Dial.lme, test = FALSE)Model df AIC BIC logLikfm3Dial.gls 1 13 643 680 -308fm2Dial.lme 2 18 655 707 -310> fm1Wheat2 <- gls(yield ~ variety - 1, Wheat2)> Variogram(fm1Wheat2, form = ~ latitude + longitude)variog dist n.pairs1 0.370 4.30 11432 0.396 5.61 12593 0.470 8.39 12634 0.508 9.32 12415 0.545 10.52 12426 0.640 12.75 12417 0.612 13.39 12838 0.657 14.76 12529 0.738 16.18 122110 0.728 17.37 126111 0.751 18.46 128812 0.875 20.24 125413 0.805 21.63 125614 0.871 22.67 118215 0.868 24.62 125716 0.859 26.24 126417 0.971 28.56 123518 0.993 30.79 122619 1.096 34.59 126320 1.341 39.36 1234> plot(Variogram(fm1Wheat2, form = ~ latitude + longitude,+ maxDist = 32), xlim = c(0,32))> fm2Wheat2 <- update(fm1Wheat2, corr = corSpher(c(28, 0.2),+ form = ~ latitude + longitude,+ nugget = TRUE))> fm2Wheat2Generalized least squares fit by REMLModel: yield ~ variety - 1Data: Wheat2Log-restricted-likelihood: -534Coefficients:varietyARAPAHOE varietyBRULE varietyBUCKSKIN26.7 25.8 34.8varietyCENTURA varietyCENTURK78 varietyCHEYENNE25.1 26.3 24.7varietyCODY varietyCOLT varietyGAGE22.5 25.2 24.3varietyHOMESTEAD varietyKS831374 varietyLANCER21.7 26.9 23.3varietyLANCOTA varietyNE83404 varietyNE8340621.3 24.0 25.3varietyNE83407 varietyNE83432 varietyNE8349825.2 21.8 28.7varietyNE83T12 varietyNE84557 varietyNE8555622.1 21.8 28.0varietyNE85623 varietyNE86482 varietyNE8650123.9 25.0 25.0varietyNE86503 varietyNE86507 varietyNE8650927.2 27.5 22.4varietyNE86527 varietyNE86582 varietyNE8660625.9 22.6 26.8varietyNE86607 varietyNE86T666 varietyNE8740325.9 16.8 21.5varietyNE87408 varietyNE87409 varietyNE8744624.3 26.3 22.2varietyNE87451 varietyNE87457 varietyNE8746324.2 23.5 23.2varietyNE87499 varietyNE87512 varietyNE8751322.2 22.6 21.8varietyNE87522 varietyNE87612 varietyNE8761319.5 27.4 27.6varietyNE87615 varietyNE87619 varietyNE8762723.8 28.5 18.5varietyNORKAN varietyREDLAND varietyROUGHRIDER22.1 28.0 25.7varietySCOUT66 varietySIOUXLAND varietyTAM10726.9 25.7 22.8varietyTAM200 varietyVONA18.8 24.8Correlation Structure: Spherical spatial correlationFormula: ~latitude + longitudeParameter estimate(s):range nugget27.457 0.209Degrees of freedom: 224 total; 168 residualResidual standard error: 7.41> fm3Wheat2 <- update(fm1Wheat2,+ corr = corRatio(c(12.5, 0.2),+ form = ~ latitude + longitude, nugget = TRUE))> fm3Wheat2Generalized least squares fit by REMLModel: yield ~ variety - 1Data: Wheat2Log-restricted-likelihood: -533Coefficients:varietyARAPAHOE varietyBRULE varietyBUCKSKIN26.5 26.3 35.0varietyCENTURA varietyCENTURK78 varietyCHEYENNE24.9 26.7 24.4varietyCODY varietyCOLT varietyGAGE23.4 25.2 24.5varietyHOMESTEAD varietyKS831374 varietyLANCER21.5 26.5 23.0varietyLANCOTA varietyNE83404 varietyNE8340621.2 24.6 25.7varietyNE83407 varietyNE83432 varietyNE8349825.5 21.8 29.1varietyNE83T12 varietyNE84557 varietyNE8555621.6 21.3 27.9varietyNE85623 varietyNE86482 varietyNE8650123.7 24.4 24.9varietyNE86503 varietyNE86507 varietyNE8650927.3 27.4 22.2varietyNE86527 varietyNE86582 varietyNE8660625.0 23.3 27.3varietyNE86607 varietyNE86T666 varietyNE8740325.7 17.3 21.8varietyNE87408 varietyNE87409 varietyNE8744624.7 26.3 22.1varietyNE87451 varietyNE87457 varietyNE8746324.4 23.6 23.4varietyNE87499 varietyNE87512 varietyNE8751321.9 22.7 21.6varietyNE87522 varietyNE87612 varietyNE8761319.6 28.3 27.7varietyNE87615 varietyNE87619 varietyNE8762724.0 28.7 19.1varietyNORKAN varietyREDLAND varietyROUGHRIDER22.7 27.7 25.6varietySCOUT66 varietySIOUXLAND varietyTAM10726.3 25.7 22.5varietyTAM200 varietyVONA18.7 25.0Correlation Structure: Rational quadratic spatial correlationFormula: ~latitude + longitudeParameter estimate(s):range nugget13.461 0.194Degrees of freedom: 224 total; 168 residualResidual standard error: 8.85> anova(fm2Wheat2, fm3Wheat2)Model df AIC BIC logLikfm2Wheat2 1 59 1186 1370 -534fm3Wheat2 2 59 1183 1368 -533> anova(fm1Wheat2, fm3Wheat2)Model df AIC BIC logLik Test L.Ratio p-valuefm1Wheat2 1 57 1355 1533 -620fm3Wheat2 2 59 1183 1368 -533 1 vs 2 176 <.0001> plot(Variogram(fm3Wheat2, resType = "n"))> plot(fm3Wheat2, resid(., type = "n") ~ fitted(.), abline = 0)> qqnorm(fm3Wheat2, ~ resid(., type = "n"))> fm4Wheat2 <- update(fm3Wheat2, model = yield ~ variety)> anova(fm4Wheat2)Denom. DF: 168numDF F-value p-value(Intercept) 1 30.40 <.0001variety 55 1.85 0.0015> anova(fm3Wheat2, L = c(-1, 0, 1))Denom. DF: 168F-test for linear combination(s)varietyARAPAHOE varietyBUCKSKIN-1 1numDF F-value p-value1 1 7.7 0.0062> # cleanup>> summary(warnings())No warnings======ch06.R======> #-*- R -*->> # initialization>> library(nlme)> options(width = 65,+ ## reduce platform dependence in printed output when testing+ digits = if(nzchar(Sys.getenv("R_TESTS"))) 3 else 5)> options(contrasts = c(unordered = "contr.helmert", ordered = "contr.poly"))> pdf(file = "ch06.pdf")> # Chapter 6 Nonlinear Mixed-Effects Models:> # Basic Concepts and Motivating Examples>> # 6.2 Indomethicin Kinetics>> plot(Indometh)> fm1Indom.nls <- nls(conc ~ SSbiexp(time, A1, lrc1, A2, lrc2),+ data = Indometh)> summary(fm1Indom.nls)Formula: conc ~ SSbiexp(time, A1, lrc1, A2, lrc2)Parameters:Estimate Std. Error t value Pr(>|t|)A1 2.773 0.253 10.95 4e-16 ***lrc1 0.886 0.222 3.99 0.00018 ***A2 0.607 0.267 2.27 0.02660 *lrc2 -1.092 0.409 -2.67 0.00966 **---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 0.174 on 62 degrees of freedomNumber of iterations to convergence: 0Achieved convergence tolerance: 3.3e-07> plot(fm1Indom.nls, Subject ~ resid(.), abline = 0)> (fm1Indom.lis <- nlsList(conc ~ SSbiexp(time, A1, lrc1, A2, lrc2),+ data = Indometh))Call:Model: conc ~ SSbiexp(time, A1, lrc1, A2, lrc2) | SubjectData: IndomethCoefficients:A1 lrc1 A2 lrc21 2.03 0.579 0.192 -1.7884 2.20 0.242 0.255 -1.6032 2.83 0.801 0.499 -1.6355 3.57 1.041 0.291 -1.5076 3.00 1.088 0.969 -0.8733 5.47 1.750 1.676 -0.412Degrees of freedom: 66 total; 42 residualResidual standard error: 0.0756> plot(intervals(fm1Indom.lis))> ## IGNORE_RDIFF_BEGIN> (fm1Indom.nlme <- nlme(fm1Indom.lis,+ random = pdDiag(A1 + lrc1 + A2 + lrc2 ~ 1),+ control = list(tolerance = 0.0001)))Nonlinear mixed-effects model fit by maximum likelihoodModel: conc ~ SSbiexp(time, A1, lrc1, A2, lrc2)Data: IndomethLog-likelihood: 54.6Fixed: list(A1 ~ 1, lrc1 ~ 1, A2 ~ 1, lrc2 ~ 1)A1 lrc1 A2 lrc22.828 0.774 0.461 -1.344Random effects:Formula: list(A1 ~ 1, lrc1 ~ 1, A2 ~ 1, lrc2 ~ 1)Level: SubjectStructure: DiagonalA1 lrc1 A2 lrc2 ResidualStdDev: 0.571 0.158 0.112 7.32e-06 0.0815Number of Observations: 66Number of Groups: 6> ## IGNORE_RDIFF_END> fm2Indom.nlme <- update(fm1Indom.nlme,+ random = pdDiag(A1 + lrc1 + A2 ~ 1))> anova(fm1Indom.nlme, fm2Indom.nlme)Model df AIC BIC logLik Test L.Ratio p-valuefm1Indom.nlme 1 9 -91.2 -71.5 54.6fm2Indom.nlme 2 8 -93.2 -75.7 54.6 1 vs 2 0.00871 0.926> (fm3Indom.nlme <- update(fm2Indom.nlme, random = A1+lrc1+A2 ~ 1))Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!Nonlinear mixed-effects model fit by maximum likelihoodModel: conc ~ SSbiexp(time, A1, lrc1, A2, lrc2)Data: IndomethLog-likelihood: 58.5Fixed: list(A1 ~ 1, lrc1 ~ 1, A2 ~ 1, lrc2 ~ 1)A1 lrc1 A2 lrc22.815 0.829 0.561 -1.141Random effects:Formula: list(A1 ~ 1, lrc1 ~ 1, A2 ~ 1)Level: SubjectStructure: General positive-definite, Log-Cholesky parametrizationStdDev CorrA1 0.6904 A1 lrc1lrc1 0.1790 0.932A2 0.1537 0.471 0.118Residual 0.0781Number of Observations: 66Number of Groups: 6> fm4Indom.nlme <-+ update(fm3Indom.nlme,+ random = pdBlocked(list(A1 + lrc1 ~ 1, A2 ~ 1)))> ## IGNORE_RDIFF_BEGIN> anova(fm3Indom.nlme, fm4Indom.nlme)Model df AIC BIC logLik Test L.Ratio p-valuefm3Indom.nlme 1 11 -94.9 -70.9 58.5fm4Indom.nlme 2 9 -98.2 -78.4 58.1 1 vs 2 0.789 0.674> ## IGNORE_RDIFF_END> anova(fm2Indom.nlme, fm4Indom.nlme)Model df AIC BIC logLik Test L.Ratio p-valuefm2Indom.nlme 1 8 -93.2 -75.7 54.6fm4Indom.nlme 2 9 -98.2 -78.4 58.1 1 vs 2 6.97 0.0083> plot(fm4Indom.nlme, id = 0.05, adj = -1)> qqnorm(fm4Indom.nlme)> plot(augPred(fm4Indom.nlme, level = 0:1))> summary(fm4Indom.nlme)Nonlinear mixed-effects model fit by maximum likelihoodModel: conc ~ SSbiexp(time, A1, lrc1, A2, lrc2)Data: IndomethAIC BIC logLik-98.2 -78.4 58.1Random effects:Composite Structure: BlockedBlock 1: A1, lrc1Formula: list(A1 ~ 1, lrc1 ~ 1)Level: SubjectStructure: General positive-definiteStdDev CorrA1 0.720 A1lrc1 0.149 1Block 2: A2Formula: A2 ~ 1 | SubjectA2 ResidualStdDev: 0.213 0.0782Fixed effects: list(A1 ~ 1, lrc1 ~ 1, A2 ~ 1, lrc2 ~ 1)Value Std.Error DF t-value p-valueA1 2.783 0.327 57 8.51 0lrc1 0.898 0.111 57 8.11 0A2 0.658 0.143 57 4.61 0lrc2 -1.000 0.150 57 -6.67 0Correlation:A1 lrc1 A2lrc1 0.602A2 -0.058 0.556lrc2 -0.109 0.570 0.702Standardized Within-Group Residuals:Min Q1 Med Q3 Max-3.459 -0.437 0.110 0.504 3.057Number of Observations: 66Number of Groups: 6> # 6.3 Growth of Soybean Plants>> head(Soybean)Grouped Data: weight ~ Time | PlotPlot Variety Year Time weight1 1988F1 F 1988 14 0.1062 1988F1 F 1988 21 0.2613 1988F1 F 1988 28 0.6664 1988F1 F 1988 35 2.1105 1988F1 F 1988 42 3.5606 1988F1 F 1988 49 6.230> plot(Soybean, outer = ~ Year * Variety)> (fm1Soy.lis <- nlsList(weight ~ SSlogis(Time, Asym, xmid, scal),+ data = Soybean,+ ## in R >= 3.4.3, more iterations are needed for "1989P5"+ ## due to a change of initial values in SSlogis();+ ## control is passed to getInitial() only since R 4.1.0+ control = list(maxiter = 60)))Warning: 1 error caught in nls(y ~ 1/(1 + exp((xmid - x)/scal)), data = xy, start = list(xmid = aux[[1L]], scal = aux[[2L]]), algorithm = "plinear", ...): step factor 0.000488281 reduced below 'minFactor' of 0.000976562Call:Model: weight ~ SSlogis(Time, Asym, xmid, scal) | PlotData: SoybeanCoefficients:Asym xmid scal1988F4 15.15 52.8 5.181988F2 19.75 56.6 8.411988F1 20.34 57.4 9.601988F7 19.87 56.2 8.071988F5 30.65 64.1 11.261988F8 22.78 59.3 9.001988F6 23.29 59.6 9.721988F3 23.70 56.4 7.641988P1 17.30 48.8 6.361988P5 17.70 51.3 6.811988P4 24.01 57.8 11.741988P8 28.25 63.0 10.951988P7 27.49 61.5 10.181988P3 24.94 56.3 8.321988P2 36.66 66.6 11.921988P6 163.70 105.0 17.931989F6 8.51 55.3 8.861989F5 9.67 51.3 7.211989F4 11.25 53.8 6.491989F1 11.25 56.6 6.071989F2 11.23 52.2 7.021989F7 10.07 51.4 5.501989F8 10.61 48.0 5.961989F3 18.42 66.1 9.221989P7 15.47 46.3 5.391989P4 18.18 57.2 8.401989P6 20.50 58.2 10.611989P5 17.89 54.1 6.051989P1 21.68 59.7 9.971989P3 22.28 53.4 7.901989P2 28.30 67.2 12.521989P8 NA NA NA1990F2 19.46 66.3 13.161990F3 19.87 58.3 12.801990F4 27.44 70.3 14.561990F5 18.72 51.3 7.761990F1 19.79 55.7 9.621990F8 20.29 55.5 7.771990F7 19.84 54.7 6.791990F6 21.20 54.6 9.261990P8 18.51 52.4 8.581990P7 19.16 54.8 10.851990P3 19.20 49.7 9.321990P1 18.45 47.9 6.611990P6 17.69 50.2 6.631990P5 19.54 51.2 7.291990P2 25.79 62.4 11.661990P4 26.13 61.2 10.97Degrees of freedom: 404 total; 263 residualResidual standard error: 1.04> ## IGNORE_RDIFF_BEGIN> (fm1Soy.nlme <- nlme(fm1Soy.lis))Nonlinear mixed-effects model fit by maximum likelihoodModel: weight ~ SSlogis(Time, Asym, xmid, scal)Data: SoybeanLog-likelihood: -740Fixed: list(Asym ~ 1, xmid ~ 1, scal ~ 1)Asym xmid scal19.3 55.0 8.4Random effects:Formula: list(Asym ~ 1, xmid ~ 1, scal ~ 1)Level: PlotStructure: General positive-definite, Log-Cholesky parametrizationStdDev CorrAsym 5.20 Asym xmidxmid 4.20 0.721scal 1.40 0.711 0.959Residual 1.12Number of Observations: 412Number of Groups: 48> ## IGNORE_RDIFF_END> fm2Soy.nlme <- update(fm1Soy.nlme, weights = varPower())Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :Iteration 6, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)> anova(fm1Soy.nlme, fm2Soy.nlme)Model df AIC BIC logLik Test L.Ratio p-valuefm1Soy.nlme 1 10 1500 1540 -740fm2Soy.nlme 2 11 746 790 -362 1 vs 2 756 <.0001> plot(ranef(fm2Soy.nlme, augFrame = TRUE),+ form = ~ Year * Variety, layout = c(3,1))> soyFix <- fixef(fm2Soy.nlme)> options(contrasts = c("contr.treatment", "contr.poly"))> ## IGNORE_RDIFF_BEGIN> (fm3Soy.nlme <-+ update(fm2Soy.nlme,+ fixed = Asym + xmid + scal ~ Year,+ start = c(soyFix[1], 0, 0, soyFix[2], 0, 0, soyFix[3], 0, 0)))Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :Iteration 6, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)Nonlinear mixed-effects model fit by maximum likelihoodModel: weight ~ SSlogis(Time, Asym, xmid, scal)Data: SoybeanLog-likelihood: -326Fixed: Asym + xmid + scal ~ YearAsym.(Intercept) Asym.Year1989 Asym.Year199020.208 -6.303 -3.465xmid.(Intercept) xmid.Year1989 xmid.Year199054.099 -2.480 -4.848scal.(Intercept) scal.Year1989 scal.Year19908.051 -0.932 -0.662Random effects:Formula: list(Asym ~ 1, xmid ~ 1, scal ~ 1)Level: PlotStructure: General positive-definite, Log-Cholesky parametrizationStdDev CorrAsym.(Intercept) 2.71e+00 As.(I) xm.(I)xmid.(Intercept) 8.34e-12 0.992scal.(Intercept) 1.08e-01 0.999 0.993Residual 2.16e-01Variance function:Structure: Power of variance covariateFormula: ~fitted(.)Parameter estimates:power0.95Number of Observations: 412Number of Groups: 48> ## IGNORE_RDIFF_END> anova(fm3Soy.nlme)numDF denDF F-value p-valueAsym.(Intercept) 1 356 2057 <.0001Asym.Year 2 356 103 <.0001xmid.(Intercept) 1 356 11420 <.0001xmid.Year 2 356 9 1e-04scal.(Intercept) 1 356 7967 <.0001scal.Year 2 356 11 <.0001> # The following line is not in the book but needed to fit the model> fm4Soy.nlme <-+ nlme(weight ~ SSlogis(Time, Asym, xmid, scal),+ data = Soybean,+ fixed = list(Asym ~ Year*Variety, xmid ~ Year + Variety, scal ~ Year),+ random = Asym ~ 1,+ start = c(17, 0, 0, 0, 0, 0, 52, 0, 0, 0, 7.5, 0, 0),+ weights = varPower(0.95), control = list(verbose = TRUE))> # FIXME: An update doesn't work for the fixed argument when fixed is a list> ## p. 293-4 :> summary(fm4Soy.nlme)Nonlinear mixed-effects model fit by maximum likelihoodModel: weight ~ SSlogis(Time, Asym, xmid, scal)Data: SoybeanAIC BIC logLik616 681 -292Random effects:Formula: Asym ~ 1 | PlotAsym.(Intercept) ResidualStdDev: 1.04 0.218Variance function:Structure: Power of variance covariateFormula: ~fitted(.)Parameter estimates:power0.943Fixed effects: list(Asym ~ Year * Variety, xmid ~ Year + Variety, scal ~ Year)Value Std.Error DF t-value p-valueAsym.(Intercept) 19.4 0.954 352 20.4 0.0000Asym.Year1989 -8.8 1.072 352 -8.2 0.0000Asym.Year1990 -3.7 1.177 352 -3.1 0.0018Asym.VarietyP 1.6 1.038 352 1.6 0.1189Asym.Year1989:VarietyP 5.6 1.171 352 4.8 0.0000Asym.Year1990:VarietyP 0.1 1.176 352 0.1 0.9004xmid.(Intercept) 54.8 0.755 352 72.6 0.0000xmid.Year1989 -2.2 0.972 352 -2.3 0.0218xmid.Year1990 -5.0 0.974 352 -5.1 0.0000xmid.VarietyP -1.3 0.414 352 -3.1 0.0019scal.(Intercept) 8.1 0.147 352 54.8 0.0000scal.Year1989 -0.9 0.201 352 -4.4 0.0000scal.Year1990 -0.7 0.212 352 -3.2 0.0016Correlation:As.(I) As.Y1989 As.Y1990 Asy.VP A.Y1989:Asym.Year1989 -0.831Asym.Year1990 -0.736 0.646Asym.VarietyP -0.532 0.374 0.304Asym.Year1989:VarietyP 0.339 -0.403 -0.249 -0.662Asym.Year1990:VarietyP 0.318 -0.273 -0.447 -0.627 0.533xmid.(Intercept) 0.729 -0.595 -0.523 -0.144 0.007xmid.Year1989 -0.488 0.603 0.394 -0.021 0.133xmid.Year1990 -0.489 0.433 0.661 -0.016 0.020xmid.VarietyP -0.337 0.127 0.052 0.572 -0.114scal.(Intercept) 0.432 -0.381 -0.345 0.023 -0.029scal.Year1989 -0.311 0.369 0.252 -0.025 0.090scal.Year1990 -0.296 0.263 0.398 -0.023 0.022A.Y1990: xm.(I) x.Y198 x.Y199 xmd.VPAsym.Year1989Asym.Year1990Asym.VarietyPAsym.Year1989:VarietyPAsym.Year1990:VarietyPxmid.(Intercept) -0.011xmid.Year1989 0.021 -0.705xmid.Year1990 0.054 -0.706 0.545xmid.VarietyP -0.057 -0.308 0.006 0.015scal.(Intercept) -0.031 0.817 -0.629 -0.628 -0.022scal.Year1989 0.023 -0.593 0.855 0.459 0.002scal.Year1990 0.048 -0.563 0.437 0.840 0.004sc.(I) s.Y198Asym.Year1989Asym.Year1990Asym.VarietyPAsym.Year1989:VarietyPAsym.Year1990:VarietyPxmid.(Intercept)xmid.Year1989xmid.Year1990xmid.VarietyPscal.(Intercept)scal.Year1989 -0.731scal.Year1990 -0.694 0.507Standardized Within-Group Residuals:Min Q1 Med Q3 Max-2.628 -0.608 -0.124 0.570 3.919Number of Observations: 412Number of Groups: 48> plot(augPred(fm4Soy.nlme))# Fig 6.14, p. 295> # 6.4 Clinical Study of Phenobarbital Kinetics>> (fm1Pheno.nlme <-+ nlme(conc ~ phenoModel(Subject, time, dose, lCl, lV),+ data = Phenobarb, fixed = lCl + lV ~ 1,+ random = pdDiag(lCl + lV ~ 1), start = c(-5, 0),+ na.action = NULL, naPattern = ~ !is.na(conc)))Nonlinear mixed-effects model fit by maximum likelihoodModel: conc ~ phenoModel(Subject, time, dose, lCl, lV)Data: PhenobarbLog-likelihood: -505Fixed: lCl + lV ~ 1lCl lV-5.093 0.343Random effects:Formula: list(lCl ~ 1, lV ~ 1)Level: SubjectStructure: DiagonallCl lV ResidualStdDev: 0.44 0.45 2.79Number of Observations: 155Number of Groups: 59> fm1Pheno.ranef <- ranef(fm1Pheno.nlme, augFrame = TRUE)> # (These plots used to encounter difficulties, now fine):> plot(fm1Pheno.ranef, form = lCl ~ Wt + ApgarInd)> plot(fm1Pheno.ranef, form = lV ~ Wt + ApgarInd)> options(contrasts = c("contr.treatment", "contr.poly"))> if(FALSE)## This fit just "ping-pongs" until max.iterations error+ fm2Pheno.nlme <-+ update(fm1Pheno.nlme,+ fixed = list(lCl ~ Wt, lV ~ Wt + ApgarInd),+ start = c(-5.0935, 0, 0.34259, 0, 0),+ control = list(pnlsTol = 1e-4, maxIter = 500,+ msVerbose = TRUE, opt = "nlm"))> ##summary(fm2Pheno.nlme)> ##fm3Pheno.nlme <-> ## update(fm2Pheno.nlme,> ## fixed = lCl + lV ~ Wt,> ## start = fixef(fm2Pheno.nlme)[-5])> ##plot(fm3Pheno.nlme, conc ~ fitted(.), abline = c(0,1))>> # cleanup>> summary(warnings())No warnings======ch08.R======> #-*- R -*->> # initialization>> library(nlme)> library(lattice)> options(width = 65,+ ## reduce platform dependence in printed output when testing+ digits = if(nzchar(Sys.getenv("R_TESTS"))) 3 else 5)> options(contrasts = c(unordered = "contr.helmert", ordered = "contr.poly"))> pdf(file = "ch08.pdf")> # Chapter 8 Fitting Nonlinear Mixed-Effects Models>> # 8.1 Fitting Nonlinear Models in S with nls and nlsList>> ## outer = ~1 is used to display all five curves in one panel> plot(Orange, outer = ~1)> logist <-+ function(x, Asym, xmid, scal) Asym/(1 + exp(-(x - xmid)/scal))> logist <- deriv(~Asym/(1+exp(-(x-xmid)/scal)),+ c("Asym", "xmid", "scal"), function(x, Asym, xmid, scal){})> Asym <- 180; xmid <- 700; scal <- 300> logist(Orange$age[1:7], Asym, xmid, scal)[1] 22.6 58.9 84.6 132.1 153.8 162.7 171.0attr(,"gradient")Asym xmid scal[1,] 0.126 -0.0659 0.1279[2,] 0.327 -0.1321 0.0951[3,] 0.470 -0.1495 0.0179[4,] 0.734 -0.1172 -0.1188[5,] 0.854 -0.0746 -0.1321[6,] 0.904 -0.0522 -0.1169[7,] 0.950 -0.0286 -0.0841> fm1Oran.nls <- nls(circumference ~ logist(age, Asym, xmid, scal),+ data = Orange, start = c(Asym = 170, xmid = 700, scal = 500))> summary(fm1Oran.nls)Formula: circumference ~ logist(age, Asym, xmid, scal)Parameters:Estimate Std. Error t value Pr(>|t|)Asym 192.7 20.2 9.52 7.5e-11 ***xmid 728.8 107.3 6.79 1.1e-07 ***scal 353.5 81.5 4.34 0.00013 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 23.4 on 32 degrees of freedomNumber of iterations to convergence: 5Achieved convergence tolerance: 4.39e-06> plot(fm1Oran.nls)> plot(fm1Oran.nls, Tree ~ resid(.), abline = 0)> Orange.sortAvg <- sortedXyData("age", "circumference", Orange)> Orange.sortAvgx y1 118 31.02 484 57.83 664 93.24 1004 134.25 1231 145.66 1372 173.47 1582 175.8> NLSstClosestX(Orange.sortAvg, 130)[1] 969> logistInit <- function(mCall, LHS, data) {+ xy <- sortedXyData(mCall[["x"]], LHS, data)+ if(nrow(xy) < 3) {+ stop("Too few distinct input values to fit a logistic")+ }+ Asym <- max(abs(xy[,"y"]))+ if (Asym != max(xy[,"y"])) Asym <- -Asym # negative asymptote+ xmid <- NLSstClosestX(xy, 0.5 * Asym)+ scal <- NLSstClosestX(xy, 0.75 * Asym) - xmid+ value <- c(Asym, xmid, scal)+ names(value) <- mCall[c("Asym", "xmid", "scal")]+ value+ }> logist <- selfStart(logist, initial = logistInit)> class(logist)[1] "selfStart"> logist <- selfStart(~ Asym/(1 + exp(-(x - xmid)/scal)),+ initial = logistInit, parameters = c("Asym", "xmid", "scal"))> getInitial(circumference ~ logist(age, Asym, xmid, scal), Orange)Warning in getInitial.selfStart(func, data, mCall = as.list(match.call(func, :selfStart initializing functions should have a final '...' argument since R 4.1.0Asym xmid scal176 637 347> nls(circumference ~ logist(age, Asym, xmid, scal), Orange)Warning in getInitial.selfStart(func, data, mCall = as.list(match.call(func, :selfStart initializing functions should have a final '...' argument since R 4.1.0Nonlinear regression modelmodel: circumference ~ logist(age, Asym, xmid, scal)data: OrangeAsym xmid scal193 729 354residual sum-of-squares: 17480Number of iterations to convergence: 4Achieved convergence tolerance: 8.63e-07> getInitial(circumference ~ SSlogis(age,Asym,xmid,scal), Orange)Asym xmid scal193 729 354> nls(circumference ~ SSlogis(age, Asym, xmid, scal), Orange)Nonlinear regression modelmodel: circumference ~ SSlogis(age, Asym, xmid, scal)data: OrangeAsym xmid scal193 729 354residual sum-of-squares: 17480Number of iterations to convergence: 0Achieved convergence tolerance: 2.2e-06> fm1Oran.lis <-+ nlsList(circumference ~ SSlogis(age, Asym, xmid, scal) | Tree,+ data = Orange)> fm1Oran.lis <- nlsList(SSlogis, Orange)> fm1Oran.lis.noSS <-+ nlsList(circumference ~ Asym/(1+exp(-(age-xmid)/scal)),+ data = Orange,+ start = c(Asym = 170, xmid = 700, scal = 500))> fm1Oran.lisCall:Model: circumference ~ SSlogis(age, Asym, xmid, scal) | TreeData: OrangeCoefficients:Asym xmid scal3 159 735 4011 154 627 3635 207 861 3802 219 700 3324 225 711 303Degrees of freedom: 35 total; 20 residualResidual standard error: 7.98> summary(fm1Oran.lis)Call:Model: circumference ~ SSlogis(age, Asym, xmid, scal) | TreeData: OrangeCoefficients:AsymEstimate Std. Error t value Pr(>|t|)3 159 19.2 8.26 0.0004601 154 13.6 11.34 0.0001695 207 22.0 9.41 0.0007382 219 13.4 16.39 0.0001054 225 11.8 19.03 0.000104xmidEstimate Std. Error t value Pr(>|t|)3 735 130.8 5.62 0.0020111 627 92.9 6.75 0.0012635 861 108.0 7.98 0.0013892 700 61.4 11.42 0.0004354 711 51.2 13.89 0.000358scalEstimate Std. Error t value Pr(>|t|)3 401 94.8 4.23 0.005711 363 81.2 4.46 0.005865 380 66.8 5.69 0.004872 332 49.4 6.73 0.003244 303 41.6 7.29 0.00415Residual standard error: 7.98 on 20 degrees of freedom> plot(intervals(fm1Oran.lis), layout = c(3,1))> plot(fm1Oran.lis, Tree ~ resid(.), abline = 0)> Theoph[1:4,]Grouped Data: conc ~ Time | SubjectSubject Wt Dose Time conc1 1 79.6 4.02 0.00 0.742 1 79.6 4.02 0.25 2.843 1 79.6 4.02 0.57 6.574 1 79.6 4.02 1.12 10.50> fm1Theo.lis <- nlsList(conc ~ SSfol(Dose, Time, lKe, lKa, lCl),+ data = Theoph)> fm1Theo.lisCall:Model: conc ~ SSfol(Dose, Time, lKe, lKa, lCl) | SubjectData: TheophCoefficients:lKe lKa lCl6 -2.31 0.152 -2.977 -2.28 -0.386 -2.968 -2.39 0.319 -3.0711 -2.32 1.348 -2.863 -2.51 0.898 -3.232 -2.29 0.664 -3.114 -2.44 0.158 -3.299 -2.45 2.182 -3.4212 -2.25 -0.183 -3.1710 -2.60 -0.363 -3.431 -2.92 0.575 -3.925 -2.43 0.386 -3.13Degrees of freedom: 132 total; 96 residualResidual standard error: 0.7> plot(intervals(fm1Theo.lis), layout = c(3,1))> pairs(fm1Theo.lis, id = 0.1)> # 8.2 Fitting Nonlinear Mixed-Effects Models with nlme>> ## no need to specify groups, as Orange is a groupedData object> ## random is omitted - by default it is equal to fixed> (fm1Oran.nlme <-+ nlme(circumference ~ SSlogis(age, Asym, xmid, scal),+ data = Orange,+ fixed = Asym + xmid + scal ~ 1,+ start = fixef(fm1Oran.lis)))Warning in nlme.formula(circumference ~ SSlogis(age, Asym, xmid, scal), :Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!Nonlinear mixed-effects model fit by maximum likelihoodModel: circumference ~ SSlogis(age, Asym, xmid, scal)Data: OrangeLog-likelihood: -130Fixed: Asym + xmid + scal ~ 1Asym xmid scal192 728 357Random effects:Formula: list(Asym ~ 1, xmid ~ 1, scal ~ 1)Level: TreeStructure: General positive-definite, Log-Cholesky parametrizationStdDev CorrAsym 27.05 Asym xmidxmid 24.25 -0.328scal 36.60 -0.992 0.443Residual 7.32Number of Observations: 35Number of Groups: 5> summary(fm1Oran.nlme)Nonlinear mixed-effects model fit by maximum likelihoodModel: circumference ~ SSlogis(age, Asym, xmid, scal)Data: OrangeAIC BIC logLik280 296 -130Random effects:Formula: list(Asym ~ 1, xmid ~ 1, scal ~ 1)Level: TreeStructure: General positive-definite, Log-Cholesky parametrizationStdDev CorrAsym 27.05 Asym xmidxmid 24.25 -0.328scal 36.60 -0.992 0.443Residual 7.32Fixed effects: Asym + xmid + scal ~ 1Value Std.Error DF t-value p-valueAsym 192 14.1 28 13.7 0xmid 728 34.6 28 21.0 0scal 357 30.5 28 11.7 0Correlation:Asym xmidxmid 0.277scal -0.193 0.665Standardized Within-Group Residuals:Min Q1 Med Q3 Max-1.819 -0.522 0.174 0.518 1.645Number of Observations: 35Number of Groups: 5> summary(fm1Oran.nls)Formula: circumference ~ logist(age, Asym, xmid, scal)Parameters:Estimate Std. Error t value Pr(>|t|)Asym 192.7 20.2 9.52 7.5e-11 ***xmid 728.8 107.3 6.79 1.1e-07 ***scal 353.5 81.5 4.34 0.00013 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 23.4 on 32 degrees of freedomNumber of iterations to convergence: 5Achieved convergence tolerance: 4.39e-06> pairs(fm1Oran.nlme)> fm2Oran.nlme <- update(fm1Oran.nlme, random = Asym ~ 1)> anova(fm1Oran.nlme, fm2Oran.nlme)Model df AIC BIC logLik Test L.Ratio p-valuefm1Oran.nlme 1 10 280 296 -130fm2Oran.nlme 2 5 273 281 -132 1 vs 2 3.19 0.671> plot(fm1Oran.nlme)> ## level = 0:1 requests fixed (0) and within-group (1) predictions> plot(augPred(fm2Oran.nlme, level = 0:1),+ layout = c(5,1))> qqnorm(fm2Oran.nlme, abline = c(0,1))> (fm1Theo.nlme <- nlme(fm1Theo.lis))Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :Iteration 2, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!Nonlinear mixed-effects model fit by maximum likelihoodModel: conc ~ SSfol(Dose, Time, lKe, lKa, lCl)Data: TheophLog-likelihood: -173Fixed: list(lKe ~ 1, lKa ~ 1, lCl ~ 1)lKe lKa lCl-2.433 0.451 -3.214Random effects:Formula: list(lKe ~ 1, lKa ~ 1, lCl ~ 1)Level: SubjectStructure: General positive-definite, Log-Cholesky parametrizationStdDev CorrlKe 0.131 lKe lKalKa 0.638 0.012lCl 0.251 0.995 -0.089Residual 0.682Number of Observations: 132Number of Groups: 12> ## IGNORE_RDIFF_BEGIN> try( intervals(fm1Theo.nlme, which="var-cov") ) ## could fail: Non-positive definite...Approximate 95% confidence intervalsRandom Effects:Level: Subjectlower est. uppersd(lKe) 0.0574 0.1310 0.299sd(lKa) 0.3845 0.6378 1.058sd(lCl) 0.1557 0.2512 0.405cor(lKe,lKa) -0.9302 0.0116 0.933cor(lKe,lCl) -0.9950 0.9948 1.000cor(lKa,lCl) -0.7711 -0.0892 0.688Within-group standard error:lower est. upper0.596 0.682 0.780> ## IGNORE_RDIFF_END> (fm2Theo.nlme <- update(fm1Theo.nlme,+ random = pdDiag(lKe + lKa + lCl ~ 1)))Nonlinear mixed-effects model fit by maximum likelihoodModel: conc ~ SSfol(Dose, Time, lKe, lKa, lCl)Data: TheophLog-likelihood: -177Fixed: list(lKe ~ 1, lKa ~ 1, lCl ~ 1)lKe lKa lCl-2.455 0.466 -3.227Random effects:Formula: list(lKe ~ 1, lKa ~ 1, lCl ~ 1)Level: SubjectStructure: DiagonallKe lKa lCl ResidualStdDev: 1.93e-05 0.644 0.167 0.709Number of Observations: 132Number of Groups: 12> fm3Theo.nlme <-+ update(fm2Theo.nlme, random = pdDiag(lKa + lCl ~ 1))> anova(fm1Theo.nlme, fm3Theo.nlme, fm2Theo.nlme)Model df AIC BIC logLik Test L.Ratio p-valuefm1Theo.nlme 1 10 367 395 -173fm3Theo.nlme 2 6 366 383 -177 1 vs 2 7.4 0.116fm2Theo.nlme 3 7 368 388 -177 2 vs 3 0.0 0.949> plot(fm3Theo.nlme)> qqnorm(fm3Theo.nlme, ~ ranef(.))> CO2Grouped Data: uptake ~ conc | PlantPlant Type Treatment conc uptake1 Qn1 Quebec nonchilled 95 16.02 Qn1 Quebec nonchilled 175 30.43 Qn1 Quebec nonchilled 250 34.84 Qn1 Quebec nonchilled 350 37.25 Qn1 Quebec nonchilled 500 35.36 Qn1 Quebec nonchilled 675 39.27 Qn1 Quebec nonchilled 1000 39.78 Qn2 Quebec nonchilled 95 13.69 Qn2 Quebec nonchilled 175 27.310 Qn2 Quebec nonchilled 250 37.111 Qn2 Quebec nonchilled 350 41.812 Qn2 Quebec nonchilled 500 40.613 Qn2 Quebec nonchilled 675 41.414 Qn2 Quebec nonchilled 1000 44.315 Qn3 Quebec nonchilled 95 16.216 Qn3 Quebec nonchilled 175 32.417 Qn3 Quebec nonchilled 250 40.318 Qn3 Quebec nonchilled 350 42.119 Qn3 Quebec nonchilled 500 42.920 Qn3 Quebec nonchilled 675 43.921 Qn3 Quebec nonchilled 1000 45.522 Qc1 Quebec chilled 95 14.223 Qc1 Quebec chilled 175 24.124 Qc1 Quebec chilled 250 30.325 Qc1 Quebec chilled 350 34.626 Qc1 Quebec chilled 500 32.527 Qc1 Quebec chilled 675 35.428 Qc1 Quebec chilled 1000 38.729 Qc2 Quebec chilled 95 9.330 Qc2 Quebec chilled 175 27.331 Qc2 Quebec chilled 250 35.032 Qc2 Quebec chilled 350 38.833 Qc2 Quebec chilled 500 38.634 Qc2 Quebec chilled 675 37.535 Qc2 Quebec chilled 1000 42.436 Qc3 Quebec chilled 95 15.137 Qc3 Quebec chilled 175 21.038 Qc3 Quebec chilled 250 38.139 Qc3 Quebec chilled 350 34.040 Qc3 Quebec chilled 500 38.941 Qc3 Quebec chilled 675 39.642 Qc3 Quebec chilled 1000 41.443 Mn1 Mississippi nonchilled 95 10.644 Mn1 Mississippi nonchilled 175 19.245 Mn1 Mississippi nonchilled 250 26.246 Mn1 Mississippi nonchilled 350 30.047 Mn1 Mississippi nonchilled 500 30.948 Mn1 Mississippi nonchilled 675 32.449 Mn1 Mississippi nonchilled 1000 35.550 Mn2 Mississippi nonchilled 95 12.051 Mn2 Mississippi nonchilled 175 22.052 Mn2 Mississippi nonchilled 250 30.653 Mn2 Mississippi nonchilled 350 31.854 Mn2 Mississippi nonchilled 500 32.455 Mn2 Mississippi nonchilled 675 31.156 Mn2 Mississippi nonchilled 1000 31.557 Mn3 Mississippi nonchilled 95 11.358 Mn3 Mississippi nonchilled 175 19.459 Mn3 Mississippi nonchilled 250 25.860 Mn3 Mississippi nonchilled 350 27.961 Mn3 Mississippi nonchilled 500 28.562 Mn3 Mississippi nonchilled 675 28.163 Mn3 Mississippi nonchilled 1000 27.864 Mc1 Mississippi chilled 95 10.565 Mc1 Mississippi chilled 175 14.966 Mc1 Mississippi chilled 250 18.167 Mc1 Mississippi chilled 350 18.968 Mc1 Mississippi chilled 500 19.569 Mc1 Mississippi chilled 675 22.270 Mc1 Mississippi chilled 1000 21.971 Mc2 Mississippi chilled 95 7.772 Mc2 Mississippi chilled 175 11.473 Mc2 Mississippi chilled 250 12.374 Mc2 Mississippi chilled 350 13.075 Mc2 Mississippi chilled 500 12.576 Mc2 Mississippi chilled 675 13.777 Mc2 Mississippi chilled 1000 14.478 Mc3 Mississippi chilled 95 10.679 Mc3 Mississippi chilled 175 18.080 Mc3 Mississippi chilled 250 17.981 Mc3 Mississippi chilled 350 17.982 Mc3 Mississippi chilled 500 17.983 Mc3 Mississippi chilled 675 18.984 Mc3 Mississippi chilled 1000 19.9> plot(CO2, outer = ~Treatment*Type, layout = c(4,1))> (fm1CO2.lis <- nlsList(SSasympOff, CO2))Call:Model: uptake ~ SSasympOff(conc, Asym, lrc, c0) | PlantData: CO2Coefficients:Asym lrc c0Qn1 38.1 -4.38 51.2Qn2 42.9 -4.67 55.9Qn3 44.2 -4.49 54.6Qc1 36.4 -4.86 31.1Qc3 40.7 -4.95 35.1Qc2 39.8 -4.46 72.1Mn3 28.5 -4.59 47.0Mn2 32.1 -4.47 56.0Mn1 34.1 -5.06 36.4Mc2 13.6 -4.56 13.1Mc3 18.5 -3.47 67.8Mc1 21.8 -5.14 -20.4Degrees of freedom: 84 total; 48 residualResidual standard error: 1.8> ## IGNORE_RDIFF_BEGIN> (fm1CO2.nlme <- nlme(fm1CO2.lis))Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!Nonlinear mixed-effects model fit by maximum likelihoodModel: uptake ~ SSasympOff(conc, Asym, lrc, c0)Data: CO2Log-likelihood: -201Fixed: list(Asym ~ 1, lrc ~ 1, c0 ~ 1)Asym lrc c032.47 -4.64 43.55Random effects:Formula: list(Asym ~ 1, lrc ~ 1, c0 ~ 1)Level: PlantStructure: General positive-definite, Log-Cholesky parametrizationStdDev CorrAsym 9.51 Asym lrclrc 0.13 -0.165c0 10.33 0.999 -0.133Residual 1.77Number of Observations: 84Number of Groups: 12> ## IGNORE_RDIFF_END> (fm2CO2.nlme <- update(fm1CO2.nlme, random = Asym + lrc ~ 1))Nonlinear mixed-effects model fit by maximum likelihoodModel: uptake ~ SSasympOff(conc, Asym, lrc, c0)Data: CO2Log-likelihood: -203Fixed: list(Asym ~ 1, lrc ~ 1, c0 ~ 1)Asym lrc c032.41 -4.56 49.34Random effects:Formula: list(Asym ~ 1, lrc ~ 1)Level: PlantStructure: General positive-definite, Log-Cholesky parametrizationStdDev CorrAsym 9.66 Asymlrc 0.20 -0.777Residual 1.81Number of Observations: 84Number of Groups: 12> anova(fm1CO2.nlme, fm2CO2.nlme)Model df AIC BIC logLik Test L.Ratio p-valuefm1CO2.nlme 1 10 423 447 -201fm2CO2.nlme 2 7 420 437 -203 1 vs 2 2.9 0.408> plot(fm2CO2.nlme,id = 0.05,cex = 0.8,adj = -0.5)> fm2CO2.nlmeRE <- ranef(fm2CO2.nlme, augFrame = TRUE)> fm2CO2.nlmeREAsym lrc Type Treatment conc uptakeQn1 6.172 0.04836 Quebec nonchilled 435 33.2Qn2 10.533 -0.17284 Quebec nonchilled 435 35.2Qn3 12.218 -0.05799 Quebec nonchilled 435 37.6Qc1 3.352 -0.07559 Quebec chilled 435 30.0Qc3 7.474 -0.19242 Quebec chilled 435 32.6Qc2 7.928 -0.18032 Quebec chilled 435 32.7Mn3 -4.073 0.03345 Mississippi nonchilled 435 24.1Mn2 -0.142 0.00565 Mississippi nonchilled 435 27.3Mn1 0.241 -0.19386 Mississippi nonchilled 435 26.4Mc2 -18.799 0.31937 Mississippi chilled 435 12.1Mc3 -13.117 0.29943 Mississippi chilled 435 17.3Mc1 -11.787 0.16676 Mississippi chilled 435 18.0> class(fm2CO2.nlmeRE)[1] "ranef.lme" "data.frame"> plot(fm2CO2.nlmeRE, form = ~ Type * Treatment)> contrasts(CO2$Type)[,1]Quebec -1Mississippi 1> contrasts(CO2$Treatment)[,1]nonchilled -1chilled 1> fm3CO2.nlme <- update(fm2CO2.nlme,+ fixed = list(Asym ~ Type * Treatment, lrc + c0 ~ 1),+ start = c(32.412, 0, 0, 0, -4.5603, 49.344))> summary(fm3CO2.nlme)Nonlinear mixed-effects model fit by maximum likelihoodModel: uptake ~ SSasympOff(conc, Asym, lrc, c0)Data: CO2AIC BIC logLik394 418 -187Random effects:Formula: list(Asym ~ 1, lrc ~ 1)Level: PlantStructure: General positive-definite, Log-Cholesky parametrizationStdDev CorrAsym.(Intercept) 2.930 As.(I)lrc 0.164 -0.906Residual 1.850Fixed effects: list(Asym ~ Type * Treatment, lrc + c0 ~ 1)Value Std.Error DF t-value p-valueAsym.(Intercept) 32.4 0.94 67 34.7 0.0000Asym.Type1 -7.1 0.60 67 -11.9 0.0000Asym.Treatment1 -3.8 0.59 67 -6.5 0.0000Asym.Type1:Treatment1 -1.2 0.59 67 -2.0 0.0462lrc -4.6 0.08 67 -54.1 0.0000c0 49.5 4.46 67 11.1 0.0000Correlation:As.(I) Asym.Ty1 Asym.Tr1 A.T1:T lrcAsym.Type1 -0.044Asym.Treatment1 -0.021 0.151Asym.Type1:Treatment1 -0.023 0.161 0.225lrc -0.660 0.202 0.113 0.132c0 -0.113 0.060 0.018 0.063 0.653Standardized Within-Group Residuals:Min Q1 Med Q3 Max-2.8929 -0.4616 -0.0328 0.5208 2.8877Number of Observations: 84Number of Groups: 12> anova(fm3CO2.nlme, Terms = 2:4)F-test for: Asym.Type, Asym.Treatment, Asym.Type:TreatmentnumDF denDF F-value p-value1 3 67 54.8 <.0001> fm3CO2.nlmeRE <- ranef(fm3CO2.nlme, aug = TRUE)> plot(fm3CO2.nlmeRE, form = ~ Type * Treatment)> fm3CO2.fix <- fixef(fm3CO2.nlme)> fm4CO2.nlme <- update(fm3CO2.nlme,+ fixed = list(Asym + lrc ~ Type * Treatment, c0 ~ 1),+ start = c(fm3CO2.fix[1:5], 0, 0, 0, fm3CO2.fix[6]))Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!> ## IGNORE_RDIFF_BEGIN> summary(fm4CO2.nlme)Nonlinear mixed-effects model fit by maximum likelihoodModel: uptake ~ SSasympOff(conc, Asym, lrc, c0)Data: CO2AIC BIC logLik388 420 -181Random effects:Formula: list(Asym ~ 1, lrc ~ 1)Level: PlantStructure: General positive-definite, Log-Cholesky parametrizationStdDev CorrAsym.(Intercept) 2.3496 As.(I)lrc.(Intercept) 0.0796 -0.92Residual 1.7920Fixed effects: list(Asym + lrc ~ Type * Treatment, c0 ~ 1)Value Std.Error DF t-value p-valueAsym.(Intercept) 32.3 0.78 64 41.2 0.0000Asym.Type1 -8.0 0.78 64 -10.3 0.0000Asym.Treatment1 -4.2 0.78 64 -5.4 0.0000Asym.Type1:Treatment1 -2.7 0.78 64 -3.5 0.0008lrc.(Intercept) -4.5 0.08 64 -55.7 0.0000lrc.Type1 0.1 0.06 64 2.4 0.0185lrc.Treatment1 0.1 0.06 64 1.8 0.0746lrc.Type1:Treatment1 0.2 0.06 64 3.3 0.0014c0 50.5 4.36 64 11.6 0.0000Correlation:As.(I) Asym.Ty1 Asym.Tr1 A.T1:T lr.(I)Asym.Type1 -0.017Asym.Treatment1 -0.010 -0.017Asym.Type1:Treatment1 -0.020 -0.006 -0.011lrc.(Intercept) -0.471 0.004 0.001 0.009lrc.Type1 -0.048 -0.548 -0.005 -0.018 0.402lrc.Treatment1 -0.031 -0.004 -0.551 -0.033 0.322lrc.Type1:Treatment1 -0.026 -0.015 -0.032 -0.547 0.351c0 -0.133 0.038 0.020 0.019 0.735lrc.Ty1 lrc.Tr1 l.T1:TAsym.Type1Asym.Treatment1Asym.Type1:Treatment1lrc.(Intercept)lrc.Type1lrc.Treatment1 0.375lrc.Type1:Treatment1 0.395 0.487c0 0.104 0.083 0.140Standardized Within-Group Residuals:Min Q1 Med Q3 Max-2.8621 -0.4944 -0.0422 0.5661 3.0405Number of Observations: 84Number of Groups: 12> ## IGNORE_RDIFF_END> fm5CO2.nlme <- update(fm4CO2.nlme, random = Asym ~ 1)> anova(fm4CO2.nlme, fm5CO2.nlme)Model df AIC BIC logLik Test L.Ratio p-valuefm4CO2.nlme 1 13 388 420 -181fm5CO2.nlme 2 11 387 414 -182 1 vs 2 2.64 0.268> CO2$type <- 2 * (as.integer(CO2$Type) - 1.5)> CO2$treatment <- 2 * (as.integer(CO2$Treatment) - 1.5)> fm1CO2.nls <- nls(uptake ~ SSasympOff(conc, Asym.Intercept ++ Asym.Type * type + Asym.Treatment * treatment ++ Asym.TypeTreatment * type * treatment, lrc.Intercept ++ lrc.Type * type + lrc.Treatment * treatment ++ lrc.TypeTreatment * type * treatment, c0), data = CO2,+ start = c(Asym.Intercept = 32.371, Asym.Type = -8.0086,+ Asym.Treatment = -4.2001, Asym.TypeTreatment = -2.7253,+ lrc.Intercept = -4.5267, lrc.Type = 0.13112,+ lrc.Treatment = 0.093928, lrc.TypeTreatment = 0.17941,+ c0 = 50.126))> anova(fm5CO2.nlme, fm1CO2.nls)Model df AIC BIC logLik Test L.Ratio p-valuefm5CO2.nlme 1 11 387 414 -182fm1CO2.nls 2 10 418 443 -199 1 vs 2 33.3 <.0001> # plot(augPred(fm5CO2.nlme, level = 0:1), ## FIXME: problem with levels> # layout = c(6,2)) ## Actually a problem with contrasts.> ## This fit just ping-pongs.> #fm1Quin.nlme <-> # nlme(conc ~ quinModel(Subject, time, conc, dose, interval,> # lV, lKa, lCl),> # data = Quinidine, fixed = lV + lKa + lCl ~ 1,> # random = pdDiag(lV + lCl ~ 1), groups = ~ Subject,> # start = list(fixed = c(5, -0.3, 2)),> # na.action = NULL, naPattern = ~ !is.na(conc), verbose = TRUE)> #fm1Quin.nlme> #fm1Quin.nlmeRE <- ranef(fm1Quin.nlme, aug = TRUE)> #fm1Quin.nlmeRE[1:3,]> # plot(fm1Quin.nlmeRE, form = lCl ~ Age + Smoke + Ethanol + ## FIXME: problem in max> # Weight + Race + Height + glyco + Creatinine + Heart,> # control = list(cex.axis = 0.7))> #fm1Quin.fix <- fixef(fm1Quin.nlme)> #fm2Quin.nlme <- update(fm1Quin.nlme,> # fixed = list(lCl ~ glyco, lKa + lV ~ 1),> # start = c(fm1Quin.fix[3], 0, fm1Quin.fix[2:1]))> fm2Quin.nlme <-+ nlme(conc ~ quinModel(Subject, time, conc, dose, interval,+ lV, lKa, lCl),+ data = Quinidine, fixed = list(lCl ~ glyco, lV + lKa ~ 1),+ random = pdDiag(diag(c(0.3,0.3)), form = lV + lCl ~ 1),+ groups = ~ Subject,+ start = list(fixed = c(2.5, 0, 5.4, -0.2)),+ na.action = NULL, naPattern = ~ !is.na(conc))> summary(fm2Quin.nlme) # wrong valuesNonlinear mixed-effects model fit by maximum likelihoodModel: conc ~ quinModel(Subject, time, conc, dose, interval, lV, lKa, lCl)Data: QuinidineAIC BIC logLik892 919 -439Random effects:Formula: list(lV ~ 1, lCl ~ 1)Level: SubjectStructure: DiagonallV lCl.(Intercept) ResidualStdDev: 0.000263 0.271 0.651Fixed effects: list(lCl ~ glyco, lV + lKa ~ 1)Value Std.Error DF t-value p-valuelCl.(Intercept) 3.12 0.0655 222 47.7 0.000lCl.glyco -0.50 0.0428 222 -11.7 0.000lV 5.27 0.0948 222 55.6 0.000lKa -0.84 0.3039 222 -2.8 0.006Correlation:lC.(I) lCl.gl lVlCl.glyco -0.880lV -0.072 0.027lKa -0.272 0.149 0.538Standardized Within-Group Residuals:Min Q1 Med Q3 Max-2.5458 -0.5342 -0.0221 0.5053 3.5016Number of Observations: 361Number of Groups: 136> options(contrasts = c("contr.treatment", "contr.poly"))> fm2Quin.fix <- fixef(fm2Quin.nlme)> ## subsequent fits don't work> #fm3Quin.nlme <- update(fm2Quin.nlme,> # fixed = list(lCl ~ glyco + Creatinine, lKa + lV ~ 1),> # start = c(fm2Quin.fix[1:2], 0.2, fm2Quin.fix[3:4]))> #summary(fm3Quin.nlme)> #fm3Quin.fix <- fixef(fm3Quin.nlme)> #fm4Quin.nlme <- update(fm3Quin.nlme,> # fixed = list(lCl ~ glyco + Creatinine + Weight, lKa + lV ~ 1),> # start = c(fm3Quin.fix[1:3], 0, fm3Quin.fix[4:5]))> #summary(fm4Quin.nlme)> ## This fit just ping-pongs> ##fm1Wafer.nlmeR <-> ## nlme(current ~ A + B * cos(4.5679 * voltage) +> ## C * sin(4.5679 * voltage), data = Wafer,> ## fixed = list(A ~ voltage + I(voltage^2), B + C ~ 1),> ## random = list(Wafer = A ~ voltage + I(voltage^2),> ## Site = pdBlocked(list(A~1, A~voltage+I(voltage^2)-1))),> ### start = fixef(fm4Wafer), method = "REML", control = list(tolerance=1e-2))> ## start = c(-4.255, 5.622, 1.258, -0.09555, 0.10434),> ## method = "REML", control = list(tolerance = 1e-2))> ##fm1Wafer.nlmeR> ##fm1Wafer.nlme <- update(fm1Wafer.nlmeR, method = "ML")>> (fm2Wafer.nlme <-+ nlme(current ~ A + B * cos(w * voltage + pi/4),+ data = Wafer,+ fixed = list(A ~ voltage + I(voltage^2), B + w ~ 1),+ random = list(Wafer = pdDiag(list(A ~ voltage + I(voltage^2), B + w ~ 1)),+ Site = pdDiag(list(A ~ voltage+I(voltage^2), B ~ 1))),+ start = c(-4.255, 5.622, 1.258, -0.09555, 4.5679)))Nonlinear mixed-effects model fit by maximum likelihoodModel: current ~ A + B * cos(w * voltage + pi/4)Data: WaferLog-likelihood: 663Fixed: list(A ~ voltage + I(voltage^2), B + w ~ 1)A.(Intercept) A.voltage A.I(voltage^2) B-4.265 5.633 1.256 -0.141w4.593Random effects:Formula: list(A ~ voltage + I(voltage^2), B ~ 1, w ~ 1)Level: WaferStructure: DiagonalA.(Intercept) A.voltage A.I(voltage^2) B wStdDev: 0.127 0.337 0.0488 0.00506 5.44e-05Formula: list(A ~ voltage + I(voltage^2), B ~ 1)Level: Site %in% WaferStructure: DiagonalA.(Intercept) A.voltage A.I(voltage^2) B ResidualStdDev: 0.0618 0.269 0.0559 4.46e-06 0.00786Number of Observations: 400Number of Groups:Wafer Site %in% Wafer10 80> plot(fm2Wafer.nlme, resid(.) ~ voltage | Wafer,+ panel = function(x, y, ...) {+ panel.grid()+ panel.xyplot(x, y)+ panel.loess(x, y, lty = 2)+ panel.abline(0, 0)+ })> ## anova(fm1Wafer.nlme, fm2Wafer.nlme, test = FALSE)> # intervals(fm2Wafer.nlme)>> # 8.3 Extending the Basic nlme Model>> #fm4Theo.nlme <- update(fm3Theo.nlme,> # weights = varConstPower(power = 0.1))> # this fit is way off> #fm4Theo.nlme> #anova(fm3Theo.nlme, fm4Theo.nlme)> #plot(fm4Theo.nlme)> ## xlim used to hide an unusually high fitted value and enhance> ## visualization of the heteroscedastic pattern> # plot(fm4Quin.nlme, xlim = c(0, 6.2))> #fm5Quin.nlme <- update(fm4Quin.nlme, weights = varPower())> #summary(fm5Quin.nlme)> #anova(fm4Quin.nlme, fm5Quin.nlme)> #plot(fm5Quin.nlme, xlim = c(0, 6.2))> var.nlme <- nlme(follicles ~ A + B * sin(2 * pi * w * Time) ++ C * cos(2 * pi * w *Time), data = Ovary,+ fixed = A + B + C + w ~ 1, random = pdDiag(A + B + w ~ 1),+ # start = c(fixef(fm5Ovar.lme), 1))+ start = c(12.18, -3.298, -0.862, 1))> ##fm1Ovar.nlme> ##ACF(fm1Ovar.nlme)> ##plot(ACF(fm1Ovar.nlme, maxLag = 10), alpha = 0.05)> ##fm2Ovar.nlme <- update(fm1Ovar.nlme, correlation = corAR1(0.311))> ##fm3Ovar.nlme <- update(fm1Ovar.nlme, correlation = corARMA(p=0, q=2))> ##anova(fm2Ovar.nlme, fm3Ovar.nlme, test = FALSE)> ##intervals(fm2Ovar.nlme)> ##fm4Ovar.nlme <- update(fm2Ovar.nlme, random = A ~ 1)> ##anova(fm2Ovar.nlme, fm4Ovar.nlme)> ##if (interactive()) fm5Ovar.nlme <- update(fm4Ovar.nlme, correlation = corARMA(p=1, q=1))> # anova(fm4Ovar.nlme, fm5Ovar.nlme)> # plot(ACF(fm5Ovar.nlme, maxLag = 10, resType = "n"),> # alpha = 0.05)> # fm5Ovar.lmeML <- update(fm5Ovar.lme, method = "ML")> # intervals(fm5Ovar.lmeML)> # fm6Ovar.lmeML <- update(fm5Ovar.lmeML, random = ~1)> # anova(fm5Ovar.lmeML, fm6Ovar.lmeML)> # anova(fm6Ovar.lmeML, fm5Ovar.nlme)> # intervals(fm5Ovar.nlme, which = "fixed")> fm1Dial.lis <-+ nlsList(rate ~ SSasympOff(pressure, Asym, lrc, c0) | QB,+ data = Dialyzer)> fm1Dial.lisCall:Model: rate ~ SSasympOff(pressure, Asym, lrc, c0) | QBData: DialyzerCoefficients:Asym lrc c0200 45.0 0.765 0.224300 62.2 0.253 0.225Degrees of freedom: 140 total; 134 residualResidual standard error: 3.8> plot(intervals(fm1Dial.lis))> fm1Dial.gnls <- gnls(rate ~ SSasympOff(pressure, Asym, lrc, c0),+ data = Dialyzer, params = list(Asym + lrc ~ QB, c0 ~ 1),+ start = c(53.6, 8.6, 0.51, -0.26, 0.225))> fm1Dial.gnlsGeneralized nonlinear least squares fitModel: rate ~ SSasympOff(pressure, Asym, lrc, c0)Data: DialyzerLog-likelihood: -383Coefficients:Asym.(Intercept) Asym.QB300 lrc.(Intercept)44.986 17.240 0.766lrc.QB300 c0-0.514 0.224Degrees of freedom: 140 total; 135 residualResidual standard error: 3.79> Dialyzer$QBcontr <- 2 * (Dialyzer$QB == 300) - 1> fm1Dial.nls <-+ nls(rate ~ SSasympOff(pressure, Asym.Int + Asym.QB * QBcontr,+ lrc.Int + lrc.QB * QBcontr, c0), data = Dialyzer,+ start = c(Asym.Int = 53.6, Asym.QB = 8.6, lrc.Int = 0.51,+ lrc.QB = -0.26, c0 = 0.225))> ## IGNORE_RDIFF_BEGIN> summary(fm1Dial.nls)Formula: rate ~ SSasympOff(pressure, Asym.Int + Asym.QB * QBcontr, lrc.Int +lrc.QB * QBcontr, c0)Parameters:Estimate Std. Error t value Pr(>|t|)Asym.Int 53.6065 0.7054 75.99 < 2e-16 ***Asym.QB 8.6201 0.6792 12.69 < 2e-16 ***lrc.Int 0.5087 0.0552 9.21 5.5e-16 ***lrc.QB -0.2568 0.0450 -5.70 7.0e-08 ***c0 0.2245 0.0106 21.13 < 2e-16 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 3.79 on 135 degrees of freedomNumber of iterations to convergence: 4Achieved convergence tolerance: 7.24e-06> ## IGNORE_RDIFF_END> logLik(fm1Dial.nls)'log Lik.' -383 (df=6)> plot(fm1Dial.gnls, resid(.) ~ pressure, abline = 0)> fm2Dial.gnls <- update(fm1Dial.gnls,+ weights = varPower(form = ~ pressure))> anova(fm1Dial.gnls, fm2Dial.gnls)Model df AIC BIC logLik Test L.Ratio p-valuefm1Dial.gnls 1 6 777 795 -383fm2Dial.gnls 2 7 748 769 -367 1 vs 2 30.8 <.0001> ACF(fm2Dial.gnls, form = ~ 1 | Subject)lag ACF1 0 1.000002 1 0.715673 2 0.504544 3 0.294815 4 0.209756 5 0.138577 6 -0.00202> plot(ACF(fm2Dial.gnls, form = ~ 1 | Subject), alpha = 0.05)> fm3Dial.gnls <-+ update(fm2Dial.gnls, corr = corAR1(0.716, form = ~ 1 | Subject))> fm3Dial.gnlsGeneralized nonlinear least squares fitModel: rate ~ SSasympOff(pressure, Asym, lrc, c0)Data: DialyzerLog-likelihood: -323Coefficients:Asym.(Intercept) Asym.QB300 lrc.(Intercept)46.911 16.400 0.542lrc.QB300 c0-0.339 0.215Correlation Structure: AR(1)Formula: ~1 | SubjectParameter estimate(s):Phi0.744Variance function:Structure: Power of variance covariateFormula: ~pressureParameter estimates:power0.572Degrees of freedom: 140 total; 135 residualResidual standard error: 3.18> intervals(fm3Dial.gnls)Approximate 95% confidence intervalsCoefficients:lower est. upperAsym.(Intercept) 43.877 46.911 49.945Asym.QB300 11.633 16.400 21.167lrc.(Intercept) 0.435 0.542 0.648lrc.QB300 -0.487 -0.339 -0.192c0 0.206 0.215 0.223Correlation structure:lower est. upperPhi 0.622 0.744 0.831Variance function:lower est. upperpower 0.443 0.572 0.702Residual standard error:lower est. upper2.59 3.13 3.77> anova(fm2Dial.gnls, fm3Dial.gnls)Model df AIC BIC logLik Test L.Ratio p-valuefm2Dial.gnls 1 7 748 769 -367fm3Dial.gnls 2 8 661 685 -323 1 vs 2 89.4 <.0001> # restore two fitted models> fm2Dial.lme <-+ lme(rate ~(pressure + I(pressure^2) + I(pressure^3) + I(pressure^4))*QB,+ Dialyzer, ~ pressure + I(pressure^2),+ weights = varPower(form = ~ pressure))> fm2Dial.lmeML <- update(fm2Dial.lme, method = "ML")> fm3Dial.gls <-+ gls(rate ~(pressure + I(pressure^2) + I(pressure^3) + I(pressure^4))*QB,+ Dialyzer, weights = varPower(form = ~ pressure),+ corr = corAR1(0.771, form = ~ 1 | Subject))> fm3Dial.glsML <- update(fm3Dial.gls, method = "ML")> anova( fm2Dial.lmeML, fm3Dial.glsML, fm3Dial.gnls, test = FALSE)Model df AIC BIC logLikfm2Dial.lmeML 1 18 652 705 -308fm3Dial.glsML 2 13 648 686 -311fm3Dial.gnls 3 8 661 685 -323> # cleanup>> summary(warnings())No warnings>> proc.time()user system elapsed61.627 0.112 61.765