popower <- function(p, odds.ratio, n, n1, n2, alpha=.05)
{
if(missing(n))
n <- n1 + n2
else {
n1 <- n2 <- n / 2
}
p <- p[! is.na(p)]
if(abs(sum(p) - 1) > .00001)
stop('probabilities in p do not add up to 1')
z <- qnorm(1 - alpha / 2)
A <- n2 / n1
ps <- 1 - sum(p ^ 3)
V <- n1 * n2 * n / 3 / ((n + 1) ^ 2) * ps
power <- pnorm(abs(logb(odds.ratio)) * sqrt(V) - z)
eff <- ps / (1 - 1 / n / n)
structure(list(power=power, eff=eff, approx.se=1./sqrt(V)),
class='popower')
}
print.popower <- function(x, ...)
{
cat('Power:',round(x$power,3),
'\nEfficiency of design compared with continuous response:',
round(x$eff, 3),
'\nApproximate standard error of log odds ratio:',
round(x$approx.se, 4), '\n\n')
invisible()
}
posamsize <- function(p, odds.ratio, fraction=.5,
alpha=.05, power=.8)
{
p <- p[!is.na(p)]
if(abs(sum(p) - 1) > .00001)
stop('probabilities in p do not add up to 1')
A <- (1 - fraction) / fraction
log.or <- logb(odds.ratio)
z.alpha <- qnorm(1 - alpha / 2)
z.beta <- qnorm(power)
ps <- 1 - sum(p ^ 3)
n <- 3 * ((A + 1) ^ 2) * (z.alpha + z.beta) ^ 2 / A / (log.or ^ 2) / ps
eff <- ps / (1 - 1 / n / n)
structure(list(n=n, eff=eff), class='posamsize')
}
print.posamsize <- function(x, ...)
{
cat('Total sample size:',round(x$n, 1),
'\nEfficiency of design compared with continuous response:',
round(x$eff, 3),'\n\n')
invisible()
}
pomodm <- function(x=NULL, p, odds.ratio=1) {
if(length(x) && (length(x) != length(p)))
stop('p and x must have same length')
if(length(x) && any(diff(x) <= 0)) stop('x is not sorted or has duplicates')
if(abs(sum(p) - 1) > .00001)
stop('probabilities do not sum to 1')
## Compute cumulative probabilities (exceedances)
cp <- function(p) c(1, 1 - cumsum(p)[- length(p)])
## Compute individual probabilities given exceedance probabilities
ip <- function(ep) {
p <- c(-diff(ep), ep[length(ep)])
if(abs(sum(p) - 1.) > 1e-7) stop('logic error')
p
}
## Function to alter a vector of individual probabilities by a given
## odds ratio on the exceedance probabilities
pmod <- function(p, or) {
## Compute exceedence probabilities
ep <- cp(p)
## Apply odds ratio
ep <- plogis(qlogis(ep) + log(or))
## Get back to individual probabilities
ip(ep)
}
## Convert individual probabilities under the odds ratio
p <- pmod(p, odds.ratio)
if(! length(x)) return(p)
## Compute mean and weighted median
xbar <- sum(p * x)
xmed1 <- approx(cumsum(p), x, xout=0.5)$y
xmed2 <- approx(rev(cumsum(rev(p))), x, xout=0.5)$y
xmed <- (xmed1 + xmed2) / 2
c(mean=xbar, median=xmed)
}
simPOcuts <- function(n, nsim=10, odds.ratio=1, p) {
if(abs(sum(p) - 1.) > 1e-5)
stop('probabilities in p must sum to 1')
p0 <- p
p1 <- pomodm(p=p0, odds.ratio=odds.ratio)
lp <- length(p)
yval <- if(length(names(p))) names(p) else as.character(1 : lp)
or <- matrix(NA, nrow=nsim, ncol=lp - 1,
dimnames=list(paste('Simulation', 1 : nsim),
paste0('y>=', yval[-1])))
for(i in 1 : nsim) {
y0 <- sample(1 : lp, n / 2, prob=p0, replace=TRUE)
y1 <- sample(1 : lp, n / 2, prob=p1, replace=TRUE)
for(ycut in 2 : lp){
prop0 <- mean(y0 >= ycut)
prop1 <- mean(y1 >= ycut)
or[i, ycut - 1] <- (prop1 / (1. - prop1)) / (prop0 / (1. - prop0))
}
}
or
}
simRegOrd <- function(n, nsim=1000, delta=0, odds.ratio=1, sigma,
p=NULL, x=NULL, X=x, Eyx, alpha=0.05, pr=FALSE) {
if (!requireNamespace("rms", quietly = TRUE))
stop("This function requires the 'rms' package.")
if(length(x) && (length(x) != n))
stop('x must be omitted or have length n')
if((n %% 2) != 0) stop('n must be an even integer')
treat <- c(rep(0, n / 2), rep(1, n / 2))
X <- cbind(X, treat)
betas <- se <- rep(NA, nsim)
if(length(p)) {
p1 <- pomodm(p=p, odds.ratio=odds.ratio)
yordinal <- 0 : (length(p) - 1)
}
xb <- delta * treat + (if(length(x)) Eyx(x) else 0)
for(i in 1 : nsim) {
if(pr) cat('Iteration', i, '\r')
y <- xb + rnorm(n, mean=0, sd=sigma)
if(length(p)) {
## Override y with length(p) - 1 ordinal categories
yo <- ifelse(treat == 0,
sample(yordinal, n, prob=p, replace=TRUE),
sample(yordinal, n, prob=p1, replace=TRUE))
ymax <- max(y)
y <- ifelse(yo == 0, y, ceiling(ymax) + yo)
}
f <- rms::orm.fit(X, y, maxit=20)
if(! f$fail) {
betas[i] <- coef(f)[length(coef(f))] ## coef of treatment
k <- nrow(f$var)
se[i] <- sqrt(f$var[k, k])
}
}
if(pr) cat('\n')
pvals <- 1 - pchisq((betas / se) ^ 2, 1)
power <- mean(pvals < alpha, na.rm=TRUE)
list(n=n, delta=delta, sigma=sigma, power=power, betas=betas, se=se,
pvals=pvals)
}
propsPO <- function(formula, odds.ratio=NULL, ref=NULL, data=NULL,
ncol=NULL, nrow=NULL) {
v <- all.vars(formula)
d <- model.frame(formula, data=data)
y <- d[[v[1]]]
yl <- label(y, default=v[1])
y <- as.factor(y)
d[[v[1]]] <- y
x <- d[[v[2]]]
xl <- label(x, default=v[2])
s <- sn <- NULL
sl <- NULL
if(length(v) > 2) {
if(length(odds.ratio))
stop('odds ratio may not be specified when a stratification variable is included')
s <- d[[v[3]]]
sl <- label(s, default=v[3])
s <- as.factor(s)
sn <- 's'
}
names(d) <- c('y', 'x', sn)
## ggplot2 bar chart puts first category at the top
## Let's put it at the bottom
revo <- function(z) {
z <- as.factor(z)
factor(z, levels=rev(levels(as.factor(z))))
}
# For each x compute the vector of proportions of y categories
# Assume levels are in order
# Put numerators and denominators into a character string so that
# data.table [ operator can operate on one variable
# The delimiter is a single space
g <- function(y) {
tab <- table(y)
structure(paste(tab, rep(sum(tab), length(tab))), names=names(tab))
}
atxt <- function(d, strat=NULL, or=FALSE) {
if(! or) {
z <- d$prop
num <- as.numeric(sub(' .*', '', z)) # get first number
denom <- as.numeric(sub('.* ', '', z)) # get second number
d$prop <- num / denom
}
d$txt <- paste0(sl, if(length(strat)) ': ', as.character(strat),
if(length(strat)) '<br>',
xl, ': ', as.character(d$x), '<br>',
yl, ': ', as.character(d$y), '<br>',
'Proportion: ', round(d$prop, 3),
if(! or) '\u2003',
if(! or) markupSpecs$html$frac(num, denom, size=90))
d
}
d <- data.table(d)
if(! length(s)) { # stratification variable not present
p <- d[, as.list(g(y)), by=x]
pm <- melt(p, id=1, variable.name='y', value.name='prop')
pm <- atxt(pm)
plegend <- guides(fill=guide_legend(title=yl))
if(! length(odds.ratio)) {
gg <- ggplot(pm, aes(x=as.factor(x), y=prop, fill=revo(y), label=txt)) +
geom_col() + plegend + xlab(xl) + ylab('Proportion')
return(gg)
}
} else { # stratification variable present; odds ratio must be absent
p <- d[, as.list(g(y)), by=.(x, s)]
pm <- melt(p, id=c('x', 's'), variable.name='y', value.name='prop')
pm <- atxt(pm, pm$s)
plegend <- guides(fill=guide_legend(title=yl))
gg <- ggplot(pm, aes(x=as.factor(x), y=prop, fill=revo(y), label=txt)) +
facet_wrap(~ s, ncol=ncol, nrow=nrow) +
geom_col() + plegend + xlab(xl) + ylab('Proportion')
return(gg)
}
## odds.ratio is present
if(! length(ref)) ref <- p$x[1]
pfx <- as.matrix(p[x == ref, -1])
propfirstx <- as.numeric(sub(' .*', '', pfx)) /
as.numeric(sub('.* ', '', pfx))
.g. <- function(x) {
w <- pomodm(p=propfirstx, odds.ratio=odds.ratio(x))
names(w) <- levels(y)
w
}
pa <- d[, as.list(.g.(x)), by=x]
pma <- melt(pa, id=1, variable.name='y', value.name='prop')
pma <- atxt(pma, or=TRUE)
w <- rbind(cbind(type=1, pm),
cbind(type=2, pma))
w$type <- factor(w$type, 1 : 2,
c('Observed', 'Asssuming Proportional Odds'))
ggplot(w, aes(x=as.factor(x), y=prop, fill=revo(y), label=txt)) +
geom_col() +
facet_wrap(~ type, nrow=2) +
plegend + xlab(xl) + ylab('Proportion')
}
utils::globalVariables(c('.', 'prop'))
propsTrans <- function(formula, data=NULL, labels=NULL, arrow='\u2794',
maxsize=12, ncol=NULL, nrow=NULL) {
v <- all.vars(formula)
d <- model.frame(formula, data=data)
y <- as.factor(d[[v[1]]])
x <- d[[v[2]]]
xlab <- label(x, default=v[2])
id <- as.character(d[[v[3]]])
uid <- sort(unique(id))
nid <- length(uid)
times <- if(is.factor(x)) levels(x) else sort(unique(x))
nt <- length(times)
itrans <- integer(0)
Prev <- Cur <- Frac <- character(0)
prop <- frq <- numeric(0)
mu <- markupSpecs$html
arrowbr <- paste0(' ', arrow, '<br>')
arrow <- paste0(' ', arrow, ' ')
for(it in 2 : nt) {
prev <- cur <- rep(NA, nid)
names(prev) <- names(cur) <- uid
i <- x == times[it - 1]
j <- x == times[it]
prev[id[i]] <- as.character(y[i])
cur [id[j]] <- as.character(y[j])
tab <- table(prev, cur)
rowf <- rowSums(tab)
tab <- as.data.frame(tab)
# tab <- subset(tab, Freq > 0)
tab$denom <- rowf[tab$prev]
tab$prop <- tab$Freq / tab$denom
frq <- c(frq, tab$Freq)
Prev <- c(Prev, as.character(tab$prev))
Cur <- c(Cur, as.character(tab$cur))
prop <- c(prop, tab$Freq / tab$denom)
Frac <- c(Frac, paste0(round(tab$Freq / tab$denom, 3),
mu$lspace,
mu$frac(tab$Freq, tab$denom, size=90)))
k <- length(tab$Freq)
itrans <- c(itrans, rep(it, k))
}
Prev <- factor(Prev, levels(y))
Cur <- factor(Cur, levels(y))
trans <- factor(itrans, 2 : nt,
labels=paste0(xlab, ' ',
times[1 : (nt - 1)], arrow, times[2 : nt]))
transp <- factor(itrans, 2 : nt,
labels=paste0(xlab, ' ',
times[1 : (nt - 1)], arrow, times[2 : nt]))
w <- data.frame(trans, transp, Prev, Cur, prop, Frac, frq,
txt=if(! length(labels))
ifelse(Prev == Cur,
paste0('Stay at:', as.character(Prev)),
paste0(as.character(Prev),
arrowbr, as.character(Cur)))
else
ifelse(Prev == Cur,
paste0('Stay at:', labels[as.integer(Prev)]),
paste0(labels[as.integer(Prev)],
arrowbr, labels[as.integer(Cur)])))
w$txt <- paste0(w$transp, '<br>', w$txt, '<br>', w$Frac)
w <- subset(w, frq > 0)
ggplot(w, aes(x=Prev, y=Cur, size=prop, label=txt)) +
facet_wrap(~ trans, ncol=ncol, nrow=nrow) +
geom_point() + scale_size(range = c(0, maxsize)) +
xlab('Previous State') + ylab('Current State') +
guides(size = guide_legend(title='Proportion'))
# ggplot(w, aes(x=Prev, y=prop, fill=Cur)) +
# facet_wrap(~ trans, ncol=ncol, nrow=nrow) +
# geom_col() + #scale_size(range = c(0, 12)) +
# xlab('Previous State') + ylab('Proportion') +
# guides(fill = guide_legend(title='Current State'))
}
multEventChart <-
function(formula, data=NULL,
absorb=NULL, sortbylast=FALSE,
colorTitle=label(y), eventTitle='Event',
palette='OrRd',
eventSymbols=c(15, 5, 1:4, 6:10),
timeInc=min(diff(unique(x))/2)) {
v <- all.vars(formula)
d <- model.frame(formula, data=data)
y <- as.factor(d[[v[1]]])
y <- factor(y, levels=rev(levels(y)))
x <- d[[v[2]]]
xlab <- label(x, default=v[2])
id <- as.factor(d[[v[3]]])
# Optionally sort subjects by last status, assuming levels of y
# are in ascending order of badness
if(sortbylast) {
last <- tapply(1 : length(y), id,
function(i) {
times <- x[i]
status <- y[i]
as.integer(status[which.max(times)])
})
i <- order(last, levels(id), decreasing=TRUE)
id <- factor(id, levels=levels(id)[i])
} else id <- factor(id, levels=rev(levels(id)))
ab <- as.character(y) %in% absorb
nay <- setdiff(levels(y), absorb)
event <- y
event[! ab] <- NA
y[ab] <- NA
de <- subset(data.frame(id, x, y, event), ! is.na(event))
ggplot(mapping=aes(x = id, y = x, fill = y)) +
scale_fill_brewer(colorTitle, palette=palette, direction=-1, breaks=nay) +
geom_segment(aes(x = id, xend = id, y = 0, yend = x - timeInc),
lty = 3) +
geom_tile(width = timeInc) +
scale_y_continuous(breaks=min(x) : max(x)) +
geom_point(aes(x = id, y = x - timeInc, shape = event), data=de) +
scale_shape_manual(eventTitle, values=eventSymbols[1 : length(absorb)],
labels=c(absorb)) +
guides(fill=guide_legend(override.aes=list(shape=NA), order=1)) +
coord_flip() + labs(y=xlab, x='')
}
utils::globalVariables('txt')
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.