-
Notifications
You must be signed in to change notification settings - Fork 2
/
ttesti.R
executable file
·245 lines (238 loc) · 11.6 KB
/
ttesti.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
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
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
#' T-test Given Summary Statistics with Improved Layout
#'
#' Performs a one- or two-sample t-test given summary statistics.
#' In the two-sample case, the user can specify
#' whether or not equal variances should be presumed.
#'
#' If \code{obs2}, \code{mean2}, or \code{sd2} is specified, then all three must be specified
#' and a two-sample t-test is run.
#'
#' @aliases ttesti ttesti.do print.ttesti ttesti.default
#'
#' @param obs number of observations for the
#' first sample.
#' @param mean the sample mean of the first
#' sample.
#' @param sd the sample standard deviation of
#' the first sample.
#' @param obs2 number of observations for the second sample (this is optional).
#' @param mean2 if \code{obs2} is supplied, then sample mean of the second
#' sample must be supplied.
#' @param sd2 if \code{obs2} is supplied, then sample standard deviation of the
#' second sample must be supplied.
#' @param null.hypoth a number specifying the
#' null hypothesis for the mean (or difference in means if performing a
#' two-sample test). Defaults to zero.
#' @param alternative a string: one of
#' \code{"less"}, \code{"two.sided"}, or \code{"greater"} specifying the form
#' of the test. Defaults to a two-sided test.
#' @param conf.level confidence level of the test.
#' Defaults to 0.95.
#' @param var.eq a logical value, either
#' \code{TRUE} or \code{FALSE} (default), specifying whether or not equal
#' variances should be presumed in a two-sample t-test.
#' @param more.digits a numeric value
#' specifying whether or not to display more or fewer digits in the output.
#' Non-integers are automatically rounded down.
#'
#' @return a list of class \code{ttesti}. The print method lays out the information in an easy-to-read
#' format.
#' \item{tab}{A formatted table of descriptive and inferential statistics (number of observations,
#' mean, standard error of the mean estimate, standard deviation),
#' along with a confidence interval for the mean.}
#' \item{df}{Degrees of freedom for the t-test.}
#' \item{p}{P-value for the t-test.}
#' \item{tstat}{Test statistic for the t-test.}
#' \item{par}{A vector of information about the type of test (null hypothesis, alternative hypothesis, etc.)}
#' \item{twosamp}{A logical value indicating whether a two-sample test was performed.}
#' \item{call}{The call made to the \code{ttesti} function.}
#'
#' @examples
#'
#' # t-test given sample descriptives
#' ttesti(24, 175, 35, null.hypoth=230)
#'
#' # two-sample test
#' ttesti(10, -1.6, 1.5, 30, -.7, 2.1)
#'
#' @export ttesti
ttesti <- function(obs,
mean,
sd,
obs2=NA,
mean2=NA,
sd2=NA,
null.hypoth = 0,
conf.level=.95,
alternative="two.sided",
var.eq = FALSE,
more.digits = 0){
ttesti.do <- function(obs, mean, sd,
obs2, mean2, sd2,
null.hypoth, conf.level,
alternative, var.eq, more.digits, ...) {
if (length(obs) > 1 || length(obs2) > 1 ||
length(mean) > 1 || length(mean2) > 1 ||
length(sd) > 1 || length(sd2) > 1){
stop("'obs', 'mean', and 'sd' must be variables of length 1.")
}
if (!(obs%%1 == 0) || obs<= 0 || (!is.na(obs2) && (!(obs2%%1 == 0) || obs2 <= 0))){
stop("Number of observations must be a positive integer.")
}
if (!is.scalar(mean) || (!is.na(mean2) && !is.scalar(mean2))){
stop("Mean must be scalar.")
}
if (!is.scalar(sd) || (!is.na(sd2) && !is.scalar(sd2))){
stop("SD must be scalar.")
}
if(is.na(obs2) && (!is.na(mean2) || !is.na(sd2))) {
stop("A second number of observations must be entered for two sample test")
}
if(!is.na(obs2) && (is.na(mean2) || is.na(sd2))){
stop("SD and mean for the second sample must be entered for two sample test")
}
# check that alternative is one of "less", "two.sided", or "greater"
if (!(alternative %in% c("less", "two.sided", "greater"))) {
stop("'alternative' must be either 'less', 'two.sided', or 'greater'")
}
# make sure var.eq is either true or false
if (!is.logical(var.eq)) {
stop("Please specify a logical value for variable 'var.eq'")
}
# check to make sure additional digit request is numeric
if (!is.scalar(more.digits)) {
stop("Argument 'more.digits' must be numeric")
}
if (!is.scalar(null.hypoth)){
#if (!is.numeric(null.hypoth) || length(null.hypoth) > 1 || !is.finite(null.hypoth)){
stop("Null must be a scalar.")
}
if (!(is.scalar(conf.level) && (conf.level > 0) && (conf.level < 1))){
stop("'conf.level' must be a single number between 0 and 1")
}
more.digits <- max(floor(more.digits), 0)
digits <- 3 + more.digits
cl <- conf.level+(1-conf.level)/2
twosamp <- FALSE
if(is.na(obs2)){
stdErr <- sd/sqrt(obs)
CIlower <- mean-stats::qt(conf.level+(1-conf.level)/2, obs-1)*stdErr
CIupper <- mean+stats::qt(conf.level+(1-conf.level)/2, obs-1)*stdErr
tstat <- (mean-null.hypoth)/(stdErr)
printMat <- matrix(c(obs,
format(mean, digits = digits),
format(stdErr, digits = digits),
format(sd, digits = digits),
paste("[",format(CIlower, digits = digits),", ",
format(CIupper, digits=digits), "]", sep="")), nrow=1)
colnames(printMat) <- c("Obs", "Mean", "Std. Error", "Std. Dev.", paste0(conf.level*100, "% CI"))
rownames(printMat) <- "var1"
df <- obs - 1
if(alternative == "two.sided"){
p <- 2*stats::pt(-abs(tstat), obs-1)
} else if (alternative == "less"){
p <- stats::pt(tstat, obs - 1)
} else {
p <- 1-stats::pt(tstat, obs-1)
}
} else{
twosamp <- TRUE
if(var.eq){
stdErr1 <- sd/sqrt(obs)
stdErr2 <- sd2/sqrt(obs2)
CIlower1 <- mean-stats::qt(conf.level+(1-conf.level)/2, obs-1)*stdErr1
CIupper1 <- mean+stats::qt(conf.level+(1-conf.level)/2, obs-1)*stdErr1
CIlower2 <- mean2-stats::qt(conf.level+(1-conf.level)/2, obs2-1)*stdErr2
CIupper2 <- mean2+stats::qt(conf.level+(1-conf.level)/2, obs2-1)*stdErr2
stdErrDiff <- sqrt((sd^2*(obs-1)+sd2^2*(obs2-1))/(obs+obs2-2))*sqrt(1/obs+1/obs2)
tstat <- (mean-mean2)/(stdErrDiff)
CIlowerDiff <- (mean-mean2)-stats::qt(conf.level+(1-conf.level)/2, obs+obs2-2)*stdErrDiff
CIupperDiff <- (mean-mean2)+stats::qt(conf.level+(1-conf.level)/2, obs+obs2-2)*stdErrDiff
printMat <- matrix(c(obs,
format(mean,digits = digits),
format(stdErr1, digits = digits),
format(sd, digits = digits),
paste("[",format(CIlower1, digits = digits),", ",
format(CIupper1, digits=digits), "]", sep="")), nrow=1)
colnames(printMat) <- c("Obs", "Mean", "Std. Error", "Std. Dev.", paste0(conf.level*100, "% CI"))
printMat <- rbind(printMat, c(obs2,
format(mean2,digits = digits),
format(stdErr2, digits = digits),
format(sd2, digits = digits),
paste("[",format(CIlower2, digits = digits),", ",
format(CIupper2, digits = digits), "]", sep="")))
printMat <- rbind(printMat, c(obs+obs2,
format(mean-mean2,digits = digits),
format(stdErrDiff, digits = digits), "<NA>",
paste("[",format(CIlowerDiff, digits = digits),", ",
format(CIupperDiff, digits = digits), "]", sep="")))
rownames(printMat) <- c("var1", "var2", "diff")
df <- obs+obs2-2
if(alternative == "two.sided"){
p <- 2*stats::pt(-abs(tstat), obs+obs2-2)
} else if (alternative == "less"){
p <- stats::pt(tstat, obs + obs2 - 2)
} else {
p <- 1-stats::pt(tstat, obs+obs2-2)
}
} else {
stdErr1 <- sd/sqrt(obs)
stdErr2 <- sd2/sqrt(obs2)
CIlower1 <- mean-stats::qt(conf.level+(1-conf.level)/2, obs-1)*stdErr1
CIupper1 <- mean+stats::qt(conf.level+(1-conf.level)/2, obs-1)*stdErr1
CIlower2 <- mean2-stats::qt(conf.level+(1-conf.level)/2, obs2-1)*stdErr2
CIupper2 <- mean2+stats::qt(conf.level+(1-conf.level)/2, obs2-1)*stdErr2
stdErrDiff <- sqrt(sd^2/(obs)+sd2^2/(obs2))
tstat <- (mean-mean2)/(stdErrDiff)
df <- (sd^2/obs+sd2^2/obs2)^2/((sd^2/obs)^2/(obs-1)+(sd2^2/obs2)^2/(obs2-1))
CIlowerDiff <- (mean-mean2)-stats::qt(conf.level+(1-conf.level)/2, obs+obs2-2)*stdErrDiff
CIupperDiff <- (mean-mean2)+stats::qt(conf.level+(1-conf.level)/2, obs+obs2-2)*stdErrDiff
printMat <- matrix(c(obs,
format(mean, digits =digits),
format(stdErr1, digits =digits),
format(sd, digits =digits),
paste("[",format(CIlower1, digits =digits),", ",
format(CIupper1, digits =digits), "]", sep="")), nrow=1)
colnames(printMat) <- c("Obs", "Mean", "Std. Error", "Std. Dev.", paste0(conf.level*100, "% CI"))
printMat <- rbind(printMat,c(obs2,
format(mean2, digits =digits),
format(stdErr2, digits =digits),
format(sd2, digits =digits),
paste("[",format(CIlower2, digits = digits),", ",
format(CIupper2, digits =digits), "]", sep="")))
printMat <- rbind(printMat,
c(obs+obs2, format(mean-mean2, digits =digits),
format(stdErrDiff, digits =digits), "<NA>",
paste("[",format(CIlowerDiff, digits =digits),", ",
format(CIupperDiff, digits =digits), "]", sep="")))
rownames(printMat) <- c("var1", "var2", "diff")
if(alternative == "two.sided"){
p <- 2*stats::pt(-abs(tstat), df)
} else if (alternative == "less"){
p <- stats::pt(tstat, df)
} else {
p <- 1-stats::pt(tstat, df)
}
}
}
tstat <- as.numeric(format(tstat, digits = max(digits +
1, 3)))
p <- as.numeric(format(p, digits = max(digits,
digits + 3)))
df <- as.numeric(format(df, digits = max(digits,
digits + 3)))
par <- c(null.hypoth = null.hypoth, alternative = alternative,
var.eq = var.eq, conf.level = conf.level, twosamp = twosamp,
digits = digits)
invisible(list(tab = printMat, df = df, p = p, tstat = tstat,
par = par))
}
myargs <- c(deparse(substitute(obs)), deparse(substitute(obs2)))
ttesti.obj <- ttesti.do(obs = obs, mean = mean, sd = sd,
obs2 = obs2, mean2 = mean2, sd2 = sd2,
null.hypoth = null.hypoth, conf.level = conf.level,
alternative = alternative, var.eq = var.eq, more.digits = more.digits)
ttesti.obj$call <- match.call()
class(ttesti.obj) <- "ttesti"
return(ttesti.obj)
}