.Rhistory

h = exp( b)
lengths <- c(breaks, maxDuration) - c(1, breaks)
H = cumsum( lengths*h)
S = exp(-H)
output <- list(S, PWC)
return(output)
}
plotPWC <- function(breaks, S, maxDuration, subset){
xgrid <- c(1,breaks, maxDuration) #returns the cut points
ygrid <- S #returns the corresponding hazards
graph <- data.frame(x = seq(1,  maxDuration, by = 0.1))
for(i in 1:length(graph$x)){
for(j in 2:length(xgrid)){
if(graph$x[i] < xgrid[j] & graph$x[i] >= xgrid[j-1]){
graph$y[i] <- ygrid[[j-1]]
}
}
}
p <- ggplot(graph[ graph$x < subset, ]) +
geom_line(aes(x = x, y = y)) +
xlab("duration") +
ylab("survival")
return(p)
}
breaks1 <- c(2, 5, 10, 20, 50, 500)
breaks2 <- c(2, 3,4,5, 8, 16, 32, 64, 128, 256, 512)
controls <- c("gender", "urban", 'public',"catholic", "protest", "special")
S1 <- PWC(breaks1, df)[[1]]
plotPWC(breaks1, S1, maxDuration, subset = 50)
# Setup code - loading the data or packages
rm(list=ls())
variables <- c("gender", "marstat", "contract", "schooltype", "urban", "birthyr", "classize")
# packages ---------------------------------------------------------------------
library(stargazer)
library(survival)
library(haven)
#library(xtable)
#library(RCT)
library(foreign)
library(tidyverse)
library(ggfortify)
library(survminer)
library(splines2)
library(casebase)
# for plodurationing the smooth hazard rate
library(visreg)
# for computing piecewise constant hazard rates
library(pch)
#for unobserved heterogeneity models
library(parfm)
# for computing the cumulative hazards
library(zoo)
# reading in data --------------------------------------------------------------
dat <- "C:\\Users\\kiera\\OneDrive\\Documents\\Tinbergen MPhil\\Applied Microeconometrics\\Applied_micro"
setwd(dat)
df <- as.data.frame(read_dta("FlowSpells.dta"), stringsAsFactors=F)
df <- df[df[["sptype"]] == 2, ] #only focus on the sickness periods
# Hazard and survival functions ------------------------------------------------
hazard <- function(df, duration, censored, t){
numerator <- sum((1-df[[censored]])*ifelse(t < df[[duration]] & df[[duration]] <= t + 1, 1, 0))
denominator <- ifelse(any(df[[duration]] > t), sum(df[[duration]] > t), 1) #can't divide by zero
return(numerator/denominator)
}
survival <- function(survivalTable, t){
return(ifelse(t > 1, prod(1 - survivalTable$hazard[1:(t-1)]), 1))
}
# Question 1 -------------------------------------------------------------------
# functions ---------------------------------------------------------------------
# functions that compute the hazard and survival functions
hazard <- function(df, duration, censored, t){
# computing the numerator in a vectorised way
numerator <- sum((1-df[[censored]])*ifelse(t < df[[duration]] & df[[duration]] <= t + 1, 1, 0))
denominator <- ifelse(any(df[[duration]] > t), sum(df[[duration]] > t), 1)
return(numerator/denominator)
}
survival <- function(survivalTable, t){
# this uses that the survival function is the complement of the cumulative
# hazard
return(ifelse(t > 1, prod(1 - survivalTable$hazard[1:(t-1)]), 1))
}
# function to compute the survival table from the dataset
survTableGen <- function(df, col = FALSE, grp = FALSE){
duration <- 1:max(df$splength)
survivalTable <- data.frame(duration)
if(col == FALSE){ #if no column specified, select all data
survivalTable$hazard <- sapply(duration, hazard, df = df,
duration = "splength", censored = "rcensor")
survivalTable$survival <- sapply(duration, survival, survivalTable = survivalTable)
} else { #if column specified, only select data from the selected group
subset <- df[df[[col]]==grp,]
survivalTable$hazard <- sapply(duration, hazard, df = subset,
duration = "splength", censored = "rcensor")
survivalTable$survival <- sapply(duration, survival, survivalTable = survivalTable)
survivalTable[[col]] <- grp
}
return(survivalTable)
}
# function that computes hazard rates and survival function by groups defined
# by continuous variables
sur_haz_grp <- function(col, df){
groups <- unique(df[[col]])
out_dfs <- lapply(groups, survTableGen, col = col, df = df)
out_df_merged <- bind_rows(out_dfs)
hazard <- ggplot(out_df_merged, aes(x =duration, y=hazard, colour = as.factor(get(col)))) +
geom_line() +
theme(legend.title=element_blank())
survival <- ggplot(out_df_merged, aes(x =duration, y=survival, colour = as.factor(get(col)))) +
geom_line() +
theme(legend.title=element_blank())
# finally testing for difference between the survival rates
return(list(hazard, survival))
}
# running code ----------------------------------------------------------------
survivalTable <- survTableGen(df)
ggplot(data = survivalTable, mapping = aes(x =duration, y= hazard)) +
geom_line() + geom_smooth()
ggplot(data = survivalTable, mapping = aes(x =duration, y=survival)) +
geom_line()
# now lets do plots for the first two weeks, and the first year
ggplot(data = survivalTable[survivalTable[["duration"]]<=14,], mapping = aes(x =duration, y=hazard)) +
geom_line() + geom_smooth()
ggplot(data = survivalTable[survivalTable[["duration"]]<=14,], mapping = aes(x =duration, y=survival)) +
geom_line()
ggplot(data = survivalTable[survivalTable[["duration"]]<=365,], mapping = aes(x =duration, y=hazard)) +
geom_line() + geom_smooth()
ggplot(data = survivalTable[survivalTable[["duration"]]<=365,], mapping = aes(x =duration, y=survival)) +
geom_line()
# now lets partition by different columns and compute the hazard rates
# partitions: type of school, gender, urbanisation
df[["school_type"]] <- NA
df[["school_type"]][df[["public"]]==1] <- "Public"
df[["school_type"]][df[["catholic"]]==1] <- "Catholic"
df[["school_type"]][df[["protest"]]==1] <- "Protestant"
df[["school_type"]][df[["special"]]==1] <- "Special"
df[["area_type"]] <- NA
df[["area_type"]][df[["urban"]]==1] <- "Urban_1"
df[["area_type"]][df[["urban"]]==2] <- "Urban_2"
df[["area_type"]][df[["urban"]]==3] <- "Urban_3"
df[["area_type"]][df[["urban"]]==4] <- "Urban_4"
df[["area_type"]][df[["urban"]]==5] <- "Urban_5"
df[["gender"]] <- ifelse(df[["gender"]] == 1, "male", "female")
#choose to cut off at 300 because from there the pictures are uninformative
sur_haz_grp("gender", df[df[["splength"]] < 300, ])
sur_haz_grp("area_type", df[df[["splength"]] < 300, ])
sur_haz_grp("school_type", df[!is.na(df[["school_type"]]) & df[["splength"]] < 300, ])
# finally lets test for a difference between the survival rates
survdiff(Surv(splength, rcensor) ~ gender, data=df)
survdiff(Surv(splength, rcensor) ~ urban, data=df)
survdiff(Surv(splength, rcensor) ~ school_type, data=df)
# Question 2 - parametric survival models --------------------------------------
# simple weibull and exponential models
# In order to use the package, we have to make a status indicator
df$status <- ifelse(df$rcensor == 1, 0, 1)
exp_mod <- survreg(Surv(splength, status) ~ public, data=df, dist="exponential")
summary(exp_mod)
weib_mod <- survreg(Surv(splength, status) ~ public, data=df, dist="weibull")
summary(weib_mod)
paste("The scale parameter alpha is: ", 1/weib_mod$scale)
#Note the scale parameter should be intepreted as 1/(rweibull shape)
#The rweibull function is: https://en.wikipedia.org/wiki/Weibull_distribution
#There, we have shape parameter k, which corresponds to our alpha.
#In the example above, the scale is 1.64, so k = 1/1.64
#Thus, k < 1, i.e. decreasing hazard over time (which makes sense)
# now we add some more regressors
weib_mod_2 <- survreg(Surv(splength, status) ~ gender+ public, data=df, dist="weibull")
print(1/weib_mod_2$scale) #no change
summary(weib_mod_2) #negative effect of being male on sickness, public insignificant
weib_mod_3 <- survreg(Surv(splength, status) ~ gender + urban + public,
data=df, dist="weibull")
print(1/weib_mod_3$scale) #again no difference
summary(weib_mod_3) #effect of male still the same, urban/public is insignificant
weib_mod_4 <- survreg(Surv(splength, status) ~ gender + urban + public + catholic +
protest + special, data=df, dist="weibull")
summary(weib_mod_4) #only male and special schools differ a lot in duration
print(1/weib_mod_3$scale) #still no change
# now lets plot the models separately by gender and see if the parameters are
# significantly different using the log-rank test
weib_mod_gen_1 <- survreg(Surv(splength, status) ~  urban + public + catholic +
protest + special, data=df[df[["gender"]]=="male",], dist="weibull")
summary(weib_mod_gen_1)
weib_mod_gen_1$scale
weib_mod_gen_2 <- survreg(Surv(splength, status) ~ urban + public + catholic +
protest + special, data=df[df[["gender"]]== "female",], dist="weibull")
summary(weib_mod_gen_2)
#weib_diff <- survdiff(Surv(splength, status) ~ gender, data=df, dist="weibull")
#I am not really sure how to test difference between the two,
#but I think just looking at the different scale values already gives us
#a reason why to prefer the gender-separate models
#We also look at 'special', because these seem to have different sickness behaviour as well
#
weib_mod_special_1 <- survreg(Surv(splength, status) ~  urban+gender, data=df[df[["special"]]=="1",], dist="weibull")
summary(weib_mod_special_1)
print(1/weib_mod_special_1$scale)
weib_mod_special_2 <- survreg(Surv(splength, status) ~  urban + gender, data=df[df[["special"]]=="0",], dist="weibull")
print(1/weib_mod_special_2$scale)
#indeed, their shape parameters are far apart
# Code for question 3
# Question 3 - piecewise constant ----------------------------------------------
#We follow: https://data.princeton.edu/wws509/r/recidivism
# Fitting piece-wise linear models
#PWC <- function(breaks, df, controls = NULL){
#splitted <- survSplit(Surv(splength, status) ~ ., data = df,
#cut = breaks, episode = "interval", start = "start")
#splitted <- mutate(splitted, exposure = splength - start,
#interval = factor(interval,
#labels = paste("(", c(1,breaks), ",",
#c(breaks,998), "]", sep=""))) %>%
#rename(events = status)
#if(is.null(controls)){
#PWC <- glm(events ~ interval + offset(log(exposure)), data=splitted, family=poisson)
#} else {
#dep_vars <- paste(c("interval", "offset(log(exposure))", controls), collapse="+")
#formula_string <- as.formula(paste("events", dep_vars, sep="~"))
#PWC <- glm(formula_string, data=splitted, family=poisson)
#}
#maxDuration <- max(df[["splength"]])
#b = coef(PWC)[1: length(unique(splitted$interval))]
#h = exp( b)
#lengths <- c(breaks, maxDuration) - c(1, breaks)
#H = cumsum( lengths*h)
#S = exp(-H)
#output <- list(S, PWC)
#return(output)
#}
#plotPWC <- function(breaks, S, maxDuration, subset){
#xgrid <- c(1,breaks, maxDuration) #returns the cut points
#ygrid <- S #returns the corresponding hazards
#graph <- data.frame(x = seq(1,  maxDuration, by = 0.1))
#for(i in 1:length(graph$x)){
#for(j in 2:length(xgrid)){
#if(graph$x[i] < xgrid[j] & graph$x[i] >= xgrid[j-1]){
#graph$y[i] <- ygrid[[j-1]]
#}
#}
#}
#p <- ggplot(graph[ graph$x < subset, ]) +
#geom_line(aes(x = x, y = y)) +
#xlab("duration") +
#ylab("survival")
#return(p)
#}
#breaks1 <- c(2, 5, 10, 20, 50, 500)
#breaks2 <- c(2, 3,4,5, 8, 16, 32, 64, 128, 256, 512)
#controls <- c("gender", "urban", 'public',"catholic", "protest", "special")
#S1 <- PWC(breaks1, df)[[1]]
#plotPWC(breaks1, S1, maxDuration, subset = 50)
#S2 <- PWC(breaks2, df)[[1]]
#plotPWC(breaks2, S2, maxDuration, subset = 50)
#S3 <- PWC(breaks2, df, controls)[[2]]
#summary(S3) #once again, only Males and Special have significant signs, but now reversed!
# breaks3 <- c(0, 1,2,3,4,5, 10, 15, 20, 50, 100, 500, 998)
# S4 <- pchreg(Surv(splength, status) ~ gender + urban + public + catholic +
#          protest + special, breaks3, data = df)
# S4$beta
# Previous code - piecewise constant ----------------------------------------------
piece_mod_1 <- pchreg(Surv(splength, rcensor) ~ 1, breaks = 4, data=df)
summary(piece_mod_1)
piece_mod_2 <- pchreg(Surv(splength, rcensor) ~ 1, breaks = 20, data=df)
summary(piece_mod_2)
piece_mod_gen_1 <- pchreg(Surv(splength, rcensor) ~ gender, breaks = 20, data=df[df[["gender"]]==1,])
# Setup code - loading the data or packages
rm(list=ls())
variables <- c("gender", "marstat", "contract", "schooltype", "urban", "birthyr", "classize")
# packages ---------------------------------------------------------------------
library(stargazer)
library(survival)
library(haven)
#library(xtable)
#library(RCT)
library(foreign)
library(tidyverse)
library(ggfortify)
library(survminer)
library(splines2)
library(casebase)
# for plodurationing the smooth hazard rate
library(visreg)
# for computing piecewise constant hazard rates
library(pch)
#for unobserved heterogeneity models
library(parfm)
# for computing the cumulative hazards
library(zoo)
# reading in data --------------------------------------------------------------
dat <- "C:\\Users\\kiera\\OneDrive\\Documents\\Tinbergen MPhil\\Applied Microeconometrics\\Applied_micro"
setwd(dat)
df <- as.data.frame(read_dta("FlowSpells.dta"), stringsAsFactors=F)
df <- df[df[["sptype"]] == 2, ] #only focus on the sickness periods
# Hazard and survival functions ------------------------------------------------
hazard <- function(df, duration, censored, t){
numerator <- sum((1-df[[censored]])*ifelse(t < df[[duration]] & df[[duration]] <= t + 1, 1, 0))
denominator <- ifelse(any(df[[duration]] > t), sum(df[[duration]] > t), 1) #can't divide by zero
return(numerator/denominator)
}
survival <- function(survivalTable, t){
return(ifelse(t > 1, prod(1 - survivalTable$hazard[1:(t-1)]), 1))
}
# Question 1 -------------------------------------------------------------------
# functions ---------------------------------------------------------------------
# functions that compute the hazard and survival functions
hazard <- function(df, duration, censored, t){
# computing the numerator in a vectorised way
numerator <- sum((1-df[[censored]])*ifelse(t < df[[duration]] & df[[duration]] <= t + 1, 1, 0))
denominator <- ifelse(any(df[[duration]] > t), sum(df[[duration]] > t), 1)
return(numerator/denominator)
}
survival <- function(survivalTable, t){
# this uses that the survival function is the complement of the cumulative
# hazard
return(ifelse(t > 1, prod(1 - survivalTable$hazard[1:(t-1)]), 1))
}
# function to compute the survival table from the dataset
survTableGen <- function(df, col = FALSE, grp = FALSE){
duration <- 1:max(df$splength)
survivalTable <- data.frame(duration)
if(col == FALSE){ #if no column specified, select all data
survivalTable$hazard <- sapply(duration, hazard, df = df,
duration = "splength", censored = "rcensor")
survivalTable$survival <- sapply(duration, survival, survivalTable = survivalTable)
} else { #if column specified, only select data from the selected group
subset <- df[df[[col]]==grp,]
survivalTable$hazard <- sapply(duration, hazard, df = subset,
duration = "splength", censored = "rcensor")
survivalTable$survival <- sapply(duration, survival, survivalTable = survivalTable)
survivalTable[[col]] <- grp
}
return(survivalTable)
}
# function that computes hazard rates and survival function by groups defined
# by continuous variables
sur_haz_grp <- function(col, df){
groups <- unique(df[[col]])
out_dfs <- lapply(groups, survTableGen, col = col, df = df)
out_df_merged <- bind_rows(out_dfs)
hazard <- ggplot(out_df_merged, aes(x =duration, y=hazard, colour = as.factor(get(col)))) +
geom_line() +
theme(legend.title=element_blank())
survival <- ggplot(out_df_merged, aes(x =duration, y=survival, colour = as.factor(get(col)))) +
geom_line() +
theme(legend.title=element_blank())
# finally testing for difference between the survival rates
return(list(hazard, survival))
}
# running code ----------------------------------------------------------------
survivalTable <- survTableGen(df)
ggplot(data = survivalTable, mapping = aes(x =duration, y= hazard)) +
geom_line() + geom_smooth()
ggplot(data = survivalTable, mapping = aes(x =duration, y=survival)) +
geom_line()
# now lets do plots for the first two weeks, and the first year
ggplot(data = survivalTable[survivalTable[["duration"]]<=14,], mapping = aes(x =duration, y=hazard)) +
geom_line() + geom_smooth()
ggplot(data = survivalTable[survivalTable[["duration"]]<=14,], mapping = aes(x =duration, y=survival)) +
geom_line()
ggplot(data = survivalTable[survivalTable[["duration"]]<=365,], mapping = aes(x =duration, y=hazard)) +
geom_line() + geom_smooth()
ggplot(data = survivalTable[survivalTable[["duration"]]<=365,], mapping = aes(x =duration, y=survival)) +
geom_line()
# now lets partition by different columns and compute the hazard rates
# partitions: type of school, gender, urbanisation
df[["school_type"]] <- NA
df[["school_type"]][df[["public"]]==1] <- "Public"
df[["school_type"]][df[["catholic"]]==1] <- "Catholic"
df[["school_type"]][df[["protest"]]==1] <- "Protestant"
df[["school_type"]][df[["special"]]==1] <- "Special"
df[["area_type"]] <- NA
df[["area_type"]][df[["urban"]]==1] <- "Urban_1"
df[["area_type"]][df[["urban"]]==2] <- "Urban_2"
df[["area_type"]][df[["urban"]]==3] <- "Urban_3"
df[["area_type"]][df[["urban"]]==4] <- "Urban_4"
df[["area_type"]][df[["urban"]]==5] <- "Urban_5"
df[["gender"]] <- ifelse(df[["gender"]] == 1, "male", "female")
#choose to cut off at 300 because from there the pictures are uninformative
sur_haz_grp("gender", df[df[["splength"]] < 300, ])
sur_haz_grp("area_type", df[df[["splength"]] < 300, ])
sur_haz_grp("school_type", df[!is.na(df[["school_type"]]) & df[["splength"]] < 300, ])
# finally lets test for a difference between the survival rates
survdiff(Surv(splength, rcensor) ~ gender, data=df)
survdiff(Surv(splength, rcensor) ~ urban, data=df)
survdiff(Surv(splength, rcensor) ~ school_type, data=df)
# Question 2 - parametric survival models --------------------------------------
# simple weibull and exponential models
# In order to use the package, we have to make a status indicator
df$status <- ifelse(df$rcensor == 1, 0, 1)
exp_mod <- survreg(Surv(splength, status) ~ public, data=df, dist="exponential")
summary(exp_mod)
weib_mod <- survreg(Surv(splength, status) ~ public, data=df, dist="weibull")
summary(weib_mod)
paste("The scale parameter alpha is: ", 1/weib_mod$scale)
#Note the scale parameter should be intepreted as 1/(rweibull shape)
#The rweibull function is: https://en.wikipedia.org/wiki/Weibull_distribution
#There, we have shape parameter k, which corresponds to our alpha.
#In the example above, the scale is 1.64, so k = 1/1.64
#Thus, k < 1, i.e. decreasing hazard over time (which makes sense)
# now we add some more regressors
weib_mod_2 <- survreg(Surv(splength, status) ~ gender+ public, data=df, dist="weibull")
print(1/weib_mod_2$scale) #no change
summary(weib_mod_2) #negative effect of being male on sickness, public insignificant
weib_mod_3 <- survreg(Surv(splength, status) ~ gender + urban + public,
data=df, dist="weibull")
print(1/weib_mod_3$scale) #again no difference
summary(weib_mod_3) #effect of male still the same, urban/public is insignificant
weib_mod_4 <- survreg(Surv(splength, status) ~ gender + urban + public + catholic +
protest + special, data=df, dist="weibull")
summary(weib_mod_4) #only male and special schools differ a lot in duration
print(1/weib_mod_3$scale) #still no change
# now lets plot the models separately by gender and see if the parameters are
# significantly different using the log-rank test
weib_mod_gen_1 <- survreg(Surv(splength, status) ~  urban + public + catholic +
protest + special, data=df[df[["gender"]]=="male",], dist="weibull")
summary(weib_mod_gen_1)
weib_mod_gen_1$scale
weib_mod_gen_2 <- survreg(Surv(splength, status) ~ urban + public + catholic +
protest + special, data=df[df[["gender"]]== "female",], dist="weibull")
summary(weib_mod_gen_2)
#weib_diff <- survdiff(Surv(splength, status) ~ gender, data=df, dist="weibull")
#I am not really sure how to test difference between the two,
#but I think just looking at the different scale values already gives us
#a reason why to prefer the gender-separate models
#We also look at 'special', because these seem to have different sickness behaviour as well
#
weib_mod_special_1 <- survreg(Surv(splength, status) ~  urban+gender, data=df[df[["special"]]=="1",], dist="weibull")
summary(weib_mod_special_1)
print(1/weib_mod_special_1$scale)
weib_mod_special_2 <- survreg(Surv(splength, status) ~  urban + gender, data=df[df[["special"]]=="0",], dist="weibull")
print(1/weib_mod_special_2$scale)
#indeed, their shape parameters are far apart
# Code for question 3
# Question 3 - piecewise constant ----------------------------------------------
#We follow: https://data.princeton.edu/wws509/r/recidivism
# Fitting piece-wise linear models
#PWC <- function(breaks, df, controls = NULL){
#splitted <- survSplit(Surv(splength, status) ~ ., data = df,
#cut = breaks, episode = "interval", start = "start")
#splitted <- mutate(splitted, exposure = splength - start,
#interval = factor(interval,
#labels = paste("(", c(1,breaks), ",",
#c(breaks,998), "]", sep=""))) %>%
#rename(events = status)
#if(is.null(controls)){
#PWC <- glm(events ~ interval + offset(log(exposure)), data=splitted, family=poisson)
#} else {
#dep_vars <- paste(c("interval", "offset(log(exposure))", controls), collapse="+")
#formula_string <- as.formula(paste("events", dep_vars, sep="~"))
#PWC <- glm(formula_string, data=splitted, family=poisson)
#}
#maxDuration <- max(df[["splength"]])
#b = coef(PWC)[1: length(unique(splitted$interval))]
#h = exp( b)
#lengths <- c(breaks, maxDuration) - c(1, breaks)
#H = cumsum( lengths*h)
#S = exp(-H)
#output <- list(S, PWC)
#return(output)
#}
#plotPWC <- function(breaks, S, maxDuration, subset){
#xgrid <- c(1,breaks, maxDuration) #returns the cut points
#ygrid <- S #returns the corresponding hazards
#graph <- data.frame(x = seq(1,  maxDuration, by = 0.1))
#for(i in 1:length(graph$x)){
#for(j in 2:length(xgrid)){
#if(graph$x[i] < xgrid[j] & graph$x[i] >= xgrid[j-1]){
#graph$y[i] <- ygrid[[j-1]]
#}
#}
#}
#p <- ggplot(graph[ graph$x < subset, ]) +
#geom_line(aes(x = x, y = y)) +
#xlab("duration") +
#ylab("survival")
#return(p)
#}
#breaks1 <- c(2, 5, 10, 20, 50, 500)
#breaks2 <- c(2, 3,4,5, 8, 16, 32, 64, 128, 256, 512)
#controls <- c("gender", "urban", 'public',"catholic", "protest", "special")
#S1 <- PWC(breaks1, df)[[1]]
#plotPWC(breaks1, S1, maxDuration, subset = 50)
#S2 <- PWC(breaks2, df)[[1]]
#plotPWC(breaks2, S2, maxDuration, subset = 50)
#S3 <- PWC(breaks2, df, controls)[[2]]
#summary(S3) #once again, only Males and Special have significant signs, but now reversed!
# breaks3 <- c(0, 1,2,3,4,5, 10, 15, 20, 50, 100, 500, 998)
# S4 <- pchreg(Surv(splength, status) ~ gender + urban + public + catholic +
#          protest + special, breaks3, data = df)
# S4$beta
# Previous code - piecewise constant ----------------------------------------------
piece_mod_1 <- pchreg(Surv(splength, rcensor) ~ 1, breaks = 4, data=df)
summary(piece_mod_1)
piece_mod_2 <- pchreg(Surv(splength, rcensor) ~ 1, breaks = 20, data=df)
summary(piece_mod_2)
piece_mod_gen_1 <- pchreg(Surv(splength, rcensor) ~ 1, breaks = 20, data=df[df[["gender"]]==1,])
piece_mod_1 <- pchreg(Surv(splength, rcensor) ~ 1, breaks = 4, data=df)
summary(piece_mod_1)
piece_mod_2 <- pchreg(Surv(splength, rcensor) ~ 1, breaks = 20, data=df)
summary(piece_mod_2)
piece_mod_gen_1 <- pchreg(Surv(splength, rcensor) ~ 1, breaks = 20, data=df[df[["gender"]]==1,])
View(df)
piece_mod_gen_1 <- pchreg(Surv(splength, rcensor) ~ 1, breaks = 20, data=df[df[["gender"]]=="male",])
summary(piece_mod_gen_1)
piece_mod_gen_2 <- pchreg(Surv(splength, rcensor) ~ 1, breaks = 20, data=df[df[["gender"]]=="female",])
summary(piece_mod_gen_2)
install.packages("bookdown")
library(bookdown)
bookdown::render_book("index.Rmd")
install.packages("bookdown")
install.packages("bookdown")
library(bookdown)
bookdown::render_book("index.Rmd")
install.packages("bookdown")
install.packages("bookdown")
library(bookdown)
bookdown::render_book("index.Rmd")