Conversion between a model object and a restriction matrix
model-coerce.RdTesting a small model under a large model corresponds imposing restrictions on the model matrix of the larger model and these restrictions come in the form of a restriction matrix. These functions converts a model to a restriction matrix and vice versa.
model2restriction_matrix(largeModel, smallModel, sparse = FALSE)
restriction_matrix2model(largeModel, L, REML = TRUE, ...)
make_model_matrix(X, L)
make_restriction_matrix(X, X2)Arguments
- largeModel, smallModel
Model objects of the same "type". Possible types are linear mixed effects models and linear models (including generalized linear models)
- sparse
Should the restriction matrix be sparse or dense?
- L
A restriction matrix; a full rank matrix with as many columns as
Xhas.- REML
Controls if new model object should be fitted with REML or ML.
- ...
Additional arguments; not used.
- X, X2
Model matrices. Must have same number of rows.
Value
model2restriction_matrix: A restriction matrix.
restriction_matrix2model: A model object.
Details
make_restriction_matrix Make a restriction matrix. If span(X2) is in
span(X) then the corresponding restriction matrix L is
returned.
Note
That these functions are visible is a recent addition; minor changes may occur.
References
Ulrich Halekoh, Søren Højsgaard (2014)., A Kenward-Roger Approximation and Parametric Bootstrap Methods for Tests in Linear Mixed Models - The R Package pbkrtest., Journal of Statistical Software, 58(10), 1-30., https://www.jstatsoft.org/v59/i09/
Examples
library(pbkrtest)
data("beets", package = "pbkrtest")
sug <- lm(sugpct ~ block + sow + harvest, data=beets)
sug.h <- update(sug, .~. - harvest)
sug.s <- update(sug, .~. - sow)
## Construct restriction matrices from models
L.h <- model2restriction_matrix(sug, sug.h); L.h
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] 0 0 0 0 0 0 0 -1
L.s <- model2restriction_matrix(sug, sug.s); L.s
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] 0 0 0 -0.9176629 0.2294157 0.2294157 0.2294157 0
#> [2,] 0 0 0 0.0636285 0.8907987 -0.3181424 -0.3181424 0
#> [3,] 0 0 0 0.1048285 0.1048285 0.8386279 -0.5241424 0
#> [4,] 0 0 0 -0.3779645 -0.3779645 -0.3779645 -0.7559289 0
## Construct submodels from restriction matrices
mod.h <- restriction_matrix2model(sug, L.h); mod.h
#>
#> Call:
#> lm(formula = sugpct ~ .X1 + .X2 + .X3 + .X4 + .X5 + .X6 + .X7 -
#> 1, data = structure(list(.X1 = c(0, 0, 0, 0, 0, 1, 1, 1,
#> 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0,
#> 0), .X2 = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0,
#> 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1), .X3 = c(0, 0, 0, 1, 0,
#> 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
#> 0, 0, 1, 0), .X4 = c(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0,
#> 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0), .X5 = c(0, 1,
#> 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0,
#> 0, 0, 0, 1, 0, 0, 0), .X6 = c(0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1,
#> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1), .X7 = c(1,
#> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
#> 1, 1, 1, 1, 1, 1, 1, 1), sugpct = c(17.1, 16.9, 16.6, 17, 17,
#> 17, 17, 16.7, 16.4, 16.8, 16.6, 16.9, 17.1, 16.8, 16.9, 17, 16.8,
#> 16.5, 16.7, 17, 16.8, 16.9, 17, 17, 16.5, 16.7, 16.6, 16.9, 16.9,
#> 16.4), block = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
#> 2L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
#> 3L, 3L, 3L, 3L, 3L), levels = c("block1", "block2", "block3"), class = "factor"),
#> sow = structure(c(3L, 4L, 5L, 2L, 1L, 3L, 2L, 4L, 5L, 1L,
#> 5L, 2L, 3L, 4L, 1L, 2L, 1L, 5L, 4L, 3L, 4L, 1L, 3L, 2L, 5L,
#> 1L, 4L, 3L, 2L, 5L), levels = c("sow1", "sow2", "sow3", "sow4",
#> "sow5"), class = "factor"), harvest = structure(c(1L, 1L,
#> 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
#> 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L), levels = c("harv1",
#> "harv2"), class = "factor")), class = "data.frame", row.names = c(NA,
#> 30L)))
#>
#> Coefficients:
#> .X1 .X2 .X3 .X4 .X5 .X6 .X7
#> -0.0500 -0.0800 0.1167 0.1667 -0.1000 -0.3500 16.8933
#>
mod.s <- restriction_matrix2model(sug, L.s); mod.s
#>
#> Call:
#> lm(formula = sugpct ~ .X1 + .X2 + .X3 + .X4 - 1, data = structure(list(
#> .X1 = c(-0.481123223875107, -0.481123223875107, -0.481123223875107,
#> -0.481123223875107, -0.481123223875107, -1.35446911968096,
#> -1.35446911968096, -1.35446911968096, -1.35446911968096,
#> -1.35446911968096, -0.557198119386703, -0.557198119386703,
#> -0.557198119386703, -0.557198119386703, -0.557198119386703,
#> -0.481123223875107, -0.481123223875107, -0.481123223875107,
#> -0.481123223875107, -0.481123223875107, -1.35446911968096,
#> -1.35446911968096, -1.35446911968096, -1.35446911968096,
#> -1.35446911968096, -0.557198119386703, -0.557198119386703,
#> -0.557198119386703, -0.557198119386703, -0.557198119386703
#> ), .X2 = c(-0.481123205404825, -0.481123205404825, -0.481123205404825,
#> -0.481123205404825, -0.481123205404825, -0.145528052837381,
#> -0.145528052837381, -0.145528052837381, -0.145528052837381,
#> -0.145528052837381, -1.2909974887317, -1.2909974887317, -1.2909974887317,
#> -1.2909974887317, -1.2909974887317, -0.481123205404825, -0.481123205404825,
#> -0.481123205404825, -0.481123205404825, -0.481123205404825,
#> -0.145528052837381, -0.145528052837381, -0.145528052837381,
#> -0.145528052837381, -0.145528052837381, -1.2909974887317,
#> -1.2909974887317, -1.2909974887317, -1.2909974887317, -1.2909974887317
#> ), .X3 = c(-0.732830747627998, -0.732830747627998, -0.732830747627998,
#> -0.732830747627998, -0.732830747627998, -0.379782819079606,
#> -0.379782819079606, -0.379782819079606, -0.379782819079606,
#> -0.379782819079606, -0.151181149090632, -0.151181149090632,
#> -0.151181149090632, -0.151181149090632, -0.151181149090632,
#> -0.732830747627998, -0.732830747627998, -0.732830747627998,
#> -0.732830747627998, -0.732830747627998, -0.379782819079606,
#> -0.379782819079606, -0.379782819079606, -0.379782819079606,
#> -0.379782819079606, -0.151181149090632, -0.151181149090632,
#> -0.151181149090632, -0.151181149090632, -0.151181149090632
#> ), .X4 = c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1,
#> 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1), sugpct = c(17.1,
#> 16.9, 16.6, 17, 17, 17, 17, 16.7, 16.4, 16.8, 16.6, 16.9,
#> 17.1, 16.8, 16.9, 17, 16.8, 16.5, 16.7, 17, 16.8, 16.9, 17,
#> 17, 16.5, 16.7, 16.6, 16.9, 16.9, 16.4), block = structure(c(1L,
#> 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 1L,
#> 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L), levels = c("block1",
#> "block2", "block3"), class = "factor"), sow = structure(c(3L,
#> 4L, 5L, 2L, 1L, 3L, 2L, 4L, 5L, 1L, 5L, 2L, 3L, 4L, 1L, 2L,
#> 1L, 5L, 4L, 3L, 4L, 1L, 3L, 2L, 5L, 1L, 4L, 3L, 2L, 5L), levels = c("sow1",
#> "sow2", "sow3", "sow4", "sow5"), class = "factor"), harvest = structure(c(1L,
#> 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L,
#> 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L), levels = c("harv1",
#> "harv2"), class = "factor")), class = "data.frame", row.names = c(NA,
#> 30L)))
#>
#> Coefficients:
#> .X1 .X2 .X3 .X4
#> -8.0892 -8.0910 -12.4612 -0.1133
#>
## Sanity check: The models have the same fitted values and log likelihood
plot(fitted(mod.h), fitted(sug.h))
plot(fitted(mod.s), fitted(sug.s))
logLik(mod.h)
#> 'log Lik.' 36.03257 (df=8)
logLik(sug.h)
#> 'log Lik.' 36.03257 (df=8)
logLik(mod.s)
#> 'log Lik.' 7.397588 (df=5)
logLik(sug.s)
#> 'log Lik.' 7.397588 (df=5)