R Under development (unstable) (2024-08-14 r87014) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu R 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 installed "mlbook" chapters > testdir <- system.file("mlbook", 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(1) + options(warn = 1, digits = 5) + source(file.path(testdir, f), echo = TRUE, + max.deparse.length = Inf, keep.source = TRUE) + } ====== ch04.R ====== > library(nlme) > # data(bdf) > ## Fit the null model > ## Compare with Table 4.1, p. 47 > fm1 <- lme(langPOST ~ 1, data = bdf, random = ~ 1 | schoolNR) > VarCorr(fm1) schoolNR = pdLogChol(1) Variance StdDev (Intercept) 19.633 4.4309 Residual 64.564 8.0352 > -2*c(logLik(fm1)) # deviance [1] 16253 > ## Fit model with fixed IQ term and random intercept > ## Compare with Table 4.2, p. 49 > ## From the results in Tables 4.2 and 4.4, it appears that > ## maximum likelihood fits are used, not REML fits. > fm2 <- update(fm1, langPOST ~ IQ.ver.cen) > summary(fm2) Linear mixed-effects model fit by REML Data: bdf AIC BIC logLik 15264 15287 -7627.9 Random effects: Formula: ~1 | schoolNR (Intercept) Residual StdDev: 3.0987 6.4996 Fixed effects: langPOST ~ IQ.ver.cen Value Std.Error DF t-value p-value (Intercept) 40.608 0.308186 2155 131.766 0 IQ.ver.cen 2.488 0.070081 2155 35.496 0 Correlation: (Intr) IQ.ver.cen 0.018 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -4.093938 -0.637456 0.057947 0.706070 3.144829 Number of Observations: 2287 Number of Groups: 131 > VarCorr(fm2) schoolNR = pdLogChol(1) Variance StdDev (Intercept) 9.6017 3.0987 Residual 42.2445 6.4996 > -2 * c(logLik(fm2)) # deviance [1] 15256 > ## Purely fixed-effects model for comparison > ## Compare with Table 4.3, p. 51 > fm3 <- lm(langPOST ~ IQ.ver.cen, data = bdf) > summary(fm3) Call: lm(formula = langPOST ~ IQ.ver.cen, data = bdf) Residuals: Min 1Q Median 3Q Max -28.702 -4.394 0.606 5.260 26.221 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 40.9348 0.1492 274.3 <2e-16 *** IQ.ver.cen 2.6539 0.0722 36.8 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 7.14 on 2285 degrees of freedom Multiple R-squared: 0.372, Adjusted R-squared: 0.372 F-statistic: 1.35e+03 on 1 and 2285 DF, p-value: <2e-16 > -2 * c(logLik(fm3)) # deviance [1] 15478 > ## Model with average IQ for the school > ## Compare with Table 4.4, p. 55 > fm4 <- update(fm2, langPOST ~ IQ.ver.cen + avg.IQ.ver.cen) > summary(fm4) Linear mixed-effects model fit by REML Data: bdf AIC BIC logLik 15242 15271 -7616.1 Random effects: Formula: ~1 | schoolNR (Intercept) Residual StdDev: 2.8082 6.494 Fixed effects: langPOST ~ IQ.ver.cen + avg.IQ.ver.cen Value Std.Error DF t-value p-value (Intercept) 40.741 0.286595 2155 142.155 0 IQ.ver.cen 2.415 0.071676 2155 33.690 0 avg.IQ.ver.cen 1.589 0.314772 129 5.049 0 Correlation: (Intr) IQ.vr. IQ.ver.cen 0.000 avg.IQ.ver.cen 0.077 -0.228 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -4.130749 -0.642218 0.060944 0.702342 3.152993 Number of Observations: 2287 Number of Groups: 131 > VarCorr(fm4) schoolNR = pdLogChol(1) Variance StdDev (Intercept) 7.8859 2.8082 Residual 42.1723 6.4940 > -2 * c(logLik(fm4)) # deviance [1] 15232 ====== ch05.R ====== > library(nlme) > # data(bdf) > ## Model with random slope for IQ.ver.cen > ## Compare with Table 5.1, p. 71. > fm5 <- lme(langPOST ~ IQ.ver.cen + avg.IQ.ver.cen, + data = bdf, random = ~ IQ.ver.cen, method = "ML") > summary(fm5) Linear mixed-effects model fit by maximum likelihood Data: bdf AIC BIC logLik 15228 15268 -7606.8 Random effects: Formula: ~IQ.ver.cen | schoolNR Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 2.81410 (Intr) IQ.ver.cen 0.44728 -0.652 Residual 6.43043 Fixed effects: langPOST ~ IQ.ver.cen + avg.IQ.ver.cen Value Std.Error DF t-value p-value (Intercept) 40.750 0.28610 2155 142.433 0 IQ.ver.cen 2.459 0.08324 2155 29.541 0 avg.IQ.ver.cen 1.405 0.32168 129 4.368 0 Correlation: (Intr) IQ.vr. IQ.ver.cen -0.274 avg.IQ.ver.cen 0.028 -0.214 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -4.17512 -0.63982 0.06693 0.70462 2.71089 Number of Observations: 2287 Number of Groups: 131 > VarCorr(fm5) schoolNR = pdLogChol(IQ.ver.cen) Variance StdDev Corr (Intercept) 7.91916 2.81410 (Intr) IQ.ver.cen 0.20006 0.44728 -0.652 Residual 41.35049 6.43043 > -2 * c(logLik(fm5)) # deviance [1] 15214 > ## Add centered class size and interaction > ## Compare with Table 5.2, p. 75 > fm6 <- update(fm5, langPOST ~ avg.IQ.ver.cen + IQ.ver.cen * grpSiz.cen) > summary(fm6) Linear mixed-effects model fit by maximum likelihood Data: bdf AIC BIC logLik 15226 15278 -7604.2 Random effects: Formula: ~IQ.ver.cen | schoolNR Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 2.76882 (Intr) IQ.ver.cen 0.42132 -0.658 Residual 6.43136 Fixed effects: langPOST ~ avg.IQ.ver.cen + IQ.ver.cen + grpSiz.cen + IQ.ver.cen:grpSiz.cen Value Std.Error DF t-value p-value (Intercept) 40.893 0.29249 2153 139.809 0.0000 avg.IQ.ver.cen 1.246 0.32642 129 3.818 0.0002 IQ.ver.cen 2.443 0.08233 2153 29.674 0.0000 grpSiz.cen 0.057 0.03691 2153 1.556 0.1198 IQ.ver.cen:grpSiz.cen -0.022 0.01091 2153 -1.989 0.0468 Correlation: (Intr) a.IQ.. IQ.vr. grpSz. avg.IQ.ver.cen -0.024 IQ.ver.cen -0.276 -0.195 grpSiz.cen 0.249 -0.175 -0.086 IQ.ver.cen:grpSiz.cen -0.118 0.169 0.032 -0.233 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -4.163071 -0.639694 0.063419 0.710478 2.687724 Number of Observations: 2287 Number of Groups: 131 > VarCorr(fm6) schoolNR = pdLogChol(IQ.ver.cen) Variance StdDev Corr (Intercept) 7.66639 2.76882 (Intr) IQ.ver.cen 0.17751 0.42132 -0.658 Residual 41.36242 6.43136 > -2 * c(logLik(fm6)) # deviance [1] 15208 > > proc.time() user system elapsed 0.846 1.606 0.491