Nothing
## ----chunkname, echo=-1-------------------------------------------------------
data.table::setDTthreads(2)
## ----echo = FALSE, message = FALSE--------------------------------------------
library(simstudy)
library(ggplot2)
library(scales)
library(grid)
library(gridExtra)
library(survival)
library(gee)
library(data.table)
library(ordinal)
odds <- function (p) p/(1 - p) # TODO temporary remove when added to package
plotcolors <- c("#B84226", "#1B8445", "#1C5974")
cbbPalette <- c("#B84226","#B88F26", "#A5B435", "#1B8446",
"#B87326","#B8A526", "#6CA723", "#1C5974")
ggtheme <- function(panelback = "white") {
ggplot2::theme(
panel.background = element_rect(fill = panelback),
panel.grid = element_blank(),
axis.ticks = element_line(colour = "black"),
panel.spacing =unit(0.25, "lines"), # requires package grid
panel.border = element_rect(fill = NA, colour="gray90"),
plot.title = element_text(size = 8,vjust=.5,hjust=0),
axis.text = element_text(size=8),
axis.title = element_text(size = 8)
)
}
## ----options, echo = FALSE----------------------------------------------------
options(digits = 2)
## ----threshold, fig.width = 5.25, fig.height = 3.5, echo = FALSE--------------
# preliminary libraries and plotting defaults
library(ggplot2)
library(data.table)
my_theme <- function() {
theme(panel.background = element_rect(fill = "grey90"),
panel.grid = element_blank(),
axis.ticks = element_line(colour = "black"),
panel.spacing = unit(0.25, "lines"),
plot.title = element_text(size = 12, vjust = 0.5, hjust = 0),
panel.border = element_rect(fill = NA, colour = "gray90"))
}
# create data points density curve
x <- seq(-6, 6, length = 1000)
pdf <- dlogis(x, location = 0, scale = 1)
dt <- data.table(x, pdf)
# set thresholds for Group A
# P = c(0.31, 0.29, .20, 0.20)
thresholdA <- c(-0.8, 0.4, 1.4)
pdf <- dlogis(thresholdA)
grpA <- data.table(threshold = thresholdA, pdf)
aBreaks <- c(-6, grpA$threshold, 6)
# plot density with cut points
dt[, grpA := cut(x, breaks = aBreaks, labels = F, include.lowest = TRUE)]
p1 <- ggplot(data = dt, aes(x = x, y = pdf)) +
geom_line() +
geom_area(aes(x = x, y = pdf, group = grpA, fill = factor(grpA))) +
geom_hline(yintercept = 0, color = "grey50") +
annotate("text", x = -5, y = .28, label = "unexposed", size = 5) +
scale_fill_manual(values = c("#d0d7d1", "#bbc5bc", "#a6b3a7", "#91a192", "#7c8f7d"),
labels = c("strongly disagree", "disagree", "agree", "strongly agree"),
name = "Frequency") +
scale_x_continuous(breaks = thresholdA) +
scale_y_continuous(limits = c(0, 0.3), name = "Density") +
my_theme() +
theme(legend.position = c(.85, .7),
legend.background = element_rect(fill = "grey90"),
legend.key = element_rect(color = "grey90"))
p1
## ----plotB, fig.width = 5.25, fig.height = 3.5, echo = FALSE------------------
pA= plogis(c(thresholdA, Inf)) - plogis(c(-Inf, thresholdA))
probs <- data.frame(pA)
rownames(probs) <- c("P(Resp = 1)", "P(Resp = 2)",
"P(Resp = 3)", "P(Resp = 4)")
probA <- data.frame(
cprob = plogis(thresholdA),
codds = plogis(thresholdA)/(1-plogis(thresholdA)),
lcodds = log(plogis(thresholdA)/(1-plogis(thresholdA)))
)
rownames(probA) <- c("P(Grp < 2)", "P(Grp < 3)", "P(Grp < 4)")
thresholdB <- thresholdA + 0.7
pdf <- dlogis(thresholdB)
grpB <- data.table(threshold = thresholdB, pdf)
bBreaks <- c(-6, grpB$threshold, 6)
pB = plogis(c(thresholdB, Inf)) - plogis(c(-Inf, thresholdB))
probs <- data.frame(pA, pB)
rownames(probs) <- c("P(Resp = 1)", "P(Resp = 2)",
"P(Resp = 3)", "P(Resp = 4)")
# Plot density for group B
dt[, grpB := cut(x, breaks = bBreaks, labels = F, include.lowest = TRUE)]
p2 <- ggplot(data = dt, aes(x = x, y = pdf)) +
geom_line() +
geom_area(aes(x = x, y = pdf, group = grpB, fill = factor(grpB))) +
geom_hline(yintercept = 0, color = "grey5") +
geom_segment(data=grpA,
aes(x=threshold, xend = threshold, y=0, yend=pdf),
size = 0.3, lty = 2, color = "#857284") +
annotate("text", x = -5, y = .28, label = "exposed", size = 5) +
scale_fill_manual(values = c("#d0d7d1", "#bbc5bc", "#a6b3a7", "#91a192", "#7c8f7d"),
name = "Frequency") +
scale_x_continuous(breaks = thresholdB) +
scale_y_continuous(limits = c(0.0, 0.3), name = "Density") +
my_theme() +
theme(legend.position = "none")
p2
## ----acuts--------------------------------------------------------------------
baseprobs <- c(0.31, 0.29, .20, 0.20)
defA <- defData(varname = "exposed", formula = "1;1", dist = "trtAssign")
defA <- defData(defA, varname = "z", formula = "-0.7*exposed", dist = "nonrandom")
set.seed(130)
dT_1_cat <- genData(25000, defA)
dX <- genOrdCat(dT_1_cat, adjVar = "z", baseprobs, catVar = "r")
## ----ordinal------------------------------------------------------------------
library(ordinal)
clmFit <- clm(r ~ exposed, data = dX)
summary(clmFit)
## -----------------------------------------------------------------------------
(logOdds.unexp <- log(odds(cumsum(dX[exposed == 0, prop.table(table(r))])))[1:3])
## -----------------------------------------------------------------------------
(logOdds.expos <- log(odds(cumsum(dX[exposed == 1, prop.table(table(r))])))[1:3])
## -----------------------------------------------------------------------------
logOdds.expos - logOdds.unexp
## ----echo=FALSE, fig.width=6, fig.height=3.5----------------------------------
getCumProbs <- function(coefs) {
cumprob0 <- data.table(
cumprob = c(1/(1 + exp(-coefs[which(rownames(coefs) != "exposed")])), 1),
r = factor(1 : nrow(coefs)),
exposed = 0
)
cumprob1 <- data.table(
cumprob = c(1/(1 + exp(-coefs[which(rownames(coefs) != "exposed")] +
coefs["exposed", 1])), 1),
r = factor(1 : nrow(coefs)),
exposed = 1
)
rbind(cumprob0, cumprob1)[]
}
fitPlot <- function(dx) {
clmFit <- clm(r ~ exposed, data = dx)
coefs <- coef(summary(clmFit))
cumModProbs <- getCumProbs(coefs)
cumObsProbs <- dx[, .N, keyby = .(exposed, r)]
cumObsProbs[, cumprob := cumsum(N)/sum(N) , keyby = exposed]
ggplot(data = cumObsProbs, aes(x = r, y = cumprob, color = factor(exposed))) +
geom_line(data = cumModProbs, alpha = 1, aes(group=exposed)) +
geom_point(size = 1.25) +
ylab("cumulative probability") +
xlab("ordinal category") +
theme(panel.grid = element_blank(),
legend.title = element_blank()) +
scale_color_manual(values = c("#7c8e8f", "#8f7c8e"), labels = c("Not exposed", "Exposed"))
}
fitPlot(dX)
## ----plotNP, fig.width = 5.25, fig.height = 3.5, echo = FALSE-----------------
pA <- plogis(c(thresholdA, Inf)) - plogis(c(-Inf, thresholdA))
probs <- data.frame(pA)
rownames(probs) <- c("P(Resp = 1)", "P(Resp = 2)",
"P(Resp = 3)", "P(Resp = 4)")
probA <- data.frame(
cprob = plogis(thresholdA),
codds = plogis(thresholdA)/(1-plogis(thresholdA)),
lcodds = log(plogis(thresholdA)/(1-plogis(thresholdA)))
)
rownames(probA) <- c("P(Grp < 2)", "P(Grp < 3)", "P(Grp < 4)")
thresholdB <- thresholdA + c(0.4, 0.0, 1.7)
pdf <- dlogis(thresholdB)
grpB <- data.table(threshold = thresholdB, pdf)
bBreaks <- c(-6, grpB$threshold, 6)
pB = plogis(c(thresholdB, Inf)) - plogis(c(-Inf, thresholdB))
probs <- data.frame(pA, pB)
rownames(probs) <- c("P(Resp = 1)", "P(Resp = 2)",
"P(Resp = 3)", "P(Resp = 4)")
# Plot density for group B
dt[, grpB := cut(x, breaks = bBreaks, labels = F, include.lowest = TRUE)]
p3 <- ggplot(data = dt, aes(x = x, y = pdf)) +
geom_line() +
geom_area(aes(x = x, y = pdf, group = grpB, fill = factor(grpB))) +
geom_hline(yintercept = 0, color = "grey5") +
geom_segment(data=grpA,
aes(x=threshold, xend = threshold, y=0, yend=pdf),
size = 0.3, lty = 2, color = "#857284") +
annotate("text", x = -4, y = .28, label = "non-proportional", size = 5) +
scale_fill_manual(values = c("#d0d7d1", "#bbc5bc", "#a6b3a7", "#91a192", "#7c8f7d"),
name = "Frequency") +
scale_x_continuous(breaks = thresholdB) +
scale_y_continuous(limits = c(0.0, 0.3), name = "Density") +
my_theme() +
theme(legend.position = "none")
p3
## -----------------------------------------------------------------------------
baseprobs <- c(0.31, 0.29, .20, 0.20)
npAdj <- c(-0.4, 0.0, -1.7, 0)
dX <- genOrdCat(dT_1_cat, baseprobs = baseprobs, adjVar = "z",
catVar = "r", npVar = "exposed", npAdj = npAdj)
## -----------------------------------------------------------------------------
(logOdds.unexp <- log(odds(cumsum(dX[exposed == 0, prop.table(table(r))])))[1:3])
## -----------------------------------------------------------------------------
(logOdds.expos <- log(odds(cumsum(dX[exposed == 1, prop.table(table(r))])))[1:3])
## -----------------------------------------------------------------------------
logOdds.expos - logOdds.unexp
## ----echo=FALSE, fig.width=6, fig.height=3.5, warning=FALSE-------------------
fitPlot(dX)
## -----------------------------------------------------------------------------
baseprobs <- matrix(c(0.2, 0.1, 0.7,
0.7, 0.2, 0.1,
0.5, 0.2, 0.3,
0.4, 0.2, 0.4,
0.6, 0.2, 0.2),
nrow = 5, byrow = TRUE)
# generate the data
set.seed(333)
dT_5_cat <- genData(10000)
dX <- genOrdCat(dT_5_cat, adjVar = NULL, baseprobs = baseprobs,
prefix = "q", rho = 0.15, corstr = "cs", asFactor = FALSE)
## -----------------------------------------------------------------------------
round(dX[, cor(cbind(q1, q2, q3, q4, q5))], 2)
## -----------------------------------------------------------------------------
dM <- melt(dX, id.vars = "id")
dProp <- dM[ , prop.table(table(value)), by = variable]
dProp[, response := rep(seq(3), 5)]
# observed probabilities
dcast(dProp, variable ~ response, value.var = "V1", fill = 0)
# specified probabilities
baseprobs
## -----------------------------------------------------------------------------
dX <- genOrdCat(dT_5_cat, adjVar = NULL, baseprobs = baseprobs,
prefix = "q", rho = 0.40, corstr = "ar1", asFactor = FALSE)
# correlation
round(dX[, cor(cbind(q1, q2, q3, q4, q5))], 2)
dM <- melt(dX, id.vars = "id")
dProp <- dM[ , prop.table(table(value)), by = variable]
dProp[, response := rep(seq(3), 5)]
# probabilities
dcast(dProp, variable ~ response, value.var = "V1", fill = 0)
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.