Construct the Design Matrix of a VLM Object

Creates a design matrix. Two types can be returned: a large one (class "vlm" or one that inherits from this such as "vglm") or a small one (such as returned if it were of class "lm").

Usage

model.matrixvlm(object, type = c("vlm", "lm", "lm2", "bothlmlm2"),
                linpred.index = NULL, label.it = TRUE, ...)

Arguments

object

an object of a class that inherits from the vector linear model (VLM).

type

Type of design matrix returned. The first is the default. The value "vlm" is the VLM model matrix corresponding to the formula argument. The value "lm" is the LM model matrix corresponding to the formula argument. The value "lm2" is the second (LM) model matrix corresponding to the form2 argument. The value "bothlmlm2" means both LM and VLM model matrices.

linpred.index

Vector of integers. The index for a linear/additive predictor, it must have values from the set 1:M. Also, if length(linpred.index) == 1 then type = "lm" must be assigned, whereas if length(linpred.index) > 1 then type = "vlm" must be assigned, Then it returns a subset of the VLM matrix corresponding to the linpred.indexth linear/additive predictor(s); this is a LM-type matrix when it is of unit length. Currently some attributes are returned, but these may change in value in the future because of ongoing development work.

label.it

Logical. Label the row and columns with character names? If FALSE, time and memory might be saved if the big model matrix is very large. The argument is only used when type = "vlm".

...

further arguments passed to or from other methods. These include data (which is a data frame created with model.framevlm), contrasts.arg, and xlev. See model.matrix for more information.

Details

This function creates a design matrix from object. This can be a small LM object or a big VLM object (default). The latter is constructed from the former and the constraint matrices.

This code implements smart prediction (see smartpred).

Value

The design matrix for a regression model with the specified formula and data. If type = "bothlmlm2" then a list is returned with components "X" and "Xm2".

Sometimes (especially if x = TRUE when calling vglm) the model matrix has attributes: "assign" ("lm"-type) and "vassign" ("vlm"-type) and "orig.assign.lm" ("lm"-type). These are used internally a lot for bookkeeping, especially regarding the columns of both types of model matrices. In particular, constraint matrices and variable selection relies on this information a lot. The "orig.assign.lm" is the ordinary "assign" attribute for lm and glm objects.

References

Chambers, J. M. (1992). Data for models. Chapter 3 of Statistical Models in S eds J. M. Chambers and T. J. Hastie, Wadsworth & Brooks/Cole.

Examples

# (I) Illustrates smart prediction ,,,,,,,,,,,,,,,,,,,,,,,
pneumo <- transform(pneumo, let = log(exposure.time))
fit <- vglm(cbind(normal, mild, severe) ~
            sm.poly(c(sm.scale(let)), 2),
            multinomial, data = pneumo, trace = TRUE, x = FALSE)
