Nothing
## ----echo = FALSE, warning=FALSE, message = FALSE, results = 'hide'-----------
cat("this will be hidden; use for general initializations.\n")
library(superb)
library(ggplot2)
## -----------------------------------------------------------------------------
grp1 <- c( 56, 54, 73, 46, 59, 62, 55, 53, 77, 60, 69, 66, 63,
62, 53, 82, 74, 70, 65, 70, 72, 65, 56, 58, 83)
grp2 <- c( 51, 99, 194, 123, 40, 83, 87, 117, 46, 89, 61, 81, 53,
141, 52, 53, 39, 96, 14, 81, 63, 66, 80, 113, 82)
## -----------------------------------------------------------------------------
dtaHetero <- data.frame( cbind(
id = 1:(length(grp1)+length(grp2)),
group = c(rep(1,length(grp1)), rep(2, length(grp2)) ),
score = c(grp1, grp2)
))
head(dtaHetero)
## -----------------------------------------------------------------------------
library(superb) # to make the summary plot
library(ggplot2) # for all the graphic directives
library(gridExtra) # for grid.arrange
ornate = list(
xlab("Difference"),
scale_x_discrete(labels=c("Pre-\nTreatment","Post-\nTreatment")),
ylab("Score"),
coord_cartesian( ylim = c(40,+110) ),
theme_light(base_size = 10) +
theme( plot.subtitle = element_text(size=16))
)
## ----fig.height=4, fig.width=3, fig.cap = "**Figure 1**. Plot of dtaHerero showing heterogeneous error bars."----
pt <- superbPlot(dtaHetero,
BSFactors = "group",
variables = "score",
adjustments = list(purpose = "difference"),
gamma = 0.95,
statistic = "mean",
errorbar = "CI",
plotStyle = "line",
lineParams = list(alpha = 0) #the line is made transparent
) + ornate
pt
## -----------------------------------------------------------------------------
# Welch's rectified degrees of freedom
wdf <- WelchDegreeOfFreedom(dtaHetero, "score","group")
## ----fig.height=4, fig.width=3, fig.cap = "**Figure 2**. Plot of dtaHerero with rectified degree of freedom."----
pw <- superbPlot(dtaHetero,
BSFactors = "group",
variables = "score",
adjustments = list(purpose = "difference"),
gamma = c(0.95, wdf), # new!
statistic = "mean",
errorbar = "CIwithDF", # new!
plotStyle = "halfwidthline",
lineParams = list(alpha = 0)
) + ornate
pw
## ----fig.height=4, fig.width=3, fig.cap = "**Figure 3**. Plot of dtaHerero with rectified degree of freedom and Tryon' difference-adjusted error bars."----
pwt <- superbPlot(dtaHetero,
BSFactors = "group",
variables = "score",
adjustments = list(purpose = "tryon"), #new!
gamma = c(0.95, wdf),
statistic = "mean",
errorbar = "CIwithDF",
plotStyle = "halfwidthline",
lineParams = list(alpha = 0)
)+ ornate
pwt
## -----------------------------------------------------------------------------
# get the summary statistics with superbData
t <- superbData(dtaHetero,
BSFactors = "group",
variables = "score",
adjustments = list(purpose = "difference"),
gamma = 0.95,
statistic = "mean", errorbar = "CI"
)
# keep only the summary statistics:
t2 <- t$summaryStatistics
# the length is in column "upperwidth", for lines 1 and 2,
# so lets do the mean in the square sense:
tmean2 <- sqrt( (t2$upperwidth[1]^2 + t2$upperwidth[2]^2)/2 )
## -----------------------------------------------------------------------------
tmean <- (t2$upperwidth[1] + t2$upperwidth[2] )/2
## -----------------------------------------------------------------------------
# get the summary statistics with superbData
wt <- superbData(dtaHetero,
BSFactors = "group",
variables = "score",
adjustments = list(purpose = "tryon"),
gamma = c(0.95, wdf),
statistic = "mean", errorbar = "CIwithDF"
)
wt2 <- wt$summaryStatistics
wmean <- (wt2$upperwidth[1]+wt2$upperwidth[2]) / 2
## ----fig.height=4, fig.width=9, fig.cap = "**Figure 4**. All three plots with relevant markers in red."----
# showing all three plots, with reference lines in red
grid.arrange(
pt + labs(subtitle="Difference-adjusted 95% CI\n with default degree of freedom") +
geom_text( x = 1.15, y = mean(grp1)+(mean(grp2)-mean(grp1))/2, label = "power-2 mean", angle = 90) +
geom_text( x = 1.55, y = mean(grp1)+(mean(grp2)-mean(grp1))/2, label = "regular mean", angle = 90) +
geom_hline(yintercept = mean(grp1)+(mean(grp2)-mean(grp1))/2-tmean2/2, colour = "red", linewidth = 0.5, linetype=2) +
geom_hline(yintercept = mean(grp1)+(mean(grp2)-mean(grp1))/2+tmean2/2, colour = "red", linewidth = 0.5, linetype=2) +
# arrow for the power-2 mean
geom_segment(arrow = arrow(length =unit(0.4,"cm")), x=1.33, y=mean(grp1)+(mean(grp2)-mean(grp1))/2,
xend=1.33, yend=mean(grp1)+(mean(grp2)-mean(grp1))/2+tmean2/2) +
geom_segment(arrow = arrow(length =unit(0.4,"cm")), x=1.333, y=mean(grp1)+(mean(grp2)-mean(grp1))/2,
xend=1.33, yend=mean(grp1)+(mean(grp2)-mean(grp1))/2-tmean2/2) +
# arrow for the regular mean
geom_segment(arrow = arrow(length =unit(0.4,"cm")), x=1.66, y=mean(grp1)+(mean(grp2)-mean(grp1))/2,
xend=1.66, yend=mean(grp1)+(mean(grp2)-mean(grp1))/2+tmean/2) +
geom_segment(arrow = arrow(length =unit(0.4,"cm")), x=1.66, y=mean(grp1)+(mean(grp2)-mean(grp1))/2,
xend=1.66, yend=mean(grp1)+(mean(grp2)-mean(grp1))/2-tmean/2),
pw + labs(subtitle="Difference-adjusted 95% CI\nwith df from Welch"),
pwt + labs(subtitle="Tryon-adjusted 95% CI\nwith df from Welch") +
geom_hline(yintercept = mean(grp1)+wt2$upperwidth[1]/2, colour = "red", linewidth = 0.5, linetype=2) +
geom_hline(yintercept = mean(grp2)+wt2$lowerwidth[2]/2, colour = "red", linewidth = 0.5, linetype=2),
ncol=3)
## -----------------------------------------------------------------------------
t.test( grp1, grp2,
var.equal=FALSE)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.