R Under development (unstable) (2026-01-12 r89300 ucrt) -- "Unsuffered Consequences" Copyright (C) 2026 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 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. > # test_compare_sas.R > library(lmerTest) Loading required package: lme4 Loading required package: Matrix Attaching package: 'lmerTest' The following object is masked from 'package:lme4': lmer The following object is masked from 'package:stats': step > > # WRE says "using if(requireNamespace("pkgname")) is preferred, if possible." > # even in tests: > assertError <- function(expr, ...) + if(requireNamespace("tools")) tools::assertError(expr, ...) else invisible() > assertWarning <- function(expr, ...) + if(requireNamespace("tools")) tools::assertWarning(expr, ...) else invisible() > > ##################################################################### > > > # Use contrasts to get particular estimates for the summary table: > l <- list(Frequency="contr.SAS", Income="contr.SAS") > m.carrots <- lmer(Preference ~ sens2*Frequency*Income + +(1+sens2|Consumer), data=carrots, contrasts=l) fixed-effect model matrix is rank deficient so dropping 12 columns / coefficients > an.m <- anova(m.carrots) Missing cells for: Frequency3:Income1, Frequency5:Income1, Frequency4:Income2, Frequency3:Income4, Frequency4:Income4, Frequency5:Income4, sens2:Frequency3:Income1, sens2:Frequency5:Income1, sens2:Frequency4:Income2, sens2:Frequency3:Income4, sens2:Frequency4:Income4, sens2:Frequency5:Income4. Interpret type III hypotheses with care. > > TOL <- 1e-4 > TOL2 <- 1e-5 > # with 4 decimals should agree with SAS output > # numbers before decimals should agree with SAS output > stopifnot( + all.equal(an.m[,"Pr(>F)"], + c(2e-5, 0.15512, 0.06939, 0.08223, 0.52459, 0.03119, 0.48344), + tolerance = TOL), + all.equal(round(an.m$DenDF), c(83, 83, 83, 83, 83, 83, 83)) + ) > > sm <- summary(m.carrots) > stopifnot( + isTRUE(all.equal(sm$coefficients[,"Pr(>|t|)"], + c(1e-10, 0.005061, 0.6865554, 0.342613, 0.129157, + 0.088231, 0.846000, 0.354472, 0.526318, 0.020646, 0.010188, + 0.031242, 0.055356, 0.694689, 0.099382, 0.28547, + 0.977774, 0.855653, 0.427737, 0.321086, 0.417465 , 0.204385, 0.784437, + 0.681434, 0.106180, 0.149122, 0.390870, 0.273686), tolerance=TOL, + check.attributes = FALSE)) + ) > > # Takes too long to run: > # if(requireNamespace("pbkrtest", quietly = TRUE)) { > # sm.kr <- summary(m.carrots, ddf = "Kenward-Roger") > # > # ## coefficients for Sat and KR agree in this example > # # cbind(sm$coefficients[,"Pr(>|t|)"], sm.kr$coefficients[,"Pr(>|t|)"]) > # all.equal(sm$coefficients[,"Pr(>|t|)"], sm.kr$coefficients[,"Pr(>|t|)"], > # tol=TOL) > # } > > ################################################################################ > ## checking lsmeans and difflsmeans > ## compare with SAS output > m <- lmer(Informed.liking ~ Product*Information*Gender + + (1|Product:Consumer) + (1|Consumer) , data=ham) > > > lsm <- lsmeansLT(m, which = "Product") > # head(lsm) > > stopifnot( + isTRUE(all.equal(lsm[, "Estimate"], c(5.8084, 5.1012, 6.0909, 5.9256), + tol=TOL, check.attributes = FALSE)), + isTRUE(all.equal(round(lsm[, "t value"], 2), c(24.93, 21.89, 26.14, 25.43), tolerance=TOL, + check.attributes = FALSE)), + isTRUE(all.equal(lsm[, "lower"], c(5.3499, 4.6428, 5.6324, 5.4672), tolerance=TOL, + check.attributes = FALSE)), + isTRUE(all.equal(lsm[, "upper"], c(6.2668, 5.5597, 6.5493, 6.3840), tolerance=TOL, + check.attributes = FALSE)) + ) > > ################################################################################ > # Not actually 'hard-coded' tests versus SAS results... > > m.carrots <- lmer(Preference ~ 0 + sens2 + Homesize + + (1+sens2 | Consumer), data=carrots, + control=lmerControl(optimizer="bobyqa")) > summary(m.carrots) Linear mixed model fit by REML. t-tests use Satterthwaite's method [ lmerModLmerTest] Formula: Preference ~ 0 + sens2 + Homesize + (1 + sens2 | Consumer) Data: carrots Control: lmerControl(optimizer = "bobyqa") REML criterion at convergence: 3748.9 Scaled residuals: Min 1Q Median 3Q Max -3.5322 -0.5571 0.0308 0.6297 2.8552 Random effects: Groups Name Variance Std.Dev. Corr Consumer (Intercept) 0.195168 0.44178 sens2 0.002779 0.05271 0.18 Residual 1.070441 1.03462 Number of obs: 1233, groups: Consumer, 103 Fixed effects: Estimate Std. Error df t value Pr(>|t|) sens2 7.068e-02 9.545e-03 1.020e+02 7.404 3.89e-11 *** Homesize1 4.910e+00 7.056e-02 1.013e+02 69.586 < 2e-16 *** Homesize3 4.661e+00 7.850e-02 1.013e+02 59.374 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation of Fixed Effects: sens2 Homsz1 Homesize1 0.061 Homesize3 0.055 0.003 > > (an.1 <- anova(m.carrots, type=1)) Type I Analysis of Variance Table with Satterthwaite's method Sum Sq Mean Sq NumDF DenDF F value Pr(>F) sens2 60.5 60.5 1 102.01 56.539 2.211e-11 *** Homesize 8927.2 4463.6 2 101.31 4169.863 < 2.2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > (an.3 <- anova(m.carrots)) Type III Analysis of Variance Table with Satterthwaite's method Sum Sq Mean Sq NumDF DenDF F value Pr(>F) sens2 58.7 58.7 1 102.01 54.821 3.892e-11 *** Homesize 8927.2 4463.6 2 101.31 4169.863 < 2.2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > (an.lme4 <- anova(m.carrots, ddf = "lme4")) # difference in SSQ MS and F-values Analysis of Variance Table npar Sum Sq Mean Sq F value sens2 1 0.0 0.0 0.0122 Homesize 2 8927.2 4463.6 4169.8633 > # Is this a problem with lme4? > # fm <- lm(Preference ~ 0 + sens2 + Homesize, data=carrots) > # anova(fm) > # coef(summary(fm)) > # Here the F value is a little greater than the squared t-value (as expected) > > stopifnot(all.equal(an.1[, "F value"], c(56.5394, 4169.87), tolerance = TOL2), + all.equal(an.3[, "F value"], c(54.8206, 4169.87), tolerance = TOL2)) > > > ################################################################################ > # Check exmaple from GLM SAS report > > ### example from the paper GLM SAS 101 report > a <- factor(c(1,1,1,2,2,2,2,2,1,2)) > b <- factor(c(1,1,2,1,2,2,2,2,2,1)) > f=factor(c(1,2,1,2,1,2,1,2,1,2)) > y <- c(12,14,11,20,17,23,35,46,15,16) > dd <- data.frame(a=a, b=b, y=y, f=f) > > ## check type 2 is order independent > model <- lmer(y ~ a*b + (1|f), data=dd) boundary (singular) fit: see help('isSingular') > model2 <- lmer(y ~ b*a + (1|f), data=dd) boundary (singular) fit: see help('isSingular') > (an <- anova(model, type=2)) Type II Analysis of Variance Table with Satterthwaite's method Sum Sq Mean Sq NumDF DenDF F value Pr(>F) a 336.00 336.00 1 6 3.9013 0.09566 . b 114.33 114.33 1 6 1.3275 0.29308 a:b 85.75 85.75 1 6 0.9956 0.35689 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > (an2 <- anova(model2, type=2)) Type II Analysis of Variance Table with Satterthwaite's method Sum Sq Mean Sq NumDF DenDF F value Pr(>F) b 114.33 114.33 1 6 1.3275 0.29308 a 336.00 336.00 1 6 3.9013 0.09566 . b:a 85.75 85.75 1 6 0.9956 0.35689 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > stopifnot( + isTRUE(all.equal(an,an2[c(2,1,3),], check.attributes = FALSE, tolerance=TOL2)) + ) > > ## check the results are the same as from SAS proc mixed > stopifnot( + isTRUE(all.equal(an[,"F value"], c(3.90131, 1.32753, 0.99565), tolerance=TOL2)) + ) > ################################################################################ > ## Check type II and III anova tables versus SAS > > m.carrots <- lmer(Preference ~ sens2*Homesize + +(1+sens2|Consumer), data=carrots) > (ancar <- anova(m.carrots, type=2)) Type II Analysis of Variance Table with Satterthwaite's method Sum Sq Mean Sq NumDF DenDF F value Pr(>F) sens2 58.697 58.697 1 101.02 54.8339 4.042e-11 *** Homesize 5.526 5.526 1 100.99 5.1621 0.02521 * sens2:Homesize 1.103 1.103 1 101.02 1.0303 0.31251 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > > stopifnot( + isTRUE(all.equal(ancar[,"F value"], c(54.8361, 5.16138, 1.03035), tolerance = TOL)) + ) > > m <- lmer(Informed.liking ~ Product*Age + + (1|Consumer) , data=ham) > (an <- anova(m, type=2)) Type II Analysis of Variance Table with Satterthwaite's method Sum Sq Mean Sq NumDF DenDF F value Pr(>F) Product 32.498 10.8327 3 561 2.4814 0.06014 . Age 0.024 0.0235 1 79 0.0054 0.94168 Product:Age 19.442 6.4808 3 561 1.4845 0.21772 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > > stopifnot( + isTRUE(all.equal(an[,"F value"], c(2.48135, .005387, 1.48451), tolerance = TOL2)) + ) > > > fm <- lmer(Preference ~ sens2*Homesize*sens1 + (1|Product), + data=carrots) > (ant2 <- anova(fm, type=2)) Type II Analysis of Variance Table with Satterthwaite's method Sum Sq Mean Sq NumDF DenDF F value Pr(>F) sens2 21.0228 21.0228 1 7.97 16.4832 0.0036583 ** Homesize 17.8570 17.8570 1 1217.01 14.0010 0.0001912 *** sens1 0.6709 0.6709 1 7.98 0.5260 0.4889911 sens2:Homesize 1.5068 1.5068 1 1216.99 1.1814 0.2772780 sens2:sens1 0.1372 0.1372 1 7.98 0.1076 0.7513796 Homesize:sens1 0.4275 0.4275 1 1217.00 0.3352 0.5627339 sens2:Homesize:sens1 1.3512 1.3512 1 1216.99 1.0595 0.3035426 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > (ant3 <- anova(fm, type=3)) Type III Analysis of Variance Table with Satterthwaite's method Sum Sq Mean Sq NumDF DenDF F value Pr(>F) sens2 21.5710 21.5710 1 8.02 16.9130 0.0033602 ** Homesize 17.8570 17.8570 1 1217.01 14.0010 0.0001912 *** sens1 0.6136 0.6136 1 8.03 0.4811 0.5074890 sens2:Homesize 1.5068 1.5068 1 1216.99 1.1814 0.2772780 sens2:sens1 0.0946 0.0946 1 8.02 0.0742 0.7922017 Homesize:sens1 0.4275 0.4275 1 1217.00 0.3352 0.5627339 sens2:Homesize:sens1 1.3512 1.3512 1 1216.99 1.0595 0.3035426 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > > stopifnot( + isTRUE(all.equal(ant2[,"F value"], + c(16.4842, 14.0010, .526076, 1.18144, + .107570, .335177, 1.05946), tolerance = TOL)), + isTRUE(all.equal(ant3[,"F value"], + c(16.9140, 14.0010,.481148, 1.18144, + .074201, .335177, 1.05946), tolerance = TOL)) + ) > > ################################################################################ > > > > proc.time() user system elapsed 2.73 0.26 2.93