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
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
|
# glmnet.R: plotmo functions for glmnet and glmnetUtils objects
plotmo.prolog.glmnet <- function(object, object.name, trace, ...) # invoked when plotmo starts
{
# save (possibly user specified) s for use by plot_glmnet and predict.glmnet
s <- dota("predict.s", ...) # get predict.s from dots, NA if not in dots
if(is.na(s))
s <- dota("s", ...) # get s from dots, NA if not in dots
if(is.na(s))
s <- 0 # unspecified, default to match plotmo.predict.glmnet
check.numeric.scalar(s)
attr(object, "plotmo.s") <- s
object
}
plotmo.predict.glmnet <- function(object, newdata, type, ..., TRACE)
{
s <- attr(object, "plotmo.s") # get the predict.glmnet s
stopifnot(!is.null(s)) # uninitialized?
# newx for predict.glmnet must be a matrix not a dataframe,
# so here we use plotmo.predict.defaultm (not plotmo.predict.default)
yhat <- plotmo.predict.defaultm(object, newdata, type=type, force.s=s,
..., TRACE=TRACE)
if(length(dim(yhat) == 2) && NCOL(yhat) == 1)
colnames(yhat) <- paste0("s=", signif(s,2)) # so s=.12 appears in plot title
yhat
}
plotmo.predict.glmnet.formula <- function(object, newdata, type, ..., TRACE) # glmnetUtils package
{
# same as plotmo.predict.glmnet but doesn't convert newx to a matrix
s <- attr(object, "plotmo.s") # get the predict.glmnet s
stopifnot(!is.null(s)) # uninitialized?
yhat <- plotmo.predict.default(object, newdata, type=type, force.s=s,
..., TRACE=TRACE)
if(length(dim(yhat) == 2) && NCOL(yhat) == 1)
colnames(yhat) <- paste0("s=", signif(s,2)) # so s=.12 appears in plot title
yhat
}
plotmo.singles.glmnet <- function(object, x, nresponse, trace, all1, ...)
{
# return the indices of the 25 biggest coefs, but exclude zero coefs
s <- attr(object, "plotmo.s") # get the predict.glmnet s
stopifnot(!is.null(s)) # uninitialized?
lambda.index <- which.min(abs(object$lambda - s)) # index into object$lambda
trace2(trace, "plotmo.singles.glmnet: s %g lambda.index %g\n", s, lambda.index)
beta <- object$beta
if(is.list(beta)) { # multiple response model?
check.integer.scalar(nresponse, min=1, max=length(beta))
beta <- beta[[nresponse]]
}
beta <- as.vector(beta[, lambda.index]) # as.vector converts from dgCMatrix
order <- order(abs(beta), decreasing=TRUE)
max.nsingles <- if(all1) Inf else 25
# extract the biggest coefs
beta <- beta[order][1:min(max.nsingles, length(beta))]
nsingles <- sum(abs(beta) > 1e-8) # drop zero coefs
order[seq_len(nsingles)]
}
plotmo.prolog.cv.glmnet <- function(object, object.name, trace, ...) # invoked when plotmo starts
{
# cv.glmnet objects don't have their call field in the usual place,
# so fix that (tested on glmnet version 2.0-2).
# Note that getCall() doesn't work on cv.glmnet objects.
if(is.null(object[["call"]])) {
object$call <- object$glmnet.fit$call
stopifnot(!is.null(object$call), is.call(object$call))
}
# save (possibly user specified) s for use by plot_glmnet and predict.glmnet
s <- dota("predict.s", ...) # get predict.s from dots, NA if not in dots
if(is.na(s))
s <- dota("s", ...) # get s from dots, NA if not in dots
if(is.na(s))
s <- "lambda.1se" # unspecified, default to match predict.cv.glmnet
s <- match.choices(s, c("lambda.1se", "lambda.min"), "s")
attr(object, "plotmo.s") <- s
object
}
plotmo.predict.cv.glmnet <- function(object, newdata, type, ..., TRACE)
{
s <- attr(object, "plotmo.s") # get the predict.glmnet s
stopifnot(!is.null(s)) # uninitialized?
if(inherits(object, "cv.glmnet.formula")) { # glmnetUtils package
yhat <- plotmo.predict.default(object, newdata, type=type, force.s=s,
..., TRACE=TRACE)
} else { # glmnet package
# newx for predict.glmnet must be a matrix not a dataframe,
# so here we use plotmo.predict.defaultm (not plotmo.predict.default)
yhat <- plotmo.predict.defaultm(object, newdata, type=type, force.s=s,
..., TRACE=TRACE)
}
if(length(dim(yhat) == 2) && NCOL(yhat) == 1)
colnames(yhat) <- paste0("s=\"", s, "\"") # so s="lambda.1se" appears in plot title
yhat
}
plotmo.singles.cv.glmnet <- function(object, x, nresponse, trace, all1, ...)
{
# return the indices of the 25 biggest coefs, but exclude zero coefs
s <- attr(object, "plotmo.s") # get the predict.glmnet s
beta <- coef(object, s=s)
if(is.list(beta)) { # multiple response model?
check.integer.scalar(nresponse, min=1, max=length(beta))
beta <- beta[[nresponse]]
}
beta <- as.vector(beta) # as.vector converts from dgCMatrix
beta <- beta[-1] # drop intercept
order <- order(abs(beta), decreasing=TRUE)
max.nsingles <- if(all1) Inf else 25
# extract the biggest coefs
beta <- beta[order][1:min(max.nsingles, length(beta))]
nsingles <- sum(abs(beta) > 1e-8) # drop zero coefs
order[seq_len(nsingles)]
}
# glmnet family="binomial", y is a vector of 1s and 2s.
# convert 1s and 2s to 0s and 1s to match predicted values
plotmo.y.lognet <- function(object, trace, naked, expected.len, nresponse, ...)
{
# plotmo.y.default returns list(field=y, do.subset=do.subset)
list <- plotmo.y.default(object, trace, naked, expected.len)
# following is needed for glmnetUtils:glmnet.formula models (but not for glmnet xy models)
if(is.data.frame(list$field))
list$field <- list$field[[1]]
stopifnot(!is.null(list$field)) # paranoia
list$do.subset <- FALSE # glmnet doesn't support subset so don't even try
# TODO following only works correctly if default ordering of factor was used?
list$field <- as.numeric(list$field) # as.numeric needed if y is a factor
list$field - min(list$field) # convert 1s and 2s to 0s and 1s
}
# glmnet family="multinomial"
plotmo.y.multnet <- function(object, trace, naked, expected.len, nresponse, ...)
{
# plotmo.y.default returns list(field=y, do.subset=do.subset)
list <- plotmo.y.default(object, trace, naked, expected.len)
list$do.subset <- FALSE # glmnet doesn't support subset so don't even try
if(is.null(nresponse)) # plotmo uses nresponse=NULL in initial checking
nresponse <- 1
if(NCOL(list$field) > 1) # if y is multiple columns assume it's an indicator matrix
y <- list$field
else { # else convert it to an indicator matrix
# TODO following only works correctly if default ordering of factor was used?
y1 <- as.numeric(list$field) # as.numeric needed if y is a factor
stopifnot(min(y1) == 1 && max(y1) > 1) # sanity check
# convert y1 to an indicator matrix of 0s and 1s (NA_real_ to avoid type convert)
y <- matrix(NA_real_, nrow=length(y1), ncol=max(y1))
for(i in 1:max(y1))
y[,i] <- as.numeric(y1 == nresponse)
}
y
}
# glmnet family="mgaussian"
plotmo.y.mrelnet <- function(object, trace, naked, expected.len, nresponse, ...)
{
plotmo.y.multnet(object, trace, naked, expected.len, nresponse, ...)
}
|