#> Iteration 1: deviance = 3.716285
#> Iteration 2: deviance = 3.499289
#> Iteration 3: deviance = 3.487318
#> Iteration 4: deviance = 3.487262
#> Iteration 5: deviance = 3.487262
class(fit)
#> [1] "vglm"
#> attr(,"package")
#> [1] "VGAM"
fit@smart.prediction  # Data-dependent parameters
#> [[1]]
#> [[1]]$center
#> [1] 3.225777
#> 
#> [[1]]$scale
#> [1] 0.7187268
#> 
#> [[1]]$match.call
#> sm.scale.default(x = let)
#> 
#> 
#> [[2]]
#> [[2]]$degree
#> [1] 2
#> 
#> [[2]]$coefs
#> [[2]]$coefs$alpha
#> [1]  2.706169e-16 -1.002119e+00
#> 
#> [[2]]$coefs$norm2
#> [1] 1.000000 8.000000 7.000000 6.175489
#> 
#> 
#> [[2]]$raw
#> [1] FALSE
#> 
#> [[2]]$match.call
#> sm.poly(x = c(sm.scale(let)), 2)
#> 
#> 
fit@x # Not saved on the object
#> <0 x 0 matrix>
model.matrix(fit)
#>     (Intercept):1 (Intercept):2 sm.poly(c(sm.scale(let)), 2)1:1
#> 1:1             1             0                     -0.77195016
#> 1:2             0             1                      0.00000000
#> 2:1             1             0                     -0.27226249
#> 2:2             0             1                      0.00000000
#> 3:1             1             0                     -0.08294405
#> 3:2             0             1                      0.00000000
#> 4:1             1             0                      0.04649255
#> 4:2             0             1                      0.00000000
#> 5:1             1             0                      0.15028005
#> 5:2             0             1                      0.00000000
#> 6:1             1             0                      0.23692162
#> 6:2             0             1                      0.00000000
#> 7:1             1             0                      0.31703465
#> 7:2             0             1                      0.00000000
#> 8:1             1             0                      0.37642782
#> 8:2             0             1                      0.00000000
#>     sm.poly(c(sm.scale(let)), 2)1:2 sm.poly(c(sm.scale(let)), 2)2:1
#> 1:1                      0.00000000                      0.50285988
#> 1:2                     -0.77195016                      0.00000000
#> 2:1                      0.00000000                     -0.43378435
#> 2:2                     -0.27226249                      0.00000000
#> 3:1                      0.00000000                     -0.42122091
#> 3:2                     -0.08294405                      0.00000000
#> 4:1                      0.00000000                     -0.29641247
#> 4:2                      0.04649255                      0.00000000
#> 5:1                      0.00000000                     -0.12815229
#> 5:2                      0.15028005                      0.00000000
#> 6:1                      0.00000000                      0.05878623
#> 6:2                      0.23692162                      0.00000000
#> 7:1                      0.00000000                      0.26926949
#> 7:2                      0.31703465                      0.00000000
#> 8:1                      0.00000000                      0.44865441
#> 8:2                      0.37642782                      0.00000000
#>     sm.poly(c(sm.scale(let)), 2)2:2
#> 1:1                      0.00000000
#> 1:2                      0.50285988
#> 2:1                      0.00000000
#> 2:2                     -0.43378435
#> 3:1                      0.00000000
#> 3:2                     -0.42122091
#> 4:1                      0.00000000
#> 4:2                     -0.29641247
#> 5:1                      0.00000000
#> 5:2                     -0.12815229
#> 6:1                      0.00000000
#> 6:2                      0.05878623
#> 7:1                      0.00000000
#> 7:2                      0.26926949
#> 8:1                      0.00000000
#> 8:2                      0.44865441
#> attr(,"assign")
#> attr(,"assign")$`(Intercept)`
#> [1] 1 2
#> 
#> attr(,"assign")$`sm.poly(c(sm.scale(let)), 2)`
#> [1] 3 4 5 6
#> 
#> attr(,"vassign")
#> attr(,"vassign")$`(Intercept):1`
#> [1] 1
#> 
#> attr(,"vassign")$`(Intercept):2`
#> [1] 2
#> 
#> attr(,"vassign")$`sm.poly(c(sm.scale(let)), 2):1`
#> [1] 3 5
#> 
#> attr(,"vassign")$`sm.poly(c(sm.scale(let)), 2):2`
#> [1] 4 6
#> 
#> attr(,"constraints")
#> attr(,"constraints")$`(Intercept)`
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> attr(,"constraints")$`sm.poly(c(sm.scale(let)), 2)1`
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> attr(,"constraints")$`sm.poly(c(sm.scale(let)), 2)2`
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
model.matrix(fit, linpred.index = 1, type = "lm")
#>     (Intercept):1 sm.poly(c(sm.scale(let)), 2)1:1
#> 1:1             1                     -0.77195016
#> 2:1             1                     -0.27226249
#> 3:1             1                     -0.08294405
#> 4:1             1                      0.04649255
#> 5:1             1                      0.15028005
#> 6:1             1                      0.23692162
#> 7:1             1                      0.31703465
#> 8:1             1                      0.37642782
#>     sm.poly(c(sm.scale(let)), 2)2:1
#> 1:1                      0.50285988
#> 2:1                     -0.43378435
#> 3:1                     -0.42122091
#> 4:1                     -0.29641247
#> 5:1                     -0.12815229
#> 6:1                      0.05878623
#> 7:1                      0.26926949
#> 8:1                      0.44865441
#> attr(,"vassign")
#> attr(,"vassign")$`(Intercept):1`
#> [1] 1
#> 
#> attr(,"vassign")$`sm.poly(c(sm.scale(let)), 2):1`
#> [1] 2 3
#> 
#> attr(,"rm.vassign")
#> attr(,"rm.vassign")$`(Intercept):1`
#> [1] 1
#> 
#> attr(,"rm.vassign")$`sm.poly(c(sm.scale(let)), 2):1`
#> [1] 3 5
#> 
#> attr(,"assign")
#> attr(,"assign")$`(Intercept)`
#> [1] 1
#> 
#> attr(,"assign")$`sm.poly(c(sm.scale(let)), 2)`
#> [1] 2 3
#> 
#> attr(,"rm.assign")
#> attr(,"rm.assign")$`(Intercept)`
#> [1] 1 2
#> 
#> attr(,"rm.assign")$`sm.poly(c(sm.scale(let)), 2)`
#> [1] 3 4 5 6
#> 
model.matrix(fit, linpred.index = 2, type = "lm")
#>     (Intercept):2 sm.poly(c(sm.scale(let)), 2)1:2
#> 1:2             1                     -0.77195016
#> 2:2             1                     -0.27226249
#> 3:2             1                     -0.08294405
#> 4:2             1                      0.04649255
#> 5:2             1                      0.15028005
#> 6:2             1                      0.23692162
#> 7:2             1                      0.31703465
#> 8:2             1                      0.37642782
#>     sm.poly(c(sm.scale(let)), 2)2:2
#> 1:2                      0.50285988
#> 2:2                     -0.43378435
#> 3:2                     -0.42122091
#> 4:2                     -0.29641247
#> 5:2                     -0.12815229
#> 6:2                      0.05878623
#> 7:2                      0.26926949
#> 8:2                      0.44865441
#> attr(,"vassign")
#> attr(,"vassign")$`(Intercept):2`
#> [1] 1
#> 
#> attr(,"vassign")$`sm.poly(c(sm.scale(let)), 2):2`
#> [1] 2 3
#> 
#> attr(,"rm.vassign")
#> attr(,"rm.vassign")$`(Intercept):2`
#> [1] 2
#> 
#> attr(,"rm.vassign")$`sm.poly(c(sm.scale(let)), 2):2`
#> [1] 4 6
#> 
#> attr(,"assign")
#> attr(,"assign")$`(Intercept)`
#> [1] 1
#> 
#> attr(,"assign")$`sm.poly(c(sm.scale(let)), 2)`
#> [1] 2 3
#> 
#> attr(,"rm.assign")
#> attr(,"rm.assign")$`(Intercept)`
#> [1] 1 2
#> 
#> attr(,"rm.assign")$`sm.poly(c(sm.scale(let)), 2)`
#> [1] 3 4 5 6
#> 

