Miscellaneous S Examples

• Graphical Demo of Numerical Differentiation
• Simulate Difference in Life Expectancy for Mixed Exponential/Weibull Distribution
• Blocked Randomization with Random Block Sizes

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