Miscellaneous S Examples
Graphical Demo of Numerical Differentiation
lim <- function(f, low, high) {
xs <- seq(low, high, length=500)
plot(xs, f(xs), type='l', xlab='x', ylab='f(x)')
w <- (high-low)/2
cat('\nClick mouse to select point at which to evaluate derivative\n')
x <- locator(1)$x
colors <- rep(1:10,length=30)
for(i in 1:30) {
x1 <- max(x - w, low)
x2 <- min(x + w, high)
slope <- (f(x2)-f(x1))/(x2-x1)
intercept <- f(x) - slope*x
# lines(c(x1,x2), intercept + slope*c(x1,x2), col=i)
abline(intercept, slope, col=colors[i])
cat(format(slope),'\n')
w <- w*.8
}
x
}
lim(sin, -2, 2)
lim(function(x)1/x, .01,20)
lim(tan, -10, 10)
lim(sqrt, 0, 100)
lim(sin, -10, 10)
Simulate Difference in Life Expectancy for Mixed Exponential/Weibull Distribution
This program plots the difference in life expectancy for various 1y survival probabilities and hazard ratios, assuming exponential survival for 2 treatments for years 0-4 and Weibull survival given survived 4 years, with gamma=2 and same 1y conditional survival as original 1y survival
hr <- seq(.1,.9,by=.1)
gam <- 1.5
i <- 0
pli <- F
if(!pli)plot(0,0,type="n",xlab="1-Year Survival for Control Subject",
ylab="Improvement in Mean Life Length (Years)",
xlim=c(0,1),ylim=c(0,4))
for(h in hr) {
i <- i+1
ps <- seq(.0001,.9999,length=100)
dif <- double(length(ps))
k <- 0
for(p in ps) {
k <- k+1
l1 <- -log(p)
l2 <- l1
times <- seq(0, 25, length=500)
S1 <- exp(-ifelse(times<=4, l1*times, l1*4+l2*(times-4)^gam))
S2 <- exp(-ifelse(times<=4, h*l1*times, h*l1*4+l2*(times-4)^gam))
if(pli) {
plot(times, S1, type="l", xlab="Years",ylab="Survival Probability")
lines(times, S2, lty=2)
legend(locator(1),c("Control","Treatment"),lty=1:2,bty="n")
title("Survival for 4-year Treatment Effect, Hazard Ratio=0.5")
}
mu1 <- sum(S1)*diff(times[1:2])
mu2 <- sum(S2)*diff(times[1:2])
dif[k] <- mu2-mu1
}
if(!pli) {
lines(ps, dif, lty=i)
maxd <- max(dif)
smax <- ps[dif==maxd]
text(smax, maxd+.1, format(h), cex=.9) #was cex=.6 for small plot
}
}
if(!pli) title("Increase in Mean Life Length for Various Hazard Ratios")
Blocked Randomization with Random Block Sizes
b <- 20 # number of blocks
s <- c(2,4) # number of subjects per block
set.seed(1) # so can reproduce results
m <- sample(s, b, replace=T) # actual number of subjects over blocks
# Get treatment assignments for each block, randomized within blocks
treats <- vector('list', b)
for(i in 1:b) treats[[i]] <-
sample(c(rep('a',m[i]/2),rep('b',m[i]/2)))
treats # print for statistician use
unlist(treats) # print blinded to block boundaries