(Check1 <- head(model.matrix(fit, type = "lm")))
#>   (Intercept) sm.poly(c(sm.scale(let)), 2)1 sm.poly(c(sm.scale(let)), 2)2
#> 1           1                   -0.77195016                    0.50285988
#> 2           1                   -0.27226249                   -0.43378435
#> 3           1                   -0.08294405                   -0.42122091
#> 4           1                    0.04649255                   -0.29641247
#> 5           1                    0.15028005                   -0.12815229
#> 6           1                    0.23692162                    0.05878623
(Check2 <- model.matrix(fit, data = head(pneumo), type = "lm"))
#>   (Intercept) sm.poly(c(sm.scale(let)), 2)1 sm.poly(c(sm.scale(let)), 2)2
#> 1           1                   -0.77195016                    0.50285988
#> 2           1                   -0.27226249                   -0.43378435
#> 3           1                   -0.08294405                   -0.42122091
#> 4           1                    0.04649255                   -0.29641247
#> 5           1                    0.15028005                   -0.12815229
#> 6           1                    0.23692162                    0.05878623
#> attr(,"assign")
#> attr(,"assign")$`(Intercept)`
#> [1] 1
#> 
#> attr(,"assign")$`sm.poly(c(sm.scale(let)), 2)`
#> [1] 2 3
#> 
all.equal(c(Check1), c(Check2))  # Should be TRUE
#> [1] TRUE

