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
118
119
120
121
122
123
124
125
|
### oh dear!
### Cox models don't have any intercept ...
model.matrix.coxph <- function(object, ...)
{
mm <- model.matrix.default(object)
at <- attributes(mm)
mm <- mm[,-1]
at$dim[2] <- at$dim[2] - 1
at$dimnames[[2]] <- at$dimnames[[2]][-1]
at$assign <- at$assign[-1]
attributes(mm) <- at
mm
}
### some methods of lmer objects
model.matrix.lmer <- function(object, ...) {
x <- object@X
if (is.null(x))
stop("models of class ", sQuote("lmer"),
" need to be fitted with argument ",
sQuote("model = TRUE"))
x
}
model.frame.lmer <- function(object, ...) {
x <- object@frame
if (is.null(x))
stop("models of class ", sQuote("lmer"),
" need to be fitted with argument ",
sQuote("model = TRUE"))
x
}
terms.lmer <- function(object, ...) {
x <- object@terms
if (is.null(x))
stop("models of class ", sQuote("lmer"),
" need to be fitted with argument ",
sQuote("model = TRUE"))
x
}
### extract coefficients, covariance matrix and
### degrees of freedom (if available) from `model'
modelparm <- function(model, coef., vcov., df, ...)
UseMethod("modelparm")
modelparm.default <- function(model, coef. = coef, vcov. = vcov,
df = NULL, ...)
{
### extract coefficients and their covariance matrix
beta <- try(coef.(model))
if (inherits(beta, "try-error"))
stop("no ", sQuote("coef"), " method for ",
sQuote("model"), " found!")
sigma <- try(vcov.(model))
if (inherits(sigma, "try-error"))
stop("no ", sQuote("vcov"), " method for ",
sQuote("model"), " found!")
sigma <- as.matrix(sigma)
if (any(length(beta) != dim(sigma)))
stop("dimensions of coefficients and covariance matrix don't match")
### determine degrees of freedom
if (is.null(df)) {
df <- 0
### check if a linear model was supplied
if (class(model)[1] %in% c("aov", "lm")) {
class(model) <- "lm"
df <- summary(model)$df[2]
}
} else {
if (df < 0) stop(sQuote("df"), " is not positive")
}
RET <- list(coef = beta, vcov = sigma, df = df)
class(RET) <- "modelparm"
RET
}
### linear mixed effects models (package lme4)
coeflmer <- function(object, ...) {
x <- object@fixef
names(x) <- rownames(vcov(object))
x
}
modelparm.lmer <- function(model, coef. = coeflmer, vcov. = vcov, df = NULL, ...)
modelparm.default(model, coef. = coef., vcov. = vcov., df = df, ...)
### survreg models (package survival)
vcovsurvreg <- function(object, ...) {
sigma <- vcov(object)
p <- length(coef(object))
return(sigma[1:p, 1:p])
}
modelparm.survreg <- function(model, coef. = coef, vcov. = vcovsurvreg, df = NULL, ...)
modelparm.default(model, coef. = coef., vcov. = vcov., df = df, ...)
### modified from package MASS
MPinv <- function (X, tol = sqrt(.Machine$double.eps))
{
if (length(dim(X)) > 2 || !(is.numeric(X) || is.complex(X)))
stop("X must be a numeric or complex matrix")
if (!is.matrix(X))
X <- as.matrix(X)
Xsvd <- svd(X)
if (is.complex(X))
Xsvd$u <- Conj(Xsvd$u)
Positive <- Xsvd$d > max(tol * Xsvd$d[1], 0)
if (all(Positive))
RET <- Xsvd$v %*% (1/Xsvd$d * t(Xsvd$u))
else if (!any(Positive))
RET <- array(0, dim(X)[2:1])
else RET <- Xsvd$v[, Positive, drop = FALSE] %*% ((1/Xsvd$d[Positive]) *
t(Xsvd$u[, Positive, drop = FALSE]))
return(list(MPinv = RET, rank = sum(Positive)))
}
|
