# 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
```
Topic revision: r1 - 11 Jul 2004, FrankHarrell

• Biostatistics Webs Copyright © 2013-2020 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding Vanderbilt Biostatistics Wiki? Send feedback