q0 <- head(predict(fit))
q1 <- head(predict(fit, newdata = pneumo))
q2 <- predict(fit, newdata = head(pneumo))
all.equal(q0, q1)  # Should be TRUE
#> [1] TRUE
all.equal(q1, q2)  # Should be TRUE
#> [1] TRUE

# (II) Attributes ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
fit2 <- vglm(cbind(normal, mild, severe) ~ let,  # x = TRUE
             multinomial, data = pneumo, trace = TRUE)
#> Iteration 1: deviance = 5.407271
#> Iteration 2: deviance = 5.34745
#> Iteration 3: deviance = 5.347382
#> Iteration 4: deviance = 5.347382
fit2@x  # "lm"-type; saved on the object; note the attributes
#>   (Intercept)      let
#> 1           1 1.757858
#> 2           1 2.708050
#> 3           1 3.068053
#> 4           1 3.314186
#> 5           1 3.511545
#> 6           1 3.676301
#> 7           1 3.828641
#> 8           1 3.941582
#> attr(,"assign")
#> attr(,"assign")$`(Intercept)`
#> [1] 1
#> 
#> attr(,"assign")$let
#> [1] 2
#> 
#> attr(,"orig.assign.lm")
#> [1] 0 1
model.matrix(fit2, type = "lm")  # Note the attributes
#>   (Intercept)      let
#> 1           1 1.757858
#> 2           1 2.708050
#> 3           1 3.068053
#> 4           1 3.314186
#> 5           1 3.511545
#> 6           1 3.676301
#> 7           1 3.828641
#> 8           1 3.941582
#> attr(,"assign")
#> attr(,"assign")$`(Intercept)`
#> [1] 1
#> 
#> attr(,"assign")$let
#> [1] 2
#> 
#> attr(,"orig.assign.lm")
#> [1] 0 1
model.matrix(fit2, type = "vlm")  # Note the attributes
#>     (Intercept):1 (Intercept):2    let:1    let:2
#> 1:1             1             0 1.757858 0.000000
#> 1:2             0             1 0.000000 1.757858
#> 2:1             1             0 2.708050 0.000000
#> 2:2             0             1 0.000000 2.708050
#> 3:1             1             0 3.068053 0.000000
#> 3:2             0             1 0.000000 3.068053
#> 4:1             1             0 3.314186 0.000000
#> 4:2             0             1 0.000000 3.314186
#> 5:1             1             0 3.511545 0.000000
#> 5:2             0             1 0.000000 3.511545
#> 6:1             1             0 3.676301 0.000000
#> 6:2             0             1 0.000000 3.676301
#> 7:1             1             0 3.828641 0.000000
#> 7:2             0             1 0.000000 3.828641
#> 8:1             1             0 3.941582 0.000000
#> 8:2             0             1 0.000000 3.941582
#> attr(,"assign")
#> attr(,"assign")$`(Intercept)`
#> [1] 1 2
#> 
#> attr(,"assign")$let
#> [1] 3 4
#> 
#> attr(,"vassign")
#> attr(,"vassign")$`(Intercept):1`
#> [1] 1
#> 
#> attr(,"vassign")$`(Intercept):2`
#> [1] 2
#> 
#> attr(,"vassign")$`let:1`
#> [1] 3
#> 
#> attr(,"vassign")$`let:2`
#> [1] 4
#> 
#> attr(,"constraints")
#> attr(,"constraints")$`(Intercept)`
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> attr(,"constraints")$let
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> attr(,"orig.assign.lm")
#> [1] 0 1