Comments (3)
bootstrap-t method for one-sample test
trimcibt <-
function(x,
tr = .2,
alpha = .05,
nboot = 599,
side = TRUE,
plotit = FALSE,
op = 1,
nullval = 0,
SEED = TRUE,
prCRIT = FALSE,
pr = TRUE,
xlab = "") {
#
# Compute a 1-alpha confidence interval for the trimmed mean
# using a bootstrap percentile t method.
#
# The default amount of trimming is tr=.2
# side=T, for true, indicates the symmetric two-sided method
#
# Side=F yields an equal-tailed confidence interval
#
# NOTE: p.value is reported when side=T only.
#
# prCRIT=TRUE, if side=FALSE, lower and upper critical values are returned.
# plotit=TRUE, plots an estimate of the null distribution of the Tukey--McLaughlin test statistic.
# op=1, use adaptive kernel density estimator when plotit=TRUE
# op=2, uses the expected frequency curve, which will have a lower execution if nboot is large.
#
x = WRS2:::elimna(x)
side <- as.logical(side)
p.value <- NA
if (SEED)
set.seed(2) # set seed of random number generator so that
# results can be duplicated.
test <- (mean(x, tr) - nullval) / WRS2::trimse(x, tr)
data <-
matrix(sample(x, size = length(x) * nboot, replace = TRUE), nrow = nboot)
data <- data - mean(x, tr)
top <- apply(data, 1, mean, tr)
bot <- apply(data, 1, WRS2::trimse, tr)
tval <- top / bot
if (plotit) {
if (op == 1)
akerd(tval, xlab = xlab)
if (op == 2)
rdplot(tval, xlab = xlab)
}
if (side)
tval <- abs(tval)
tval <- sort(tval)
icrit <- round((1 - alpha) * nboot)
ibot <- round(alpha * nboot / 2)
itop <-
nboot - ibot #altered code very slightly to correspond to recent versions of my books.
if (!side) {
if (prCRIT)
print(paste("Lower crit=", tval[ibot], "Upper crit=", tval[itop]))
if (prCRIT)
print(paste(".025 Lower Type I=", mean(tval <= 0 - 1.96)))
if (prCRIT)
print(paste(".05 Lower Type I=", mean(tval <= 0 - 1.645)))
if (prCRIT)
print(paste(".025 Upper Type I=", mean(tval >= 1.96)))
if (prCRIT)
print(paste(".05 Upper Type I=", mean(tval >= 1.645)))
trimcibt <- mean(x, tr) - tval[itop] * WRS2::trimse(x, tr)
trimcibt[2] <- mean(x, tr) - tval[ibot] * WRS2::trimse(x, tr)
if (pr)
print("NOTE: p.value is computed only when side=T")
}
if (side) {
if (prCRIT)
print(paste("Symmetric Crit.val=", tval[icrit]))
trimcibt <- mean(x, tr) - tval[icrit] * WRS2::trimse(x, tr)
trimcibt[2] <- mean(x, tr) + tval[icrit] * WRS2::trimse(x, tr)
p.value <- (sum(abs(test) <= abs(tval))) / nboot
}
list(
estimate = mean(x, tr),
ci = trimcibt,
test.stat = test,
p.value = p.value,
n = length(x)
)
}
trimcibt(mtcars$wt, nullval = 3)
#> $estimate
#> [1] 3.197
#>
#> $ci
#> [1] 2.841377 3.552623
#>
#> $test.stat
#> [1] 1.179181
#>
#> $p.value
#> [1] 0.2537563
#>
#> $n
#> [1] 32
Created on 2020-04-04 by the reprex package (v0.3.0.9001)
from statsexpressions.
Simpler version of the same function:
trimcibt <- function(x,
tr = 0.2,
alpha = 0.05,
nboot = 200,
nullval = 0,
...) {
test <- (mean(x, tr) - nullval) / WRS2::trimse(x, tr)
data <- matrix(sample(x, size = length(x) * nboot, replace = TRUE), nrow = nboot) - mean(x, tr)
top <- apply(data, 1, mean, tr)
bot <- apply(data, 1, WRS2::trimse, tr)
tval <- sort(abs(top / bot))
icrit <- round((1 - alpha) * nboot)
ibot <- round(alpha * nboot / 2)
itop <- nboot - ibot
trimcibt <- mean(x, tr) - tval[icrit] * WRS2::trimse(x, tr)
trimcibt[2] <- mean(x, tr) + tval[icrit] * WRS2::trimse(x, tr)
p.value <- (sum(abs(test) <= abs(tval))) / nboot
list(
estimate = mean(x, tr),
ci = trimcibt,
test.stat = test,
p.value = p.value,
n = length(x)
)
}
trimcibt(mtcars$wt, nullval = 3)
#> $estimate
#> [1] 3.197
#>
#> $ci
#> [1] 2.857201 3.536799
#>
#> $test.stat
#> [1] 1.179181
#>
#> $p.value
#> [1] 0.275
#>
#> $n
#> [1] 32
Created on 2021-01-09 by the reprex package (v0.3.0.9001)
from statsexpressions.
nboot
argument needs to be added.
(This is painfully slow function)
library(WRS2)
wmcpAKP <- function(x, tr = .2) {
# Compute Algina et al measure of effect size for pairs of J dependent groups
if (is.matrix(x) || is.data.frame(x)) x <- WRS2:::listm(x)
J <- length(x)
C <- (J^2 - J) / 2
A <- matrix(NA, nrow = C, ncol = 5)
dimnames(A) <- list(NULL, c("Group", "Group", "Effect.Size", "ci.low", "ci.up"))
ic <- 0
for (j in 1:J) {
for (k in 1:J) {
if (j < k) {
ic <- ic + 1
A[ic, 1] <- j
A[ic, 2] <- k
A[ic, 3] <- WRS2::dep.effect(x[[j]], x[[k]], tr = tr)[5]
A[ic, 4] <- WRS2::dep.effect(x[[j]], x[[k]], tr = tr)[21]
A[ic, 5] <- WRS2::dep.effect(x[[j]], x[[k]], tr = tr)[25]
}
}
}
A
}
wAKP.avg <- function(x, tr = .2) {
#
# Have J dependent groups. For each pair of groups, compute
# AKP type measure of effect size and average the results
#
# For tr=0, get Cohen d type measure
#
#
a <- wmcpAKP(x, tr = tr)
data.frame(
"Effect.Size" = mean(a[, 3]),
"ci.low" = mean(a[, 4]),
"ci.up" = mean(a[, 5])
)
}
before <- c(190, 210, 300, 240, 280, 170, 280, 250, 240, 220)
now <- c(170, 280, 250, 240, 190, 260, 180, 200, 100, 200)
after <- c(210, 210, 340, 190, 260, 180, 200, 220, 230, 200)
df <- data.frame(before, now, after)
wAKP.avg(df)
#> Effect.Size ci.low ci.up
#> 1 0.2867495 -0.7193793 1.424093
Created on 2021-01-14 by the reprex package (v0.3.0.9001)
from statsexpressions.
Related Issues (20)
- Guidance for using statsExpressions to customize plots?
- In case of `NA`, show it in the label instead of a blank space
- revdep issues with insight HOT 1
- Use `safe_model_parameters()` in `tidy_model_parameters()` instead?
- Avoid duplicating examples
- Don't include CIs in `centrality_description()` by default
- Exact p-value calculation using Wilcoxon rank sum test HOT 1
- DRY test dataset
- Stop converting to factors
- Revep check failure for {insight} HOT 4
- i wonder if `pairwise_comparisons` can do Tukey's HSD? HOT 1
- Use `{roxyglobals}` HOT 1
- Rename `stats_type_switch()` to `extract_stats_type()`
- wilcoxon test default parameters modification HOT 1
- Fix failures on R-devel
- Necessity of patch-level version bumps in dependencies HOT 1
- Introduce parameters to pass extra arguments
- ggbarstats does not show subtitle when paired = TRUE (for more than two groups) HOT 4
- _performance_ revdep-check issues HOT 2
- Switch to MIT licence
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from statsexpressions.