Extract design matrix and data specials from a model.frame
Usage
model.design(
terms,
data,
xlev = NULL,
dropIntercept = FALSE,
maxOrder = 1,
unspecialsDesign = TRUE,
specialsFactor = FALSE,
specialsDesign = FALSE
)Arguments
- terms
terms object as obtained either with function
termsorstrip.terms.- data
A data set in which terms are defined.
- xlev
a named list of character vectors giving the full set of levels to be assumed for the factors. Can have less elements, in which case the other levels are learned from the
data.- dropIntercept
If TRUE drop intercept term from the design matrix
- maxOrder
An error is produced if special variables are involved in interaction terms of order higher than max.order.
- unspecialsDesign
A logical value: if
TRUEapplymodel.matrixto unspecial covariates. IfFALSEextract unspecial covariates from data.- specialsFactor
A character vector containing special variables which should be coerced into a single factor. If
TRUEall specials are treated in this way, ifFALSEnone of the specials is treated in this way.- specialsDesign
A character vector containing special variables which should be transformed into a design matrix via
model.matrix. IfTRUEall specials are treated in this way.
Value
A list which contains - the design matrix with the levels of the variables stored in attribute 'levels' - separate data.frames which contain the values of the special variables.
Details
The function separates special terms from the unspecial terms and returns
a list of design matrices, one for unspecial terms and one for each special.
Some special specials cannot or should not be evaluated in
data. E.g., y~a+dummy(x)+strata(v) the function strata can and should be evaluated,
but in order to have model.frame also evaluate dummy(x) one would be to define
and export the function dummy. Still the term dummy(x) can be used
to identify a special treatment of the variable x. To deal with this case,
one can specify stripSpecials="dummy". In addition, the data
should include variables strata(z) and x, not dummy(x).
See examples.
The function untangle.specials of the survival function does a similar job.
See also
EventHistory.frame model.frame terms model.matrix .getXlevels
Examples
# specials that are evaluated. here ID needs to be defined
set.seed(8)
d <- data.frame(y=rnorm(5),x=factor(c("a","b","b","a","c")),z=c(2,2,7,7,7),v=sample(letters)[1:5])
d$z <- factor(d$z,levels=c(1:8))
ID <- function(x)x
f <- formula(y~x+ID(z))
t <- terms(f,special="ID",data=d)
mda <- model.design(terms(t),data=d,specialsFactor=TRUE)
mda$ID
#> ID(z)
#> 1 2
#> 2 2
#> 3 7
#> 4 7
#> 5 7
mda$design
#> (Intercept) xb xc
#> 1 1 0 0
#> 2 1 1 0
#> 3 1 1 0
#> 4 1 0 0
#> 5 1 0 1
#> attr(,"assign")
#> [1] 0 1 1
#> attr(,"contrasts")
#> attr(,"contrasts")$x
#> [1] "contr.treatment"
#>
#> attr(,"levels")
#> attr(,"levels")$x
#> [1] "a" "b" "c"
#>
##
mdb <- model.design(terms(t),data=d,specialsFactor=TRUE,unspecialsDesign=FALSE)
mdb$ID
#> ID(z)
#> 1 2
#> 2 2
#> 3 7
#> 4 7
#> 5 7
mdb$design
#> x
#> 1 a
#> 2 b
#> 3 b
#> 4 a
#> 5 c
# set x-levels
attr(mdb$ID,"levels")
#> $`ID(z)`
#> [1] "1" "2" "3" "4" "5" "6" "7" "8"
#>
attr(model.design(terms(t),data=d,xlev=list("ID(z)"=1:10),
specialsFactor=TRUE)$ID,"levels")
#> $`ID(z)`
#> [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10"
#>
# special specials (avoid define function SP)
f <- formula(y~x+SP(z)+factor(v))
t <- terms(f,specials="SP",data=d)
st <- strip.terms(t,specials="SP",arguments=NULL)
md2a <- model.design(st,data=d,specialsFactor=TRUE,specialsDesign="SP")
md2a$SP
#> z2 z3 z4 z5 z6 z7 z8
#> 1 1 0 0 0 0 0 0
#> 2 1 0 0 0 0 0 0
#> 3 0 0 0 0 0 1 0
#> 4 0 0 0 0 0 1 0
#> 5 0 0 0 0 0 1 0
#> attr(,"levels")
#> attr(,"levels")$z
#> [1] "1" "2" "3" "4" "5" "6" "7" "8"
#>
#> attr(,"arguments")
#> attr(,"arguments")$z
#> NULL
#>
#> attr(,"arguments.terms")
#> list()
#> attr(,"matrix.terms")
#> attr(,"matrix.terms")$z
#> [1] "z1" "z2" "z3" "z4" "z5" "z6" "z7" "z8"
#>
md2b <- model.design(st,data=d,specialsFactor=TRUE,specialsDesign=FALSE)
md2b$SP
#> z
#> 1 2
#> 2 2
#> 3 7
#> 4 7
#> 5 7
# special function with argument
f2 <- formula(y~x+treat(z,power=2)+treat(v,power=-1))
t2 <- terms(f2,special="treat")
st2 <- strip.terms(t2,specials="treat",arguments=list("treat"=list("power")))
model.design(st2,data=d,specialsFactor=FALSE)
#> $design
#> (Intercept) xb xc
#> 1 1 0 0
#> 2 1 1 0
#> 3 1 1 0
#> 4 1 0 0
#> 5 1 0 1
#> attr(,"assign")
#> [1] 0 1 1
#> attr(,"contrasts")
#> attr(,"contrasts")$x
#> [1] "contr.treatment"
#>
#> attr(,"levels")
#> attr(,"levels")$x
#> [1] "a" "b" "c"
#>
#>
#> $treat
#> z v
#> 1 2 a
#> 2 2 r
#> 3 7 c
#> 4 7 x
#> 5 7 f
#>
model.design(st2,data=d,specialsFactor=TRUE)
#> $design
#> (Intercept) xb xc
#> 1 1 0 0
#> 2 1 1 0
#> 3 1 1 0
#> 4 1 0 0
#> 5 1 0 1
#> attr(,"assign")
#> [1] 0 1 1
#> attr(,"contrasts")
#> attr(,"contrasts")$x
#> [1] "contr.treatment"
#>
#> attr(,"levels")
#> attr(,"levels")$x
#> [1] "a" "b" "c"
#>
#>
#> $treat
#> z, v
#> 1 2, a
#> 2 2, r
#> 3 7, c
#> 4 7, x
#> 5 7, f
#>
model.design(st2,data=d,specialsDesign=TRUE)
#> $design
#> (Intercept) xb xc
#> 1 1 0 0
#> 2 1 1 0
#> 3 1 1 0
#> 4 1 0 0
#> 5 1 0 1
#> attr(,"assign")
#> [1] 0 1 1
#> attr(,"contrasts")
#> attr(,"contrasts")$x
#> [1] "contr.treatment"
#>
#> attr(,"levels")
#> attr(,"levels")$x
#> [1] "a" "b" "c"
#>
#>
#> $treat
#> z2 z3 z4 z5 z6 z7 z8 vc vf vr vx
#> 1 1 0 0 0 0 0 0 0 0 0 0
#> 2 1 0 0 0 0 0 0 0 0 1 0
#> 3 0 0 0 0 0 1 0 1 0 0 0
#> 4 0 0 0 0 0 1 0 0 0 0 1
#> 5 0 0 0 0 0 1 0 0 1 0 0
#> attr(,"levels")
#> attr(,"levels")$z
#> [1] "1" "2" "3" "4" "5" "6" "7" "8"
#>
#> attr(,"levels")$v
#> [1] "a" "c" "f" "r" "x"
#>
#> attr(,"arguments")
#> attr(,"arguments")$z
#> NULL
#>
#> attr(,"arguments")$v
#> NULL
#>
#> attr(,"arguments.terms")
#> list()
#> attr(,"matrix.terms")
#> attr(,"matrix.terms")$z
#> [1] "z1" "z2" "z3" "z4" "z5" "z6" "z7" "z8"
#>
#> attr(,"matrix.terms")$v
#> [1] "va" "vc" "vf" "vr" "vx"
#>
#>
library(survival)
data(pbc)
t3 <- terms(Surv(time,status!=0)~factor(edema)*age+strata(I(log(bili)>1))+strata(sex),
specials=c("strata","cluster"))
st3 <- strip.terms(t3,specials=c("strata"),arguments=NULL)
md3 <- model.design(terms=st3,data=pbc[1:4,])
md3$strata
#> I(log(bili) > 1) sex
#> 1 TRUE f
#> 2 FALSE f
#> 3 FALSE m
#> 4 FALSE f
md3$cluster
#> NULL
f4 <- Surv(time,status)~const(age)+prop(factor(edema),power=2)+
strata(sex,test=0)+prop(bili,power=1)+tp(albumin)
t4 <- terms(f4,specials=c("prop","timevar","strata","tp","const"))
st4 <- strip.terms(t4,
specials=c("prop","timevar"),
unspecials="prop",
alias.names=list("timevar"="strata","prop"=c("const","tp")),
arguments=list("prop"=list("power"=0),"timevar"=list("test"=0)))
formula(st4)
#> Surv(time, status) ~ age + factor(edema) + sex + bili + albumin
#> <environment: 0x6402e034b638>
md4 <- model.design(st4,data=pbc[1:4,],specialsDesign=TRUE)
md4$prop
#> age factor(edema)0.5 factor(edema)1 bili albumin
#> 1 58.76523 0 1 14.5 2.60
#> 2 56.44627 0 0 1.1 4.14
#> 3 70.07255 1 0 1.4 3.48
#> 4 54.74059 1 0 1.8 2.54
#> attr(,"levels")
#> attr(,"levels")$`factor(edema)`
#> [1] "0" "0.5" "1"
#>
#> attr(,"arguments")
#> attr(,"arguments")$age
#> attr(,"arguments")$age$power
#> [1] 0
#>
#>
#> attr(,"arguments")$`factor(edema)`
#> attr(,"arguments")$`factor(edema)`$power
#> [1] "2"
#>
#>
#> attr(,"arguments")$bili
#> attr(,"arguments")$bili$power
#> [1] "1"
#>
#>
#> attr(,"arguments")$albumin
#> attr(,"arguments")$albumin$power
#> [1] 0
#>
#>
#> attr(,"arguments.terms")
#> attr(,"arguments.terms")$power
#> age factor(edema)0 factor(edema)0.5 factor(edema)1
#> "0" "2" "2" "2"
#> bili albumin
#> "1" "0"
#>
#> attr(,"matrix.terms")
#> attr(,"matrix.terms")$age
#> [1] "age"
#>
#> attr(,"matrix.terms")$`factor(edema)`
#> [1] "factor(edema)0" "factor(edema)0.5" "factor(edema)1"
#>
#> attr(,"matrix.terms")$bili
#> [1] "bili"
#>
#> attr(,"matrix.terms")$albumin
#> [1] "albumin"
#>
md4$timevar
#> sexf
#> 1 1
#> 2 1
#> 3 0
#> 4 1
#> attr(,"levels")
#> attr(,"levels")$sex
#> [1] "m" "f"
#>
#> attr(,"arguments")
#> attr(,"arguments")$sex
#> attr(,"arguments")$sex$test
#> [1] "0"
#>
#>
#> attr(,"arguments.terms")
#> attr(,"arguments.terms")$test
#> sexm sexf
#> "0" "0"
#>
#> attr(,"matrix.terms")
#> attr(,"matrix.terms")$sex
#> [1] "sexm" "sexf"
#>