1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
|
#
# Determine which 'canonical variables' depend on a supplied covariate
#
# $Revision: 1.11 $ $Date: 2021/03/29 08:11:44 $
#
model.depends <- function(object) {
## supplied covariates
fo <- formula(object)
if(length(as.list(fo)) == 3) {
## formula has a response: strip it
fo <- fo[-2]
}
covars <- variablesinformula(fo)
mm <- model.matrix(object)
depends <- matrix(FALSE, ncol(mm), length(covars),
dimnames=list(colnames(mm), covars))
## model term labels
tt <- terms(object)
lab <- attr(tt, "term.labels")
## map from canonical covariates to term labels
ass <- attr(mm, "assign") %orifnull% object[["assign"]]
## determine which canonical covariate depends on which supplied covariate
if(length(ass) == ncol(mm)) {
for(i in seq(along=ass)) {
if(ass[i] == 0) # 0 is the intercept term
depends[i,] <- FALSE
else {
turm <- lab[ass[i]]
depends[i, ] <- covars %in% all.vars(parse(text=turm))
}
}
} else {
warning("model.depends: unable to determine the dependence structure",
call.=FALSE)
}
## detect offsets
if(!is.null(oo <- attr(tt, "offset")) && ((noo <- length(oo)) > 0)) {
# entries of 'oo' index the list of variables in terms object
vv <- attr(tt, "variables")
offdep <- matrix(FALSE, noo, length(covars))
offnms <- character(noo)
for(i in seq_len(noo)) {
offseti <- languageEl(vv, oo[i] + 1)
offdep[i, ] <- covars %in% all.vars(offseti)
offnms[i] <- deparse(offseti)
}
rownames(offdep) <- offnms
colnames(offdep) <- covars
attr(depends, "offset") <- offdep
}
return(depends)
}
model.is.additive <- function(object) {
dep <- model.depends(object)
hit <- t(dep) %*% dep
diag(hit) <- 0
ok <- all(hit == 0)
return(ok)
}
model.covariates <- function(object, fitted=TRUE, offset=TRUE) {
md <- model.depends(object)
nm <- colnames(md)
keep <- rep.int(FALSE, length(nm))
# variables used in formula with coefficients
if(fitted) keep <- apply(md, 2, any)
# variables used in offset
if(offset) {
oo <- attr(md, "offset")
if(!is.null(oo))
keep <- keep | apply(oo, 2, any)
}
return(nm[keep])
}
has.offset.term <- function(object) {
# model terms
tt <- terms(object)
oo <- attr(tt, "offset")
return(!is.null(oo) && (length(oo) > 0))
}
has.offset <- function(object) {
has.offset.term(object) || !is.null(model.offset(model.frame(object)))
}
check.separable <- function(dmat, covname, isconstant, fatal=TRUE) {
#' Determine whether the effect of 'covname' is separable from other terms.
#' dmat = model.depends(model)
#' Find covariates entangled with 'covname' in the model
relevant <- dmat[, covname]
othercov <- (colnames(dmat) != covname)
conflict <- dmat[relevant, othercov, drop=FALSE]
if(!any(conflict)) return(TRUE)
#' names of entangled covariates
entangled <- colnames(conflict)[matcolany(conflict)]
#' not problematic if constant
if(is.null(names(isconstant))) names(isconstant) <- colnames(dmat)
ok <- unlist(isconstant[entangled])
conflict[ , ok] <- FALSE
if(!any(conflict)) return(TRUE)
#' there are conflicts
conflictterms <- matrowany(conflict)
conflictcovs <- matcolany(conflict)
whinge <- paste("The covariate", sQuote(covname),
"cannot be separated from the",
ngettext(sum(conflictcovs), "covariate", "covariates"),
commasep(sQuote(colnames(conflict)[conflictcovs])),
"in the model",
ngettext(sum(conflictterms), "term", "terms"),
commasep(sQuote(rownames(conflict)[conflictterms])))
if(fatal) stop(whinge, call.=FALSE)
warning(whinge, call.=FALSE)
return(FALSE)
}
|
