R version 4.4.0 alpha (2024-03-26 r86209 ucrt) Copyright (C) 2024 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. > require(confintROB) Loading required package: confintROB > require(lme4) Loading required package: lme4 Loading required package: Matrix > require(robustlmm) Loading required package: robustlmm > > test <- function(object, ...) { + doTest <- function(method, ...) { + cat( + "Running test for object of class ", + class(object), + " with arguments method = ", + method, + "\n" + ) + set.seed(1234) + result <- + confintROB( + object = object, + level = .95, + method = method, + nsim = 10, + boot.type = "parametric", + clusterID = "Subject", + ... + ) + print(result, digits = 2) + } + + for (method in c("boot", "BCa", "Wald")) { + doTest(method, ...) + } + } > > test.wild <- function(object, ...) { + cat("Running test.wild for object of class ", class(object), "\n") + set.seed(123) + y <- getME(object, "y") + X <- as.matrix(getME(object, "X")) + id <- getME(object, "flist")[[1]] + bet <- unname(fixef(object)) + result <- + confintROB:::createWildSampleFunction(y = y, + X = X, + id = id, + bet = bet)(1) + print(result, digits = 5) + } > > > control <- lmerControl(check.conv.grad = "ignore") > > model.ds.ML <- + lmer(Yield ~ (1 | Batch), + Dyestuff, + REML = FALSE, + control = control) > print(summary(model.ds.ML), digits = 2) Linear mixed model fit by maximum likelihood ['lmerMod'] Formula: Yield ~ (1 | Batch) Data: Dyestuff Control: control AIC BIC logLik deviance df.resid 333.3 337.5 -163.7 327.3 27 Scaled residuals: Min 1Q Median 3Q Max -1.43 -0.80 0.15 0.77 1.80 Random effects: Groups Name Variance Std.Dev. Batch (Intercept) 1388 37 Residual 2451 50 Number of obs: 30, groups: Batch, 6 Fixed effects: Estimate Std. Error t value (Intercept) 1528 18 86 > test(model.ds.ML) Running test for object of class lmerMod with arguments method = boot Ignoring argument 'clusterID' as it's not needed for this combination of arguments 2.5 % 97.5 % (Intercept) 1503 1535 Sigma Batch (Intercept) 17 42 Sigma Residual 40 61 attr(,"fullResults") Full results of confintROB, a list with components: "Percentile", "bootstrap_estimates" Running test for object of class lmerMod with arguments method = BCa 2.5 % 97.5 % (Intercept) 1500 1534 Sigma Batch (Intercept) 32 44 Sigma Residual 41 62 attr(,"fullResults") Full results of confintROB, a list with components: "BCa", "Percentile", "bootstrap_estimates", "biasBCa", "acc" Running test for object of class lmerMod with arguments method = Wald 2.5 % 97.5 % (Intercept) 1493 1562 > test.wild(model.ds.ML, .export = "control") Running test.wild for object of class lmerMod [[1]] 1 2 3 4 5 6 7 8 9 10 11 1516.5 1582.6 1582.6 1532.2 1494.4 1548.1 1572.9 1465.6 1581.1 1473.9 1485.0 12 13 14 15 16 17 18 19 20 21 22 1513.3 1478.7 1538.5 1507.0 1391.4 1383.1 1638.9 1424.4 1556.4 1638.9 1696.6 23 24 25 26 27 28 29 30 1506.9 1704.9 1688.4 1532.2 1573.2 1576.3 1557.4 1579.5 > > model.ds.DAStau <- + rlmer( + Yield ~ (1 | Batch), + Dyestuff, + rho.sigma.e = psi2propII(smoothPsi, k = 2.28), + rho.b = chgDefaults(smoothPsi, k = 5.14, s = 10), + rho.sigma.b = chgDefaults(smoothPsi, k = 5.14, s = 10), + init = function(...) + lmer(..., control = control) + ) > print(summary(model.ds.DAStau), digits = 2) Robust linear mixed model fit by DAStau Formula: Yield ~ (1 | Batch) Data: Dyestuff Scaled residuals: Min 1Q Median 3Q Max -1.472 -0.682 0.084 0.738 1.934 Random effects: Groups Name Variance Std.Dev. Batch (Intercept) 2023 45 Residual 2517 50 Number of obs: 30, groups: Batch, 6 Fixed effects: Estimate Std. Error t value (Intercept) 1527 21 74 Robustness weights for the residuals: 25 weights are ~= 1. The remaining 5 ones are 2 3 5 18 23 0.950 0.950 0.930 0.695 0.901 Robustness weights for the random effects: [1] 1 1 1 1 1 1 Rho functions used for fitting: Residuals: eff: smoothed Huber (k = 1.345, s = 10) sig: smoothed Huber, Proposal 2 (k = 2.28, s = 10) Random Effects, variance component 1 (Batch): eff: smoothed Huber (k = 5.14, s = 10) vcp: smoothed Huber (k = 5.14, s = 10) > test(model.ds.DAStau) Running test for object of class rlmerMod with arguments method = boot 2.5 % 97.5 % (Intercept) 1502 1534 Sigma Batch (Intercept) 23 46 Sigma Residual 41 61 attr(,"fullResults") Full results of confintROB, a list with components: "Percentile", "bootstrap_estimates" Running test for object of class rlmerMod with arguments method = BCa 2.5 % 97.5 % (Intercept) 1499 1533 Sigma Batch (Intercept) 38 49 Sigma Residual 42 63 attr(,"fullResults") Full results of confintROB, a list with components: "BCa", "Percentile", "bootstrap_estimates", "biasBCa", "acc" Running test for object of class rlmerMod with arguments method = Wald 2.5 % 97.5 % (Intercept) 1486 1567 > test.wild(model.ds.DAStau, .export = "control") Running test.wild for object of class rlmerMod [[1]] 1 2 3 4 5 6 7 8 9 10 11 1515.4 1581.5 1581.5 1531.1 1493.3 1548.6 1573.3 1466.1 1581.6 1474.3 1483.8 12 13 14 15 16 17 18 19 20 21 22 1512.2 1477.5 1537.4 1505.9 1391.8 1383.6 1639.3 1424.8 1556.8 1639.3 1697.1 23 24 25 26 27 28 29 30 1507.3 1705.3 1688.8 1531.1 1572.1 1575.2 1556.3 1578.4 > > model.ss.ML <- + lmer( + Reaction ~ Days + (Days | Subject), + data = sleepstudy, + REML = FALSE, + control = control + ) > print(summary(model.ss.ML), digits = 2) Linear mixed model fit by maximum likelihood ['lmerMod'] Formula: Reaction ~ Days + (Days | Subject) Data: sleepstudy Control: control AIC BIC logLik deviance df.resid 1763.9 1783.1 -876.0 1751.9 174 Scaled residuals: Min 1Q Median 3Q Max -3.94 -0.47 0.03 0.46 5.18 Random effects: Groups Name Variance Std.Dev. Corr Subject (Intercept) 565 23.8 Days 33 5.7 0.08 Residual 655 25.6 Number of obs: 180, groups: Subject, 18 Fixed effects: Estimate Std. Error t value (Intercept) 251.4 6.6 38 Days 10.5 1.5 7 Correlation of Fixed Effects: (Intr) Days -0.138 > test.wild(model.ss.ML, .export = "control") Running test.wild for object of class lmerMod [[1]] 1 2 3 4 5 6 7 8 9 10 11 252.56 263.84 285.72 258.84 253.86 234.96 272.03 346.13 275.73 270.23 204.55 12 13 14 15 16 17 18 19 20 21 22 169.61 159.53 155.97 154.42 161.28 150.86 150.74 154.48 168.61 284.09 303.92 23 24 25 26 27 28 29 30 31 32 33 295.96 313.80 332.93 355.37 363.08 367.49 381.29 406.85 366.04 324.67 291.07 34 35 36 37 38 39 40 41 42 43 44 286.58 281.14 293.75 259.04 314.24 286.58 359.40 310.57 299.57 320.29 343.40 45 46 47 48 49 50 51 52 53 54 55 330.61 286.83 275.01 341.17 267.67 388.06 261.73 273.74 271.95 266.08 278.28 56 57 58 59 60 61 62 63 64 65 66 299.86 227.39 310.91 338.16 402.89 231.16 244.64 269.59 272.26 290.86 282.40 67 68 69 70 71 72 73 74 75 76 77 303.16 309.67 336.25 335.35 274.40 285.23 225.22 237.11 270.15 272.01 300.11 78 79 80 81 82 83 84 85 86 87 88 335.81 336.13 397.40 257.52 254.36 283.43 290.25 319.20 334.18 356.86 380.19 89 90 91 92 93 94 95 96 97 98 99 397.02 413.26 213.35 229.54 260.37 243.53 248.36 249.13 258.35 267.51 259.99 100 101 102 103 104 105 106 107 108 109 110 274.88 226.40 210.44 218.00 237.52 224.20 248.62 261.19 297.72 337.00 355.47 111 112 113 114 115 116 117 118 119 120 121 248.35 273.34 282.36 299.73 308.37 287.63 273.74 300.92 298.20 318.49 251.95 122 123 124 125 126 127 128 129 130 131 132 238.10 273.86 284.18 306.57 303.19 330.62 360.76 343.61 344.39 269.96 239.26 133 134 135 136 137 138 139 140 141 142 143 238.46 243.07 258.89 273.31 289.26 302.46 309.93 318.81 238.60 257.79 281.72 144 145 146 147 148 149 150 151 152 153 154 286.00 279.92 295.39 324.18 310.11 331.95 332.56 208.68 217.30 217.96 214.05 155 156 157 158 159 160 161 162 163 164 165 251.50 369.39 260.52 361.94 384.07 389.12 239.87 255.29 268.89 283.44 302.02 166 167 168 169 170 171 172 173 174 175 176 315.66 348.29 337.13 325.41 330.72 240.16 254.65 256.65 265.55 297.06 287.71 177 178 179 180 301.63 313.16 313.98 334.05 > > model.ss.DAStau <- + rlmer( + Reaction ~ Days + (Days | Subject), + data = sleepstudy, + rho.sigma.e = psi2propII(smoothPsi, k = 2.28), + rho.b = chgDefaults(smoothPsi, k = 5.14, s = 10), + rho.sigma.b = chgDefaults(smoothPsi, k = 5.14, s = 10), + init = function(...) + lmer(..., control = control) + ) > print(summary(model.ss.DAStau), digits = 2) Robust linear mixed model fit by DAStau Formula: Reaction ~ Days + (Days | Subject) Data: sleepstudy Scaled residuals: Min 1Q Median 3Q Max -5.83 -0.54 -0.02 0.54 6.62 Random effects: Groups Name Variance Std.Dev. Corr Subject (Intercept) 782 28.0 Days 42 6.5 -0.04 Residual 399 20.0 Number of obs: 180, groups: Subject, 18 Fixed effects: Estimate Std. Error t value (Intercept) 251.1 7.3 35 Days 10.7 1.6 7 Correlation of Fixed Effects: (Intr) Days -0.141 Robustness weights for the residuals: 154 weights are ~= 1. The remaining 26 ones are summarized as Min. 1st Qu. Median Mean 3rd Qu. Max. 0.203 0.619 0.730 0.697 0.896 0.999 Robustness weights for the random effects: All 36 weights are ~= 1. Rho functions used for fitting: Residuals: eff: smoothed Huber (k = 1.345, s = 10) sig: smoothed Huber, Proposal 2 (k = 2.28, s = 10) Random Effects, variance component 1 (Subject): eff: smoothed Huber (k = 5.14, s = 10) vcp: smoothed Huber (k = 5.14, s = 10) > test.wild(model.ss.DAStau, control = control, .export = "control") Running test.wild for object of class rlmerMod [[1]] 1 2 3 4 5 6 7 8 9 10 11 252.03 263.65 285.86 259.31 254.66 236.09 273.49 347.91 277.85 272.69 204.75 12 13 14 15 16 17 18 19 20 21 22 169.69 159.48 155.79 154.11 160.85 150.30 150.05 153.66 167.65 283.56 303.72 23 24 25 26 27 28 29 30 31 32 33 296.09 314.27 333.72 356.50 364.54 369.28 383.41 409.30 366.24 324.74 291.02 34 35 36 37 38 39 40 41 42 43 44 286.41 280.83 293.32 258.48 313.55 285.76 358.44 310.78 299.64 320.23 343.22 45 46 47 48 49 50 51 52 53 54 55 330.30 286.39 274.45 340.48 266.85 387.10 261.21 273.54 272.09 266.55 279.08 56 57 58 59 60 61 62 63 64 65 66 300.99 228.85 312.70 340.28 405.34 230.63 244.44 269.72 272.73 291.65 283.53 67 68 69 70 71 72 73 74 75 76 77 304.61 311.45 338.37 337.81 274.60 285.30 225.16 236.93 269.84 271.58 299.54 78 79 80 81 82 83 84 85 86 87 88 335.12 335.30 396.44 256.99 254.16 283.56 290.72 320.00 335.31 358.32 381.98 89 90 91 92 93 94 95 96 97 98 99 399.14 415.71 212.82 229.35 260.50 244.00 249.15 250.26 259.80 269.30 262.11 100 101 102 103 104 105 106 107 108 109 110 277.33 226.60 210.52 217.94 237.34 223.89 248.18 260.63 297.03 336.17 354.51 111 112 113 114 115 116 117 118 119 120 121 247.83 273.14 282.50 300.20 309.17 288.76 275.20 302.71 300.32 320.94 251.43 122 123 124 125 126 127 128 129 130 131 132 237.91 273.99 284.65 307.37 304.31 332.08 362.55 345.73 346.84 269.44 239.07 133 134 135 136 137 138 139 140 141 142 143 238.59 243.54 259.69 274.44 290.72 304.25 312.05 321.26 238.07 257.59 281.85 144 145 146 147 148 149 150 151 152 153 154 286.47 280.71 296.52 325.64 311.90 334.07 335.02 208.89 217.37 217.90 213.87 155 156 157 158 159 160 161 162 163 164 165 251.20 368.96 259.95 361.24 383.25 388.16 239.35 255.10 269.02 283.90 302.81 166 167 168 169 170 171 172 173 174 175 176 316.79 349.75 338.92 327.53 333.17 239.64 254.46 256.79 266.01 297.85 288.83 177 178 179 180 303.09 314.94 316.10 336.51 > > proc.time() user system elapsed 10.45 0.62 11.07