You are here:
Vanderbilt Biostatistics Wiki
>
Main Web
>
StatComp
>
RS
>
SexamplesMisc
(11 Jul 2004,
FrankHarrell
)
(raw view)
E
dit
A
ttach
---+ Miscellaneous S Examples * [[#NumderiveLimits][Graphical Demo of Numerical Differentiation]] * [[#SimdifLifeexpect][Simulate Difference in Life Expectancy for Mixed Exponential/Weibull Distribution]] * [[#RandomBlocks][Blocked Randomization with Random Block Sizes]] #NumderiveLimits ---+++ Graphical Demo of Numerical Differentiation <verbatim> 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) </verbatim> #SimdifLifeexpect ---+++ 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 <verbatim> 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") </verbatim> #RandomBlocks ---+++ Blocked Randomization with Random Block Sizes <verbatim> 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 </verbatim>
E
dit
|
A
ttach
|
P
rint version
|
H
istory
: r1
|
B
acklinks
|
V
iew topic
|
Edit
w
iki text
|
M
ore topic actions
Topic revision: r1 - 11 Jul 2004,
FrankHarrell
Main
Department Home Page
Biostatistics Graduate Program
Vanderbilt University Medical Center
Main Web
Main Web Home
Search
Recent Changes
Changes
Topic list
Biostatistics Webs
Archive
Main
Sandbox
System
Register
|
Log In
Copyright © 2013-2022 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