StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

PO
text
Body
copied!<p>I think the most straightforward solution is to use <code>model.matrix</code>. Possibly, you could achieve what you want with some fancy footwork and custom contrasts. However, if you want "type 3 esque" p-values, You probably want it for every term in your model, in which case, I think my approach with <code>model.matrix</code> is convenient anyway because you can easily implicitly loop through all models dropping one column at a time. The provision of a possible approach is not an endorsement of the statistical merits of it, but I do think you formulated a clear question and seem to know it may be unsound statistically so I see no reason not to answer it.</p> <pre><code>## initial data set.seed(10) d <- data.frame( A = rep(c("a1", "a2"), each = 50), B = c("b1", "b2"), value = rnorm(100)) options(contrasts=c('contr.sum','contr.poly')) ## create design matrix X <- model.matrix(~ A * B, data = d) ## fit models dropping one effect at a time ## change from 1:ncol(X) to 2:ncol(X) ## to avoid a no intercept model m <- lapply(1:ncol(X), function(i) { lm(value ~ 0 + X[, -i], data = d) }) ## fit (and store) the full model m$full <- lm(value ~ 0 + X, data = d) ## fit the full model in usual way to compare ## full and regular should be equivalent m$regular <- lm(value ~ A * B, data = d) ## extract and view coefficients lapply(m, coef) </code></pre> <p>This results in this final output:</p> <pre><code>[[1]] X[, -i]A1 X[, -i]B1 X[, -i]A1:B1 -0.2047465 -0.1330705 0.1133502 [[2]] X[, -i](Intercept) X[, -i]B1 X[, -i]A1:B1 -0.1365489 -0.1330705 0.1133502 [[3]] X[, -i](Intercept) X[, -i]A1 X[, -i]A1:B1 -0.1365489 -0.2047465 0.1133502 [[4]] X[, -i](Intercept) X[, -i]A1 X[, -i]B1 -0.1365489 -0.2047465 -0.1330705 $full X(Intercept) XA1 XB1 XA1:B1 -0.1365489 -0.2047465 -0.1330705 0.1133502 $regular (Intercept) A1 B1 A1:B1 -0.1365489 -0.2047465 -0.1330705 0.1133502 </code></pre> <p>That is nice so far for models using <code>lm</code>. You mentioned this is ultimately for <code>lmer()</code>, so here is an example using mixed models. I believe it may become more complex if you have more than a random intercept (i.e., effects need to be dropped from the fixed and random parts of the model).</p> <pre><code>## mixed example require(lme4) ## data is a bit trickier set.seed(10) mixed <- data.frame( ID = factor(ID <- rep(seq_along(n <- sample(3:8, 60, TRUE)), n)), A = sample(c("a1", "a2"), length(ID), TRUE), B = sample(c("b1", "b2"), length(ID), TRUE), value = rnorm(length(ID), 3) + rep(rnorm(length(n)), n)) ## model matrix as before X <- model.matrix(~ A * B, data = mixed) ## as before but allowing a random intercept by ID ## becomes trickier if you need to add/drop random effects too ## and I do not show an example of this mm <- lapply(1:ncol(X), function(i) { lmer(value ~ 0 + X[, -i] + (1 | ID), data = mixed) }) ## full model mm$full <- lmer(value ~ 0 + X + (1 | ID), data = mixed) ## full model regular way mm$regular <- lmer(value ~ A * B + (1 | ID), data = mixed) ## view all the fixed effects lapply(mm, fixef) </code></pre> <p>Which gives us...</p> <pre><code>[[1]] X[, -i]A1 X[, -i]B1 X[, -i]A1:B1 0.009202554 0.028834041 0.054651770 [[2]] X[, -i](Intercept) X[, -i]B1 X[, -i]A1:B1 2.83379928 0.03007969 0.05992235 [[3]] X[, -i](Intercept) X[, -i]A1 X[, -i]A1:B1 2.83317191 0.02058800 0.05862495 [[4]] X[, -i](Intercept) X[, -i]A1 X[, -i]B1 2.83680235 0.01738798 0.02482256 $full X(Intercept) XA1 XB1 XA1:B1 2.83440919 0.01947658 0.02928676 0.06057778 $regular (Intercept) A1 B1 A1:B1 2.83440919 0.01947658 0.02928676 0.06057778 </code></pre>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload