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
|
#' Compute Confidence Intervals
#'
#' Compute and display confidence intervals for model estimates. Methods are
#' provided for the mean of a numeric vector `ci.default`, the probability
#' of a binomial vector `ci.binom`, and for `lm`, `lme`, and
#' `mer` objects are provided.
#'
#'
#' @aliases ci ci.numeric ci.binom ci.lm ci.lme ci.estimable ci.fit_contrast
#' @param x object from which to compute confidence intervals.
#' @param confidence confidence level. Defaults to 0.95.
#' @param alpha type one error rate. Defaults to 1.0-`confidence`
#' @param \dots Arguments for methods
#' @return vector or matrix with one row per model parameter and
#' elements/columns `Estimate`, `CI lower`, `CI upper`,
#' `Std. Error`, `DF` (for lme objects only), and `p-value`.
#' @author Gregory R. Warnes \email{greg@@warnes.net}
#' @seealso [stats::confint()], [stats::lm()],
#' [stats::summary.lm()]
#' @keywords regression
#'
#' @examples
#'
#'
#' # mean and confidence interval
#' ci( rnorm(10) )
#'
#' # binomial proportion and exact confidence interval
#' b <- rbinom( prob=0.75, size=1, n=20 )
#' ci.binom(b) # direct call
#' class(b) <- 'binom'
#' ci(b) # indirect call
#'
#' # confidence intervals for regression parameteres
#' data(state)
#' reg <- lm(Area ~ Population, data=as.data.frame(state.x77))
#' ci(reg)
#'
#' @export
ci <- function(x, confidence=0.95,alpha=1-confidence,...)
UseMethod("ci")
#' @rdname ci
#' @param na.rm `logical` indicating whether missing values should be removed.
#' @exportS3Method gmodels::ci
#' @importFrom stats qt
#' @importFrom stats sd
ci.numeric <- function(x, confidence=0.95,alpha=1-confidence,na.rm=FALSE,...)
{
warning("No class or unkown class. Using default calcuation.")
est <- mean(x, na.rm=na.rm)
stderr <- sd(x, na.rm=na.rm)/sqrt(nobs(x));
ci.low <- est + qt(alpha/2, nobs(x)-1) * stderr
ci.high <- est - qt(alpha/2, nobs(x)-1) * stderr
retval <- c(
Estimate=est,
"CI lower"=ci.low,
"CI upper"=ci.high,
"Std. Error"=stderr
)
retval
}
#' @export ci.binom
#' @exportS3Method gmodels::ci
#' @importFrom gdata nobs
#' @importFrom stats qbeta
#' @importFrom stats qt
ci.binom <- function(x, confidence=0.95,alpha=1-confidence,...)
{
if( !(all(x %in% c(0,1))) ) stop("Binomial values must be either 0 or 1.")
if( all(x==0) || all(x==1) )
warning("All observed values are ", as.numeric(x[1]), ", so estimated Std. Error is 0.")
est <- mean(x, na.rm=TRUE)
n <- nobs(x)
x <- sum(x)
stderr <- sqrt(est*(1-est)/n)
ci.low <- qbeta( alpha/2, x , n + 1 - x)
ci.high <- qbeta(1- alpha/2, x+1, n-x )
retval <- cbind(Estimate=est,
"CI lower"=ci.low,
"CI upper"=ci.high,
"Std. Error"= stderr
)
retval
}
#' @exportS3Method gmodels::ci
#' @importFrom stats coef
ci.lm <- function(x,confidence=0.95,alpha=1-confidence,...)
{
x <- summary(x)
est <- coef(x)[,1] ;
ci.low <- est + qt(alpha/2, x$df[2]) * coef(x)[,2] ;
ci.high <- est - qt(alpha/2, x$df[2]) * coef(x)[,2] ;
retval <- cbind(Estimate=est,
"CI lower"=ci.low,
"CI upper"=ci.high,
"Std. Error"= coef(x)[,2],
"p-value" = coef(x)[,4])
retval
}
#' @exportS3Method gmodels::ci
#' @importFrom stats qt
ci.lme <- function(x,confidence=0.95,alpha=1-confidence,...)
{
x <- summary(x)
est <- x$tTable[,"Value"] ;
ci.low <- est + qt(alpha/2, x$tTable[,"DF"]) * x$tTable[,"Std.Error"] ;
ci.high <- est - qt(alpha/2, x$tTable[,"DF"]) * x$tTable[,"Std.Error"] ;
retval <- cbind(Estimate=est,
"CI lower"=ci.low,
"CI upper"=ci.high,
"Std. Error"= x$tTable[,"Std.Error"],
"DF" = x$tTable[,"DF"],
"p-value" = x$tTable[,"p-value"])
rownames(retval) <- rownames(x$tTable)
retval
}
## ci.mer <- function (x,
## confidence = 0.95,
## alpha = 1 - confidence,
## n.sim = 1e4,
## ...)
## {
## x.effects <- x@fixef
## n <- length(x.effects)
## retval <- gmodels::est.mer(obj = x,
## cm = diag(n),
## beta0 = rep(0, n),
## conf.int = confidence,
## show.beta0 = FALSE,
## n.sim = n.sim)
## retval <- retval[,
## c("Estimate", "Lower.CI", "Upper.CI", "Std. Error", "p value"),
## drop=FALSE
## ]
## colnames(retval)[c(2:3, 5)] <- c("CI lower", "CI upper", "p-value")
## rownames(retval) <- names(x.effects)
## retval
## }
#' @exportS3Method gmodels::ci
#' @importFrom stats qt
ci.estimable <- function(x,confidence=0.95,alpha=1-confidence,...)
{
ci.low <- x$Estimate + qt(alpha/2, x$DF) * x$"Std. Error"
ci.high <- x$Estimate - qt(alpha/2, x$DF) * x$"Std. Error"
retval <- cbind(Estimate=x$Estimate,
"CI lower"=ci.low,
"CI upper"=ci.high,
"Std. Error"= x$"Std. Error",
"p-value" = x$"Pr(>|t|)"
)
rownames(retval) <- rownames(x)
retval
}
#' @exportS3Method gmodels::ci
ci.fit_contrast <- function (x, confidence = 0.95, alpha = 1 - confidence, ...)
{
if( !all(c("lower CI", "upper CI") %in% colnames(x) ) )
stop("object does not contain confidence interval information.")
colnames(x) <- c("Estimate", "Std. Error", "Delete",
"p-value",
"CI lower", "CI upper")
x[, c(1, 5:6, 2, 4), drop=FALSE]
}
|
