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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
|
### sCorrect-summary.glht2.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: maj 2 2018 (09:20)
## Version:
## Last-Updated: jan 24 2022 (11:11)
## By: Brice Ozenne
## Update #: 228
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
#' @title Outcome of Linear Hypothesis Testing
#' @description Estimates, p-values, and confidence intevals for linear hypothesis testing, possibly adjusted for multiple comparisons.
#'
#' @param object a \code{glht2} object.
#' @param confint [logical] should confidence intervals be output
#' @param conf.level [numeric 0-1] level of the confidence intervals.
#' @param transform [function] function to backtransform the estimates, standard errors, null hypothesis, and the associated confidence intervals
#' (e.g. \code{exp} if the outcomes have been log-transformed).
#' @param seed [integer] value that will be set before adjustment for multiple comparisons to ensure reproducible results.
#' Can also be \code{NULL}: in such a case no seed is set.
#' @param rowname.rhs [logical] when naming the hypotheses, add the right-hand side (i.e. "X1-X2=0" instead of "X1-X2").
#' @param ... argument passed to \code{multcomp:::summary.glht}, e.g. argument \code{test} to choose the type of adjustment for multiple comparisons.
#'
## * summary.glht2
#' @export
summary.glht2 <- function(object, confint = TRUE, conf.level = 0.95, transform = NULL, seed = NULL, rowname.rhs = TRUE, ...){
if(!is.null(seed)){
old.seed <- get0(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
on.exit( assign(".Random.seed", old.seed, envir = .GlobalEnv, inherits = FALSE) )
set.seed(seed)
}
keep.class <- class(object)
object$test <- NULL
object$confint <- NULL
class(object) <- setdiff(keep.class, "glht2")
keep.df <- object$df
test.df <- any( (keep.df>0) * (!is.infinite(keep.df)) == 1 )
object$df <- round(stats::median(object$df))
output <- summary(object, ...)
## restaure df when possible
method.adjust <- output$test$type
if(NROW(object$linfct)==1){method.adjust <- "none"}
if(test.df && method.adjust %in% c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none","univariate")){
output$df <- keep.df
output$test$pvalues <- stats::p.adjust(2*(1-stats::pt(abs(output$test$tstat), df = keep.df)), method = method.adjust)
}
name.hypo <- rownames(output$linfct)
n.hypo <- length(name.hypo)
if(confint && method.adjust %in% c("univariate","none","bonferroni","single-step")){
if(method.adjust %in% c("none","univariate","bonferroni")){
alpha <- switch(method.adjust,
"none" = 1-conf.level,
"univariate" = 1-conf.level,
"bonferroni" = (1-conf.level)/n.hypo)
if(test.df){
q <- stats::qt(1-alpha/2, df = output$df)
}else{
q <- stats::qnorm(1-alpha/2)
}
output$confint <- data.frame(matrix(NA, ncol = 3, nrow = n.hypo,
dimnames = list(name.hypo, c("Estimate","lwr","upr"))))
output$confint$Estimate <- as.double(output$test$coef)
output$confint$lwr <- as.double(output$test$coef - q * output$test$sigma)
output$confint$upr <- as.double(output$test$coef + q * output$test$sigma)
## range(confint(output, level = 1-alpha, calpha = univariate_calpha())$confint-output$confint)
}else if(method.adjust == "single-step"){
output <- confint(output, level = conf.level, calpha = multcomp::adjusted_calpha())
}else{
output$confint <- matrix(NA, nrow = n.hypo, ncol = 3,
dimnames = list(name.hypo, c("Estimate","lwr","upr")))
}
}
if(rowname.rhs){
table2.rownames <- paste0(name.hypo, " == ", output$rhs)
}else{
table2.rownames <- name.hypo
}
output$table2 <- data.frame(matrix(NA, nrow = n.hypo, ncol = 7,
dimnames = list(table2.rownames,
c("estimate","se","df","lower","upper","statistic","p.value"))
), stringsAsFactors = FALSE)
output$table2$estimate <- output$test$coefficients
output$table2$se <- output$test$sigma
output$table2$df <- output$df
output$table2$df[output$table2$df==0] <- Inf
output$table2$lower <- output$confint[,"lwr"]
output$table2$upper <- output$confint[,"upr"]
output$table2$statistic <- output$test$tstat
output$table2$p.value <- output$test$pvalues
output$seed <- seed
## ** transformation
output$transform <- transform
output$table2 <- transformSummaryTable(output$table2,
transform = transform)
## ** export
class(output) <- append(c("summary.glht2","summary.glht"),keep.class)
return(output)
}
## * print.summary.glht2
#' @export
print.summary.glht2 <- function(x,
digits = max(3L, getOption("digits") - 2L),
digits.p.value = max(3L, getOption("digits") - 2L),
columns = c("estimate","se","df","lower","upper","statistic","p.value"),
...){
columns <- match.arg(columns, choices = c("estimate","se","df","lower","upper","statistic","p.value"), several.ok = TRUE)
type <- x$type
call <- if(isS4(x$model)){x$model@call}else{x$model$call}
alternative <- x$alternativ
type <- x$test$type
txt.type <- switch(type,
"univariate" = "(CIs/p-values not adjusted for multiple comparisons)",
"none" = "(CIs/p-values not adjusted for multiple comparisons)",
"single-step" = paste0("(CIs/p-values adjusted for multiple comparisons -- single step max-test)"),
"free" = paste0("(CIs/p-values adjusted for multiple comparisons -- step down max-test)"),
"Westfall" = paste0("(CIs/p-values adjusted for multiple comparisons -- step down max-test with logical restrictions)"),
paste0("(CIs/p-values adjusted for multiple comparisons -- ", type, " method)")
)
txt.robust <- switch(as.character(x$robust),
"TRUE" = "Robust",
"FALSE" = "Model-based"
)
## txt.correction <- switch(as.character(x$ssc),
## "Cox" = " corrected for small sample bias (Cox correction)",
## "residuals" = " corrected for small sample bias (residual correction)",
## "NA" = ""
## )
txt.alternative <- switch(alternative,
"less" = "one sided tests - inferiority",
"greater" = "one sided tests - superiority",
"two.sided" = "two sided tests")
## display
cat("\n\t", "Simultaneous Tests for General Linear Hypotheses\n\n")
if (!is.null(type)) {
cat("Multiple Comparisons of Means (",txt.alternative,") \n\n", sep = "")
}
if (!is.null(call)) {
cat("Fit: ")
print(call)
cat("Standard errors: ",txt.robust,"\n",sep="")
cat("\n")
}
cat("Linear Hypotheses:\n")
stats::printCoefmat(x$table2[,columns[columns %in% names(x$table2)],drop=FALSE], digits = digits,
has.Pvalue = "p.value" %in% columns,
P.values = "p.value" %in% columns,
eps.Pvalue = 10^{-digits.p.value})
if(NROW(x$table2)>1){
cat(txt.type,"\n")
}
error <- attr(x$test$pvalues,"error")
if(!is.null(error) && error > 1e-12 && "p.value" %in% columns){
txt.error <- paste0("Error when computing the adjusted p-value by numerical integration: ", signif(error, digits = digits))
if(!is.null(x$seed)){
txt.error <- paste0(txt.error," (seed ",x$seed,")")
}
cat(txt.error,"\n")
}
if(!is.null(x$global)){
cat("\nGlobal test: p.value=",format.pval(x$global["p.value"], digits = digits, eps = 10^(-digits.p.value)),
" (statistic=",round(x$global["statistic"], digits = digits),
", df=",round(x$global["df"], digits = digits),")\n",sep="")
}
## if(nchar(txt.correction)>0){cat("(",txt.correction,")\n",sep="")}
cat("\n")
return(invisible(x))
}
######################################################################
### sCorrect-summary.glht2.R ends here
|
