R Under development (unstable) (2024-06-18 r86781 ucrt) -- "Unsuffered Consequences" 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. > library("MuMIn") > options(na.action = "na.fail") > data(Orthodont, package = "nlme") > > fm1 <- lm(distance ~ Sex * age + age * Sex, data = Orthodont) > > dispersion <- function(object) { + wts <- weights(object) + if (is.null(wts)) + wts <- 1 + sum((wts * resid(object, type = "working")^2)[wts > 0])/df.residual(object) + } > > dd <- dredge(fm1, extra = alist(dispersion)) Fixed term is "(Intercept)" > gm <- get.models(dd, subset = 1:4) > ma <- model.avg(gm, revised = F) > > vcov(ma) (Intercept) SexFemale age SexFemale:age (Intercept) 2.1094155 -2.2423049 -0.18448478 0.17510000 SexFemale -2.2423049 5.5038392 0.19656574 -0.42979092 age -0.1844848 0.1965657 0.01677137 -0.01591818 SexFemale:age 0.1751000 -0.4297909 -0.01591818 0.03907190 > summary(ma) Call: model.avg(object = gm, revised = F) Component model call: lm(formula = distance ~ <4 unique rhs>, data = Orthodont) Component models: df logLik AICc delta weight 123 5 -239.12 488.83 0.00 0.53 12 4 -240.34 489.07 0.24 0.47 2 3 -252.79 511.81 22.98 0.00 1 3 -259.82 525.87 37.04 0.00 Term codes: Sex age Sex:age 1 2 3 Model-averaged coefficients: (full average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 16.9824 1.4524 1.4657 11.587 <2e-16 *** SexFemale -0.5432 2.3460 2.3596 0.230 0.818 age 0.7260 0.1295 0.1307 5.556 <2e-16 *** SexFemale:age -0.1616 0.2094 0.2106 0.767 0.443 (conditional average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 16.9824 1.4524 1.4657 11.587 <2e-16 *** SexFemale -0.5432 2.3460 2.3596 0.230 0.818 age 0.7260 0.1295 0.1307 5.556 <2e-16 *** SexFemale:age -0.3048 0.1977 0.2000 1.524 0.127 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > confint(ma) 2.5 % 97.5 % (Intercept) 14.1097091 19.85508711 SexFemale -5.1679459 4.08162225 age 0.4699087 0.98215438 SexFemale:age -0.6968089 0.08714982 > > predict(ma) 1 2 3 4 5 6 7 8 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 9 10 11 12 13 14 15 16 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 17 18 19 20 21 22 23 24 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 25 26 27 28 29 30 31 32 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 33 34 35 36 37 38 39 40 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 41 42 43 44 45 46 47 48 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 49 50 51 52 53 54 55 56 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 57 58 59 60 61 62 63 64 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 65 66 67 68 69 70 71 72 20.95451 22.08333 23.21214 24.34096 20.95451 22.08333 23.21214 24.34096 73 74 75 76 77 78 79 80 20.95451 22.08333 23.21214 24.34096 20.95451 22.08333 23.21214 24.34096 81 82 83 84 85 86 87 88 20.95451 22.08333 23.21214 24.34096 20.95451 22.08333 23.21214 24.34096 89 90 91 92 93 94 95 96 20.95451 22.08333 23.21214 24.34096 20.95451 22.08333 23.21214 24.34096 97 98 99 100 101 102 103 104 20.95451 22.08333 23.21214 24.34096 20.95451 22.08333 23.21214 24.34096 105 106 107 108 20.95451 22.08333 23.21214 24.34096 > predict(ma, se.fit = TRUE) $fit [1] 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 [9] 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 [17] 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 [25] 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 [33] 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 [41] 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 [49] 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 [57] 22.79065 24.24271 25.69478 27.14684 22.79065 24.24271 25.69478 27.14684 [65] 20.95451 22.08333 23.21214 24.34096 20.95451 22.08333 23.21214 24.34096 [73] 20.95451 22.08333 23.21214 24.34096 20.95451 22.08333 23.21214 24.34096 [81] 20.95451 22.08333 23.21214 24.34096 20.95451 22.08333 23.21214 24.34096 [89] 20.95451 22.08333 23.21214 24.34096 20.95451 22.08333 23.21214 24.34096 [97] 20.95451 22.08333 23.21214 24.34096 20.95451 22.08333 23.21214 24.34096 [105] 20.95451 22.08333 23.21214 24.34096 $se.fit [1] 0.4806525 0.3112188 0.3112210 0.4806567 0.4806525 0.3112188 0.3112210 [8] 0.4806567 0.4806525 0.3112188 0.3112210 0.4806567 0.4806525 0.3112188 [15] 0.3112210 0.4806567 0.4806525 0.3112188 0.3112210 0.4806567 0.4806525 [22] 0.3112188 0.3112210 0.4806567 0.4806525 0.3112188 0.3112210 0.4806567 [29] 0.4806525 0.3112188 0.3112210 0.4806567 0.4806525 0.3112188 0.3112210 [36] 0.4806567 0.4806525 0.3112188 0.3112210 0.4806567 0.4806525 0.3112188 [43] 0.3112210 0.4806567 0.4806525 0.3112188 0.3112210 0.4806567 0.4806525 [50] 0.3112188 0.3112210 0.4806567 0.4806525 0.3112188 0.3112210 0.4806567 [57] 0.4806525 0.3112188 0.3112210 0.4806567 0.4806525 0.3112188 0.3112210 [64] 0.4806567 0.5835523 0.3760119 0.3760157 0.5835597 0.5835523 0.3760119 [71] 0.3760157 0.5835597 0.5835523 0.3760119 0.3760157 0.5835597 0.5835523 [78] 0.3760119 0.3760157 0.5835597 0.5835523 0.3760119 0.3760157 0.5835597 [85] 0.5835523 0.3760119 0.3760157 0.5835597 0.5835523 0.3760119 0.3760157 [92] 0.5835597 0.5835523 0.3760119 0.3760157 0.5835597 0.5835523 0.3760119 [99] 0.3760157 0.5835597 0.5835523 0.3760119 0.3760157 0.5835597 0.5835523 [106] 0.3760119 0.3760157 0.5835597 > predict(ma, data.frame(Sex = "Male", age = 8:12)) 1 2 3 4 5 22.79065 23.51668 24.24271 24.96874 25.69478 > > rm(list = ls()) > > data(Cement, package = "MuMIn") > > nseq <- function(x, len = length(x)) seq(min(x, na.rm = TRUE), + max(x, na.rm = TRUE), length = len) > > fm1 <- glm(y ~ (X1 + X2 + X3)^2, data = Cement) > dd <- dredge(fm1) Fixed term is "(Intercept)" > > gm <- get.models(dd, subset = 1L:10L) > > summary(ma <- model.avg(gm)) Call: model.avg(object = gm) Component model call: glm(formula = y ~ <10 unique rhs>, data = Cement) Component models: df logLik AICc delta weight 12 4 -28.16 69.31 0.00 0.78 123 5 -26.95 72.48 3.16 0.16 124 5 -28.07 74.72 5.40 0.05 1236 6 -26.90 79.81 10.49 0.00 1234 6 -26.90 79.81 10.49 0.00 1235 6 -26.93 79.87 10.55 0.00 12346 7 -26.12 88.64 19.33 0.00 12356 7 -26.88 90.17 20.86 0.00 12345 7 -26.89 90.17 20.86 0.00 23 4 -40.96 94.93 25.62 0.00 Term codes: X1 X2 X3 X1:X2 X1:X3 X2:X3 1 2 3 4 5 6 Model-averaged coefficients: (full average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 5.189e+01 3.249e+00 3.602e+00 14.408 <2e-16 *** X1 1.497e+00 2.183e-01 2.446e-01 6.118 <2e-16 *** X2 6.597e-01 5.037e-02 5.749e-02 11.475 <2e-16 *** X3 4.256e-02 1.235e-01 1.322e-01 0.322 0.747 X1:X2 2.316e-04 3.065e-03 3.499e-03 0.066 0.947 X2:X3 1.025e-05 6.396e-04 7.438e-04 0.014 0.989 X1:X3 2.918e-05 3.050e-03 3.578e-03 0.008 0.993 (conditional average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 51.891051 3.248583 3.601545 14.408 <2e-16 *** X1 1.496664 0.218266 0.244626 6.118 <2e-16 *** X2 0.659732 0.050371 0.057491 11.475 <2e-16 *** X3 0.247629 0.195000 0.225527 1.098 0.272 X1:X2 0.004121 0.012294 0.014206 0.290 0.772 X2:X3 0.002465 0.009610 0.011271 0.219 0.827 X1:X3 0.007275 0.047602 0.056025 0.130 0.897 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > > vcov(ma) (Intercept) X1 X2 X3 (Intercept) 10.553293917 -4.078316e-01 -0.1042457812 -0.6488166976 X1 -0.407831558 4.763986e-02 0.0012054149 0.0314328032 X2 -0.104245781 1.205415e-03 0.0025372004 0.0004687388 X3 -0.648816698 3.143280e-02 0.0004687388 0.0380249494 X1:X2 0.041557193 -7.002098e-03 -0.0009596112 -0.0005732640 X2:X3 0.043537547 9.559883e-05 -0.0011790729 -0.0035872184 X1:X3 -0.006270502 -8.503388e-03 -0.0014709980 0.0005571689 X1:X2 X2:X3 X1:X3 (Intercept) 4.155719e-02 4.353755e-02 -6.270502e-03 X1 -7.002098e-03 9.559883e-05 -8.503388e-03 X2 -9.596112e-04 -1.179073e-03 -1.470998e-03 X3 -5.732640e-04 -3.587218e-03 5.571689e-04 X1:X2 1.511398e-04 7.343512e-04 -8.882992e-06 X2:X3 7.343512e-04 9.234858e-05 -7.080071e-06 X1:X3 -8.882992e-06 -7.080071e-06 2.265956e-03 > > summary(ma1 <- model.avg(dd[1L:10L])) Call: model.avg(object = dd[1:10]) Component model call: glm(formula = y ~ <10 unique rhs>, data = Cement) Component models: df logLik AICc delta weight 12 4 -28.16 69.31 0.00 0.78 123 5 -26.95 72.48 3.16 0.16 124 5 -28.07 74.72 5.40 0.05 1236 6 -26.90 79.81 10.49 0.00 1234 6 -26.90 79.81 10.49 0.00 1235 6 -26.93 79.87 10.55 0.00 12346 7 -26.12 88.64 19.33 0.00 12356 7 -26.88 90.17 20.86 0.00 12345 7 -26.89 90.17 20.86 0.00 23 4 -40.96 94.93 25.62 0.00 Term codes: X1 X2 X3 X1:X2 X1:X3 X2:X3 1 2 3 4 5 6 Model-averaged coefficients: (full average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 5.189e+01 3.249e+00 3.602e+00 14.408 <2e-16 *** X1 1.497e+00 2.183e-01 2.446e-01 6.118 <2e-16 *** X2 6.597e-01 5.037e-02 5.749e-02 11.475 <2e-16 *** X3 4.256e-02 1.235e-01 1.322e-01 0.322 0.747 X1:X2 2.316e-04 3.065e-03 3.499e-03 0.066 0.947 X2:X3 1.025e-05 6.396e-04 7.438e-04 0.014 0.989 X1:X3 2.918e-05 3.050e-03 3.578e-03 0.008 0.993 (conditional average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 51.891051 3.248583 3.601545 14.408 <2e-16 *** X1 1.496664 0.218266 0.244626 6.118 <2e-16 *** X2 0.659732 0.050371 0.057491 11.475 <2e-16 *** X3 0.247629 0.195000 0.225527 1.098 0.272 X1:X2 0.004121 0.012294 0.014206 0.290 0.772 X2:X3 0.002465 0.009610 0.011271 0.219 0.827 X1:X3 0.007275 0.047602 0.056025 0.130 0.897 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > summary(ma2 <- model.avg(model.sel(dd[1L:10L], rank = "AICc"))) New rank 'AICc' applied to logLik objects Call: model.avg(object = model.sel(dd[1:10], rank = "AICc")) Component model call: glm(formula = y ~ <10 unique rhs>, data = Cement) Component models: df logLik AICc delta weight 12 4 -28.16 69.31 0.00 0.78 123 5 -26.95 72.48 3.16 0.16 124 5 -28.07 74.72 5.40 0.05 1236 6 -26.90 79.81 10.49 0.00 1234 6 -26.90 79.81 10.49 0.00 1235 6 -26.93 79.87 10.55 0.00 12346 7 -26.12 88.64 19.33 0.00 12356 7 -26.88 90.17 20.86 0.00 12345 7 -26.89 90.17 20.86 0.00 23 4 -40.96 94.93 25.62 0.00 Term codes: X1 X2 X3 X1:X2 X1:X3 X2:X3 1 2 3 4 5 6 Model-averaged coefficients: (full average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 5.189e+01 3.249e+00 3.602e+00 14.408 <2e-16 *** X1 1.497e+00 2.183e-01 2.446e-01 6.118 <2e-16 *** X2 6.597e-01 5.037e-02 5.749e-02 11.475 <2e-16 *** X3 4.256e-02 1.235e-01 1.322e-01 0.322 0.747 X1:X2 2.316e-04 3.065e-03 3.499e-03 0.066 0.947 X2:X3 1.025e-05 6.396e-04 7.438e-04 0.014 0.989 X1:X3 2.918e-05 3.050e-03 3.578e-03 0.008 0.993 (conditional average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 51.891051 3.248583 3.601545 14.408 <2e-16 *** X1 1.496664 0.218266 0.244626 6.118 <2e-16 *** X2 0.659732 0.050371 0.057491 11.475 <2e-16 *** X3 0.247629 0.195000 0.225527 1.098 0.272 X1:X2 0.004121 0.012294 0.014206 0.290 0.772 X2:X3 0.002465 0.009610 0.011271 0.219 0.827 X1:X3 0.007275 0.047602 0.056025 0.130 0.897 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > all.equal(ma$avg.model, ma1$avg.mode) [1] TRUE > > predict(ma) == predict(ma, Cement) 1 2 3 4 5 6 7 8 9 10 11 12 13 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE > predict(ma, se.fit = TRUE) $fit [1] 79.82105 73.16996 105.78963 89.23059 97.01783 105.17072 104.00872 [8] 74.79054 91.31201 114.73157 80.77524 112.45681 112.22533 $se.fit [1] 1.3923364 1.2766548 0.8358514 1.2173773 1.0359279 0.8194428 1.4824953 [8] 1.3213611 1.0399669 1.8658856 1.2287624 1.0871067 1.1500574 > predict(ma, lapply(Cement, nseq)) 1 2 3 4 5 6 7 8 70.71819 75.76731 80.81961 85.87508 90.93371 95.99552 101.06050 106.12865 9 10 11 12 13 111.19997 116.27447 121.35213 126.43297 131.51697 > > proc.time() user system elapsed 1.62 0.34 1.93