require(survival)

# Create the simplest test data set 
test1 <- list(time=c(4,3,1,1,2,2,3), 
              status=c(1,1,1,0,1,1,0), 
              x=c(0,2,1,1,1,0,0), 
              sex=c(0,0,0,0,1,1,1)) 
# Fit a stratified model 
coxph(Surv(time, status) ~ x + strata(sex), test1) 

data <- data.frame(time=c(4,3,1,1,2,2,3), 
              status=c(1,1,1,0,1,1,0), 
              x=c(0,2,1,1,1,0,0), 
              sex=c(0,0,0,0,1,1,1)) 
coxph(Surv(time, status) ~ x + strata(sex), data) 

# Bootstrap of beta
index <- 1:nrow(data)
B <- 200
beta.boot <- rep(0,B)
for(b in 1:B){
    i.b <- sample(index,nrow(data),rep=TRUE)
    beta.boot[b] <- coxph(Surv(time, status) ~ x, data[i.b,])$coef
}
sd(beta.boot)
    

# Simulate data
x <- rep(c(0,1),c(150,150))
lambda <- exp(0.8*x)
y <- rexp(300,lambda)
cen <- runif(300,0,3)
data <- data.frame(time=pmin(y,cen),status=as.numeric(cen >= y),x=x)

coxph(Surv(time, status) ~ x, data)

# Bootstrap of beta
index <- 1:nrow(data)
B <- 200
beta.boot <- rep(0,B)
for(b in 1:B){
    i.b <- sample(index,nrow(data),rep=TRUE)
    beta.boot[b] <- coxph(Surv(time, status) ~ x, data[i.b,])$coef
}
sd(beta.boot)

# PPR
data(rock)
area1 <- rock$area/10000; peri1 <- rock$peri/10000
rock.ppr <- ppr(log(perm) ~ area1 + peri1 + shape,
                data = rock, nterms = 2)
rock.ppr

# bootstrap ppr
rock.data <- rock
rock.data <- cbind(rock.data,area1,peri1)
index <- 1:nrow(rock.data)
B <- 200
beta.boot <- rep(0,B)
for(b in 1:B){
    i.b <- sample(index,nrow(rock.data),rep=TRUE)
    rock.ppr.b <- ppr(log(perm) ~ area1 + peri1 + shape,
                data = rock.data[i.b,], nterms = 2)
    beta.boot[b] <- rock.ppr.b$beta[2]
}
sd(beta.boot)