-
Notifications
You must be signed in to change notification settings - Fork 1
/
summarymethods.R
145 lines (128 loc) · 4.41 KB
/
summarymethods.R
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
#' Summary method for dcalasso objects
#'
#' \code{summary.dcalasso} summarizes output of dcalasso fit
#' @param object a dcalasso fit
#' @param unpen whether to print out the unpenalized result
#' @param ... ...
#' @export
#' @author Yan Wang, Tianxi Cai, Chuan Hong
summary.dcalasso = function(object, unpen = F, ...){
stopifnot(inherits(object, "dcalasso"))
coef.u <- object$coefficients.unpen
coef.p <- object$coefficients.pen
s.err.u <- sqrt(diag(object$cov.unpen))
s.err.p <- sqrt(diag(object$cov.pen))
tvalue.u <- coef.u/s.err.u
tvalue.p <- coef.p/s.err.p
coef.table.p <- cbind(coef.p, s.err.p, tvalue.p, 2 * pnorm(-abs(tvalue.p)))
coef.table.u <- cbind(coef.u, s.err.u, tvalue.u, 2 * pnorm(-abs(tvalue.u)))
colnames(coef.table.p) = c("Penalized Est", "Std. Error", "z value", "Pr(>|z|)")
colnames(coef.table.u) = c("Unpenalized Est", "Std. Error", "z value", "Pr(>|z|)")
ans = list(coef.table.p = coef.table.p, coef.table.u = coef.table.u,
cov.pen = object$cov.pen, cov.unpen = object$cov.unpen,
K = object$K, n = object$n, n.pen = object$n.pen, family = object$family$family,
iter = object$iter, BIC = object$BIC.opt, lambda = object$lambda.opt)
class(ans) <- 'summary.dcalasso'
return(ans)
}
#' Print summary for dcalasso objects
#'
#' \code{print.summary.dcalasso} summarizes output of dcalasso fit
#'
#' @param x a summary.dcalasso object
#' @param ... ...
# @param unpen whether to print out the unpenalized result
#' @export
#' @author Yan Wang, Tianxi Cai, Chuan Hong
print.summary.dcalasso = function(x, ...){
unpen = F
stopifnot(inherits(x, "summary.dcalasso"))
cat("\nDivide-and-conquer adaptive lasso for a ", x$family," model, n=",
x$n,".\n")
cat("\nInitial estimator computed for K=",x$K, "and one-step estimation with",
x$iter,"iterations.\n")
cat("\nPenalized summary:\n")
printCoefmat(x$coef.table.p, signif.stars=T, na.print = ".")
if (unpen){
cat("\nUnpenalized summary:\n")
printCoefmat(x$coef.table.u, signif.stars=T, na.print = "NA")
}
cat("\nBIC = ", x$BIC," with lambda = ", x$lambda,"\n")
}
#' Print dcalasso objects
#'
#' \code{print.dcalasso} summarizes output of dcalasso fit
#'
#' @param x a dcalasso object
#' @param ... ...
#' @export
#' @author Yan Wang, Tianxi Cai, Chuan Hong
print.dcalasso = function(x, ...){
summary.dcalasso(x)
}
#' Plot BIC paths for dcalasso objects
#'
#' \code{plot.dcalasso} summarizes output of dcalasso fit
#'
#' @param x a dcalasso object
#' @param ... ...
#' @export
#' @author Yan Wang, Tianxi Cai, Chuan Hong
plot.dcalasso = function(x, ...){
plot(log10(x$lambda), x$BIC, xlab = 'Log10 lambda', ylab = 'BIC')
}
#' Extract variance covariance from a dcalasso objects
#'
#' \code{vcov.dcalasso} extracts variance covariance objects
#'
#' @param object a dcalasso object
#' @param unpen whether to switch to the unpenalized variance covariance
#' @param ... ...
#' @export
#' @author Yan Wang, Tianxi Cai, Chuan Hong
vcov.dcalasso = function(object, unpen=F, ...){
if (unpen){
return(object$cov.unpen)
}else{
return(object$cov.pen)
}
}
#' Extract coefficients from dcalasso objects
#'
#' \code{coef.dcalasso} extracts coefficients from dcalasso objects
#'
#' @param object a dcalasso object
#' @param unpen whether to switch to the unpenalized coefficients
#' @param ... ...
#' @export
#' @author Yan Wang, Tianxi Cai, Chuan Hong
coef.dcalasso = function(object, unpen=F, ...){
if (unpen){
return(object$coefficients.unpen)
}else{
return(object$coefficients.pen)
}
}
#' Prediction of dcalasso object
#'
#' \code{predict.dcalasso} makes prediction of a dcalasso object based on the adaptive lasso estimation.
#' @param object a dcalasso object
#' @param newdata a new data frame
#' @param type "terms", "link", "response" same as predict.glm
#' @export
#' @author Yan Wang, Tianxi Cai, Chuan Hong
predict.dcalasso = function(object, newdata, type = 'link'){
Terms = delete.response(object$Terms)
m = model.frame(Terms, newdata)
X = model.matrix(Terms, m)
X = X[,names(object$coefficients.pen)]
if (type == 'terms'){
return(list(fit = X * object$coefficients.pen,
se.fit = X * sqrt(diag(object$cov.unpen))))
}else if (type == 'link'){
return(list(fit = X %*% object$coefficients.pen,
se.fit = X %*% sqrt(diag(object$cov.unpen))))
}else if (type == 'response'){
return(fit = object$family$linkinv(X %*% object$coefficients.pen))
}
}