-
Notifications
You must be signed in to change notification settings - Fork 48
/
predict.lrm.s
126 lines (121 loc) · 4.32 KB
/
predict.lrm.s
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
126
predict.lrm <- function(object, ...,
type=c("lp","fitted","fitted.ind","mean","x","data.frame",
"terms", "cterms", "ccterms", "adjto", "adjto.data.frame",
"model.frame"),
se.fit=FALSE, codes=FALSE)
{
type <- match.arg(type)
if(type %nin% c("fitted","fitted.ind", "mean"))
return(predictrms(object,...,type=type, se.fit=se.fit))
xb <- predictrms(object, ..., type="lp", se.fit=FALSE)
rnam <- names(xb)
ns <- object$non.slopes
cnam <- names(object$coef[1:ns])
trans <- object$trans
## If orm object get cumulative probability function used
cumprob <- if(length(trans)) trans$cumprob else plogis
if(se.fit)
warning('se.fit not supported with type="fitted" or type="mean"')
if(ns == 1 & type == "mean")
stop('type="mean" makes no sense with a binary response')
if(ns == 1) return(cumprob(xb))
intcept <- object$coef[1:ns]
interceptRef <- object$interceptRef
if(!length(interceptRef)) interceptRef <- 1
xb <- xb - intcept[interceptRef]
xb <- sapply(intcept, "+", xb)
P <- cumprob(xb)
nam <- names(object$freq)
if(is.matrix(P)) dimnames(P) <- list(rnam, cnam)
else names(P) <- names(object$coef[1:ns])
if(type=="fitted") return(P)
##type="mean" or "fitted.ind"
vals <- names(object$freq)
P <- matrix(P, ncol=ns)
Peq <- cbind(1, P) - cbind(P, 0)
if(type == "fitted.ind") {
ynam <- as.character(attr(object$terms, "formula")[2])
ynam <- paste(ynam, "=", vals, sep="")
dimnames(Peq) <- list(rnam, ynam)
return(drop(Peq))
}
##type="mean"
if(codes) vals <- 1:length(object$freq)
else {
vals <- as.numeric(vals)
if(any(is.na(vals)))
stop('values of response levels must be numeric for type="mean" and codes=F')
}
m <- drop(Peq %*% vals)
names(m) <- rnam
m
}
predict.orm <- function(object, ...,
type=c("lp","fitted","fitted.ind","mean","x","data.frame",
"terms", "cterms", "ccterms", "adjto", "adjto.data.frame",
"model.frame"),
se.fit=FALSE, codes=FALSE)
{
type <- match.arg(type)
predict.lrm(object, ..., type=type, se.fit=se.fit, codes=codes)
}
Mean.lrm <- function(object, codes=FALSE, ...)
{
ns <- object$non.slopes
if(ns < 2)
stop('using this function only makes sense for >2 ordered response categories')
if(codes) vals <- 1:length(object$freq)
else {
vals <- object$yunique
if(!length(vals)) vals <- names(object$freq)
vals <- as.numeric(vals)
if(any(is.na(vals)))
stop('values of response levels must be numeric for codes=FALSE')
}
f <- function(lp=numeric(0), X=numeric(0),
intercepts=numeric(0), slopes=numeric(0),
info=numeric(0), values=numeric(0),
interceptRef=integer(0), trans=trans, conf.int=0)
{
ns <- length(intercepts)
lp <- if(length(lp)) lp - intercepts[interceptRef] else matxv(X, slopes)
xb <- sapply(intercepts, '+', lp)
P <- matrix(trans$cumprob(xb), ncol = ns)
P <- cbind(1, P) - cbind(P, 0)
m <- drop(P %*% values)
names(m) <- names(lp)
if(conf.int) {
if(! length(X)) stop('must specify X if conf.int > 0')
lb <- matrix(sapply(intercepts, '+', lp), ncol = ns)
dmean.dalpha <- t(apply(trans$deriv(lb),
1, FUN=function(x)
x * (values[2:length(values)] - values[1:ns])))
dmean.dbeta <- apply(dmean.dalpha, 1, sum) * X
dmean.dtheta <- cbind(dmean.dalpha, dmean.dbeta)
mean.var <- diag(dmean.dtheta %*% solve(info, t(dmean.dtheta)))
w <- qnorm((1 + conf.int) / 2) * sqrt(mean.var)
attr(m, 'limits') <- list(lower = m - w,
upper = m + w)
}
m
}
## If lrm fit, add information that orm fits have
family <- object$family
trans <- object$trans
if(! length(family)) {
family <- 'logistic'
trans <- probabilityFamilies$logistic
}
## Re-write first derivative so that it doesn't need the f argument
if(family == "logistic")
trans$deriv <- function(x) {p <- plogis(x); p * (1. - p)}
ir <- object$interceptRef
if(!length(ir)) ir <- 1
formals(f) <- list(lp=numeric(0), X=numeric(0),
intercepts=object$coef[1 : ns],
slopes=object$coef[- (1 : ns)],
info=object$info.matrix,
values=vals, interceptRef=ir, trans=trans,
conf.int=0)
f
}