Logit and COX PH Regression of ICU Data(Revised with GGPLOT2)
Kyung Yun
1. Introduction of Intensive Care Unit data.
The ICU data set consists of a sample of 200 subjects who were part of a much larger study on survival of patients following admission to an adult intensive care unit (ICU). The major goal of this study was to predict the probability of survival to hospital discharge of these patients and to study the risk factors associated with ICU mortality. (Hosmer Jr and Sturdivant. 2013)
2. Some variables of Intensive Care Unit data.
There are 200 observations of 21 variables.
STA : Vital Status(0=Lived, 1=Died)
AGE : Years
GENDER: 0=Male, 1=Female
RACE: 1=White, 2=Black, 3=Other
TYP : Type of Admission(0=Elective, 1=Emergency)
LOC : Level of Consciousness at ICU Admission (0=No Coma,or Stupor, 1=Deep Stupor, 2=Coma)
3. Descriptive Statistics using ggplot2
library(plyr)
library(ggplot2)
icu <- read.table("icu.txt", header=TRUE)
icu2 <- icu
icu2$GENDER <-factor(icu2$GENDER, labels = c("male", "female"))
icu2$RACE <-factor(icu2$RACE, labels = c("White", "Black", "Others"))
icu2$STA <-factor(icu2$STA, labels = c("Lived", "Died"))
qplot(GENDER, AGE, data = icu2, geom="boxplot")+ geom_point() +
facet_grid(.~RACE, scales="free",space="free") 
ds <- ddply(icu2, .(STA), summarise, mean = mean(AGE), sd = sd(AGE))
ggplot()+
geom_point(data= icu2, aes(x=STA, y=AGE)) + scale_y_continuous() + scale_x_discrete() + xlab("STATUS") +
geom_point(data=ds, aes(x=STA, y = mean, colour = 'red', size = 3)) +
geom_errorbar(data=ds, aes(x = STA, y = mean, ymin = mean - sd, ymax = mean + sd),
colour = 'red', width = 0.4) +
facet_grid(GENDER ~ RACE, margins = TRUE)
4. Estimation of Logistic Regression Models using “Zelig” Package
library(Zelig)
icu.z1 <- zelig(STA ~ AGE+GENDER+RACE+SER+CAN+CRN+INF+CPR+SYS+HRA+PRE+TYP+FRA+PO2+PH+PCO+BIC+CRE+LOC, data=icu, model="logit", x=TRUE)
icu.z2 <- zelig(STA ~ AGE+TYP, data=icu, model="logit",x=TRUE)
icu.z3 <- zelig(STA ~ AGE+TYP + AGE:TYP, data=icu, model="logit",x=TRUE)library(stargazer)
stargazer(icu.z1, icu.z2, icu.z3, type="html")| Dependent variable: | |||
| STA | |||
| (1) | (2) | (3) | |
| AGE | 0.052*** | 0.034*** | 0.074 |
| (0.017) | (0.011) | (0.088) | |
| GENDER | -0.555 | ||
| (0.501) | |||
| RACE | -0.005 | ||
| (0.528) | |||
| SER | -0.545 | ||
| (0.575) | |||
| CAN | 2.758*** | ||
| (0.980) | |||
| CRN | -0.102 | ||
| (0.762) | |||
| INF | -0.056 | ||
| (0.535) | |||
| CPR | 0.978 | ||
| (0.983) | |||
| SYS | -0.012 | ||
| (0.008) | |||
| HRA | -0.004 | ||
| (0.009) | |||
| PRE | 0.929 | ||
| (0.629) | |||
| TYP | 2.744*** | 2.454*** | 5.308 |
| (0.995) | (0.753) | (6.413) | |
| FRA | 1.152 | ||
| (0.999) | |||
| PO2 | 0.387 | ||
| (0.851) | |||
| PH | 2.415** | ||
| (1.231) | |||
| PCO | -3.173** | ||
| (1.386) | |||
| BIC | -0.794 | ||
| (0.916) | |||
| CRE | 0.233 | ||
| (1.075) | |||
| LOC | 2.706*** | ||
| (0.752) | |||
| AGE:TYP | -0.040 | ||
| (0.088) | |||
| Constant | -5.414** | -5.509*** | -8.315 |
| (2.122) | (1.034) | (6.375) | |
| Observations | 200 | 200 | 200 |
| Log Likelihood | -64.800 | -86.538 | -86.416 |
| Akaike Inf. Crit. | 169.600 | 179.076 | 180.832 |
| Note: | p<0.1; p<0.05; p<0.01 | ||
“Zelig” package is used to estimate logistic regresssion models of patient’s Vital Status(0=Lived, 1=Died) at ICU(Intensive Care Unit) as a function of explanatory variables.(Kosuke Imai and Lau 2015)
The first model uses all the variables in the dataset and shows only 6 variables are significant among 19 explanatory variables.
For the second model, only 2 explanatory variables: AGE and TYP(Type of Admission, 0=Elective, 1=Emergency) are included to explain the relationships of Patient’s age and type of admission on the response variable vital status. The coefficients for AGE is estimated as 0.034 which means that the ICU patient’s log odds of dying(Vital Status =1) increases 0.034 times higher for additonal age increment with type of admission to the ICU being equal. Exponentiating the logistic regression model removes the log for the odds and changes it from a model that is additive in the log-odds scale to one that is multiplicative in the odds scale. Therefore, we also can interpret that when type of admission to the ICU is elective(TYP=0) the ICU patient’s odds of dying(Vital Status =1) increases 1.0346(=e^0.034) times higher for additonal age increment. The third model includes AGE, TYP and their interaction term AGE:TYP to measure the interaction effect of AGE and TYP. The result shows that the interaction term of AGE and TYP is not statistically significant. Then 3 logit models are compared using “stargazer” package.(Hlavac 2014)
tempx <- as.data.frame(icu.z2$x)
tempy <- as.data.frame(icu.z2$fitted.values)
df <- cbind(tempx, tempy)
ggplot(df, aes(x=AGE, y=icu.z2$fitted.values)) + ylab("Fitted Values of Logit Regression") +
geom_point() + stat_smooth( aes(y = icu.z2$fitted.values), method="glm", family="binomial", se=T) 
5. Simulating the logistic models using the Zelig package
Next the graph visualizes the relationship of expected value(Vital Status) and age. As age increases the expected probabilty of dying increases.
library(dplyr)
library(Zelig)
library(ggplot2)
library(reshape2)
library(gridExtra)
library(stargazer)
set.seed(11111)
xrange <- min(icu$AGE):max(icu$AGE)
x.out <- setx(icu.z2, AGE = xrange)
s1 <-sim(icu.z2, x = x.out)
s.mod <- (s1$qi)
Ev <- (s.mod$ev)
Ev <- as.matrix(Ev)
E.mod <- as.data.frame(Ev)
E.mod <- melt(E.mod)
S.mod <- select(E.mod, variable, value)
S.mod <- group_by(S.mod, variable)
ggplot(S.mod, aes(variable, value)) +
geom_point(color = "lightblue", alpha = I(0.05)) +
stat_smooth(aes(group=1), method = "lm") +scale_x_discrete(breaks = c("V20", "V40", "v60"), labels =c("20","40","60"))+
xlab("AGE") + ylab("Expected Values: E(Y|X)")
6. Significance test of the interaction effect
The third logistic model icu.z3 includes the interaction term of AGE and STA. To test the significance of the interaction effect, the following “difference in differences” method is conducted. The Zelig packages versions 4 and 5 had the systematic error with the simulated expected joint probabilities or expected values. Therefore, Zelig 3.5-2 package is used instead.
library(mvtnorm)
set.seed(12345)
xh1 <- setx(icu.z3, AGE = mean(icu$AGE)+sd(icu$AGE), TYP=1)
xl1 <- setx(icu.z3, AGE = mean(icu$AGE), TYP=1)
xh0 <- setx(icu.z3, AGE = mean(icu$AGE)+sd(icu$AGE), TYP =0)
xl0 <- setx(icu.z3, AGE = mean(icu$AGE), TYP=0)
zh1 <- sim(icu.z3, x=xh1)
zl1 <- sim(icu.z3, x=xl1)
zh0 <- sim(icu.z3, x=xh0)
zl0 <- sim(icu.z3, x=xl0)
eff <- (zh1$qi$ev - zl1$qi$ev) -(zh0$qi$ev - zl0$qi$ev)
quantile(eff, c(.025, .975))## 2.5% 97.5%
## -0.1395350 0.3470681As a result, 95% Confidence interval(-0.1395, 0.3471) of the distribution of the random values for the interaction effect includes 0. Which implies the interaction term AGE:TYP has no significant effect on the vital status of patient. This result is the same as the logistic regression using “Zelig” package in (3) above.
eff<- data.frame(eff)
ggplot(eff, aes(x=eff)) +
geom_histogram(colour="black", fill="yellow") +
scale_x_continuous(limits=c(-0.6, 0.8), name="Simulation of Interaction effect with 95% C.I.(with red lines)")+
geom_vline(data=eff, aes(xintercept=quantile(eff,c(.025, .975))),
linetype="dashed", size=1, colour="red")
Histogram of the distribution of the interaction effect measures using “difference in difference methods” gives a better description of no interaction effect. 95% C.I. including 0 implies no significant interaction effect.
7. Cox Proportional Hazard model
Next, the cox proprotional hazards regression model is used to analyze the relationship between the survival time of patients(AGE) and risk factors assessing their conditons at the ICU. Cox PH model is the most common statistical methods for censored data and it does not require any particular probability distribution of survival times. This fact relieves a concern about the unknown distribution of survival times for this ICU dataset. In this Cox Proportional Hazard model, the variable STA is used for censored observation. That is, if STA=1, then the patient is dead at the ICU and the survival duration of the patient is not censored and if STA=0, then the patient is alive and the survival duration of the patient is censored.
z.out0 <- zelig(Surv(AGE, STA) ~ GENDER+RACE+SER+CAN+CRN+INF+CPR+SYS+HRA+PRE+TYP+FRA+PO2+PH+PCO+BIC+CRE+LOC, data = icu, model = "coxph")stargazer(z.out0, type = "html")| Dependent variable: | |
| AGE | |
| GENDER | -0.942** |
| (0.466) | |
| RACE | 0.609 |
| (0.405) | |
| SER | -0.289 |
| (0.497) | |
| CAN | 2.680*** |
| (0.735) | |
| CRN | -0.372 |
| (0.513) | |
| INF | -0.083 |
| (0.414) | |
| CPR | 0.428 |
| (0.620) | |
| SYS | -0.003 |
| (0.006) | |
| HRA | -0.004 |
| (0.008) | |
| PRE | 1.288** |
| (0.570) | |
| TYP | 2.238** |
| (0.871) | |
| FRA | 1.180 |
| (0.897) | |
| PO2 | 0.069 |
| (0.651) | |
| PH | 1.909** |
| (0.845) | |
| PCO | -1.621* |
| (0.863) | |
| BIC | -0.540 |
| (0.686) | |
| CRE | -0.318 |
| (0.637) | |
| LOC | 1.255*** |
| (0.331) | |
| Observations | 200 |
| R2 | 0.251 |
| Max. Possible R2 | 0.798 |
| Log Likelihood | -131.029 |
| Wald Test | 49.750*** (df = 18) |
| LR Test | 57.815*** (df = 18) |
| Score (Logrank) Test | 72.416*** (df = 18) |
| Note: | p<0.1; p<0.05; p<0.01 |
The result of Cox PH regression analysis shows that explanatory variables, GENDER,CAN,PRE,TYP,PH,PCO,LOC are statistically significant for the explanation of the response variable AGE. Next, only set of significant explanatory variables are used to set up a new Cox PH model. Then each explanatory variable is used to explain the relationship between the response variable AGE. From the comparison table, only TYP and LOC become significant.
z.out1 <- zelig(Surv(AGE, STA) ~ GENDER+CAN+PRE+TYP+PH+PCO+LOC, data = icu, model = "coxph")
z.out2 <- zelig(Surv(AGE, STA) ~ GENDER, data = icu, model = "coxph")
z.out3 <- zelig(Surv(AGE, STA) ~ CAN, data = icu, model = "coxph")
z.out4 <- zelig(Surv(AGE, STA) ~ PRE, data = icu, model = "coxph")
z.out5 <- zelig(Surv(AGE, STA) ~ TYP, data = icu, model = "coxph")
z.out6 <- zelig(Surv(AGE, STA) ~ PH, data = icu, model = "coxph")
z.out7 <- zelig(Surv(AGE, STA) ~ LOC, data = icu, model = "coxph")stargazer(z.out1, z.out2, z.out3, z.out4, z.out5, z.out6, z.out7, type = "html")| Dependent variable: | |||||||
| AGE | |||||||
| (1) | (2) | (3) | (4) | (5) | (6) | (7) | |
| GENDER | -0.781** | -0.162 | |||||
| (0.380) | (0.331) | ||||||
| CAN | 2.491*** | 0.439 | |||||
| (0.635) | (0.536) | ||||||
| PRE | 1.183** | 0.218 | |||||
| (0.498) | (0.422) | ||||||
| TYP | 2.266*** | 2.056*** | |||||
| (0.773) | (0.728) | ||||||
| PH | 1.624*** | 0.496 | |||||
| (0.605) | (0.530) | ||||||
| PCO | -1.292** | ||||||
| (0.601) | |||||||
| LOC | 1.308*** | 0.885*** | |||||
| (0.238) | (0.187) | ||||||
| Observations | 200 | 200 | 200 | 200 | 200 | 200 | 200 |
| R2 | 0.222 | 0.001 | 0.003 | 0.001 | 0.071 | 0.004 | 0.080 |
| Max. Possible R2 | 0.798 | 0.798 | 0.798 | 0.798 | 0.798 | 0.798 | 0.798 |
| Log Likelihood | -134.803 | -159.816 | -159.636 | -159.809 | -152.541 | -159.552 | -151.624 |
| Wald Test | 44.700*** (df = 7) | 0.240 (df = 1) | 0.670 (df = 1) | 0.270 (df = 1) | 7.960*** (df = 1) | 0.880 (df = 1) | 22.360*** (df = 1) |
| LR Test | 50.269*** (df = 7) | 0.242 (df = 1) | 0.602 (df = 1) | 0.256 (df = 1) | 14.793*** (df = 1) | 0.770 (df = 1) | 16.627*** (df = 1) |
| Score (Logrank) Test | 60.157*** (df = 7) | 0.241 (df = 1) | 0.682 (df = 1) | 0.269 (df = 1) | 11.149*** (df = 1) | 0.894 (df = 1) | 28.120*** (df = 1) |
| Note: | p<0.1; p<0.05; p<0.01 | ||||||
z.out11 <- zelig(Surv(AGE, STA) ~ GENDER+CAN+PRE+TYP+PH+PCO+LOC, data = icu, model = "coxph")
z.out21 <- zelig(Surv(AGE, STA) ~ TYP+LOC, data = icu, model = "coxph")
z.out31 <- zelig(Surv(AGE, STA) ~ TYP+LOC+TYP:LOC, data = icu, model = "coxph")
z.out41 <- zelig(Surv(AGE, STA) ~ TYP:LOC, data = icu, model="coxph")stargazer(z.out11, z.out21, z.out31, z.out41, type = "html")| Dependent variable: | ||||
| AGE | ||||
| (1) | (2) | (3) | (4) | |
| GENDER | -0.781** | |||
| (0.380) | ||||
| CAN | 2.491*** | |||
| (0.635) | ||||
| PRE | 1.183** | |||
| (0.498) | ||||
| TYP | 2.266*** | 1.814** | 2.485** | |
| (0.773) | (0.736) | (1.021) | ||
| PH | 1.624*** | |||
| (0.605) | ||||
| PCO | -1.292** | |||
| (0.601) | ||||
| LOC | 1.308*** | 0.732*** | 3.975*** | |
| (0.238) | (0.188) | (1.421) | ||
| TYP:LOC | -3.275** | 0.859*** | ||
| (1.432) | (0.189) | |||
| Observations | 200 | 200 | 200 | 200 |
| R2 | 0.222 | 0.125 | 0.142 | 0.074 |
| Max. Possible R2 | 0.798 | 0.798 | 0.798 | 0.798 |
| Log Likelihood | -134.803 | -146.564 | -144.623 | -152.273 |
| Wald Test | 44.700*** (df = 7) | 24.320*** (df = 2) | 22.690*** (df = 3) | 20.690*** (df = 1) |
| LR Test | 50.269*** (df = 7) | 26.746*** (df = 2) | 30.629*** (df = 3) | 15.328*** (df = 1) |
| Score (Logrank) Test | 60.157*** (df = 7) | 34.533*** (df = 2) | 36.476*** (df = 3) | 25.817*** (df = 1) |
| Note: | p<0.1; p<0.05; p<0.01 | |||
8. Simulation of Cox PH model
The depedent variable of Cox PH model is the set of Surv(AGE, STA) where AGE is the duration of ICU patients andSTA is a binary status such that STA =0 if the patient is alive(censored), and STA=1 if the patient is dead(not censored). 2 explatory variables TYP(Type of admission) and LOC(Level of Consciousness at ICU Admission, 0=No Coma,or Stupor, 1=Deep Stupor, 2=Coma) and their interaction effect term are used for the survival function.
The PH regression result also shows that each coefficient of the predictors is significant at 0.01 levels. The coefficient estimate of TYP:LOC has negative value, -3.27508. This means that increase of TYP:LOC gives negative effect for death of a patient unlike the other predictors. As a result, when one more unit of TYP:LOC is increased, a patient is 26.45(=1/0.03781) times more likely to survive.
However, interpretation of the interaction effect has problems because of the collinearity and characteristics of 2 categorical variables. In addition, the interaction term alone has postive coefficient value for the death of a patient in the model 4 in the above table although it has negative coefficient value in the model 3 in the above table. However, for the other variables the increased amount implies to have more risk of death. So the TYP=1(Emergency admission) increases the risk of death 12.00094 times higher than TYP=0(elctive admission), and LOC has much higher hazard rates than TYP. The increase of one level of LOC increases the risk of death 53.27601 higher(Patients with No Coma VS Patients with Stupor or Patients with Stupor VS Patients with Coma).
z.surv <- zelig(Surv(AGE, STA) ~ TYP+LOC+TYP:LOC, data = icu, model = "coxph")summary(z.surv)## Call:
## zelig(formula = Surv(AGE, STA) ~ TYP + LOC + TYP:LOC, model = "coxph",
## data = icu)
##
## n= 200, number of events= 40
##
## coef exp(coef) se(coef) z Pr(>|z|)
## TYP 2.48499 12.00094 1.02061 2.435 0.01490 *
## LOC 3.97549 53.27601 1.42054 2.799 0.00513 **
## TYP:LOC -3.27508 0.03781 1.43220 -2.287 0.02221 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## TYP 12.00094 0.08333 1.623582 88.7067
## LOC 53.27601 0.01877 3.291290 862.3773
## TYP:LOC 0.03781 26.44531 0.002283 0.6262
##
## Concordance= 0.737 (se = 0.048 )
## Rsquare= 0.142 (max possible= 0.798 )
## Likelihood ratio test= 30.63 on 3 df, p=1.018e-06
## Wald test = 22.69 on 3 df, p=4.683e-05
## Score (logrank) test = 36.48 on 3 df, p=5.939e-089. Porportional Hazard Assumption Check
To evaluate the Assumption of Cox Proportional Hazards model, the plot of the log-Negative Survival functin(log of cumulative hazard) vs. log-time is checked. As a result, 2 curves for the different groups for TYP look parallel which implies that hazards are proportional and the assumption is satisfied.
library(survival)
library(plyr)
coxph.mod2 <- coxph(Surv(AGE, STA == 1) ~ TYP, data = icu,
ties = 'breslow')
icu2 <- with(icu, data.frame(TYP = c(0, 1)))
tfit.add <- survfit(coxph.mod2, newdata = icu2)
df1 <- data.frame(
time = tfit.add[1, ]$time,
n.risk = tfit.add[1, ]$n.risk,
n.event = tfit.add[1, ]$n.event,
surv = tfit.add[1, ]$surv,
strata = "0",
upper = tfit.add[1, ]$upper,
lower = tfit.add[1, ]$lower,
log.surv = log(-log(tfit.add[1, ]$surv))
)
df2 <- data.frame(
time = tfit.add[2, ]$time,
n.risk = tfit.add[2, ]$n.risk,
n.event = tfit.add[2, ]$n.event,
surv = tfit.add[2, ]$surv,
strata = "1",
upper = tfit.add[2, ]$upper,
lower = tfit.add[2, ]$lower,
log.surv = log(-log(tfit.add[2, ]$surv))
)
dfpar.add <- rbind(df1, df2)
zeros <- data.frame(time = 0, surv = 1, strata = c(1, 2),
upper = 1, lower = 1)
dfpar.add <- rbind.fill(zeros, dfpar.add)
dfpar.add$strata <- factor(dfpar.add$strata, labels = c("Elective", "Emergency"))
ggplot(dfpar.add, aes(log(time), log.surv, colour = strata)) +
geom_step(size = 0.6) +
scale_color_manual("Type of Admission", values = c('blue4', 'darkorange')) +
xlab("ln(time)") + ylab("Log-log survival")
10. Predicted Survival plots
Next graph shows the comparison of the plots of predcited survival times from a Cox PH model to Kaplan-Meier Survivor functin which do not assume Proportional Hazards.
tfit.km <- survfit(Surv(AGE, STA == 1) ~ TYP, data = icu)
df3.km <- data.frame(
time = tfit.km$time,
n.risk = tfit.km$n.risk,
n.event = tfit.km$n.event,
surv = tfit.km$surv,
strata = gsub("TYP=", "", summary(tfit.km, censored = T)$strata),
upper = tfit.km$upper,
lower = tfit.km$lower
)
zeros <- data.frame(time = 0, surv = 1, strata = gsub("TYP=", "",
levels(summary(tfit.km)$strata)),
upper = 1, lower = 1)
df3.km <- rbind.fill(df3.km, zeros)
df3.km$cat <- with(df3.km, ifelse(strata == "0", "Elective, Observed",
"Emergency, Observed"))
dfpar.add$cat <- with(dfpar.add, ifelse(strata == "1", "Elective, Expected",
"Emergency, Expected"))
dfpar.obs <- rbind.fill(dfpar.add, df3.km)
ggplot(dfpar.obs, aes(time, surv, colour = cat)) +
geom_step(size = 0.6) +
scale_color_manual("", values = c('blue1', 'green4', 'red1',
'blue4')) +
xlab("time") + ylab("survival probability")
11. Evaluating the Overall Fit of the Model
The Cox and Snell residuals can be used to assess the fit of a model based on the Cox proportional hazards model. The plot shows that the estimated cumulative hazard of the residuals doesn’t fit to the modified Cox-Snell residuals. The possible reason is that the data has no variable for the survival time(obsered time). The variable AGE is just used for the survival time.
coxph.mod <- coxph(Surv(AGE, STA ==1) ~ TYP+LOC+TYP:LOC,
data = icu, ties = 'breslow')
cox.snell <- (as.numeric(icu$STA) - 1) - resid(coxph.mod,
type = "martingale")
coxph.res <- survfit(coxph(Surv(cox.snell, icu$STA == 1) ~ 1,
method = 'breslow'), type = 'aalen')
plot(coxph.res$time, -log(coxph.res$surv), type = 's',
xlab = 'Modified Cox-Snell residuals', ylab = 'Cumulative hazard')
abline(0, 1, col = 'red', lty = 2)
12. The simulated Hazard Ratios with Zelig and GGPLOt2
The next simulation of the survival function starts with defining the range of TYP variable(TYP=0, TYP=1). Then simulation shows that simuated hazard ratio distribution “s.surv\(qi\)hr” returns the estimation of the hazard ratio for the TYP variable. So the result shows that mean hazard rate is 11.5919 which implies that the patients with emergency admission has 11.5919 times higher risk of death than the patiens with elective admission. This results is very close to the hazard rate estimated in the PH model using Zelig.
The plots also visualizes the simulated survival function and estimated cumulated hazard and hazard rate over time. The plots show the typical behavior of survival function. As time goes by, the survival rate falls to 0, the cumulative hazard increases while the hazard rate being constant over time, which is required by the Cox PH model assumption.
x.survlow1 <- setx(z.surv, TYP=0)
x.survhigh1 <- setx(z.surv, TYP=1)
s.surv <- sim(z.surv, x=x.survlow1, x1=x.survhigh1)
summary(s.surv$qi$hr)## 1
## Min. : 0.6441
## 1st Qu.: 4.6084
## Median : 8.0260
## Mean : 11.5919
## 3rd Qu.: 14.0998
## Max. :209.0712plot(s.surv)
Next simulated Hazard Ratio plot was created by the GGPLOT2 and it also shows the same right-skewed distribution with mean 11.59, median=8.03. It is exactly matched with the above Hazard Ratios plot.
# Create data.frame
si.s <- (s.surv$qi)
si.s <- data.frame(si.s$hr)
hr <- melt(si.s)
ggplot(hr, aes(x=hr$value)) + geom_histogram(colour="black", fill="yellow") + xlab("Simulated Hazard Ratio")
13. Interaction Effect of Cox PH model
The interaction term LOC:TYP is checked by the simulation.
library(mvtnorm)
set.seed(123)
xh1 <- setx(z.surv, LOC=1, TYP=1)
xl1 <- setx(z.surv, LOC=0, TYP=1)
xh0 <- setx(z.surv, LOC=1, TYP=0)
xl0 <- setx(z.surv, LOC=0, TYP=0)
zh1 <- sim(z.surv, x=xh1)
zl1 <- sim(z.surv, x=xl1)
zh0 <- sim(z.surv, x=xh0)
zl0 <- sim(z.surv, x=xl0)
eff <- (zh1$qi$survival - zl1$qi$survival) -(zh0$qi$survival - zl0$qi$survival)
quantile(eff, c(.025, .975))## 2.5% 97.5%
## -2.3851564 0.8367467The above Zelig simulation shows that the 95% Confidence interval of the interaction term effect is (-2.39, 0.84). Which implies that there is no significant interaction effect of TYP:LOC term. Next GGPlot repeats the same result as the 95% confidence interval (-2.39, 0.84) of the distribution of the interaction effect measured using “diff in diff” method. The red dotted line makes the lower and upper limits of 95% CI.
# Create data.frame
Mod.s1 <- data.frame(eff)
Mod.s1 <- melt(Mod.s1)
ggplot(Mod.s1, aes(x=Mod.s1$value)) +
geom_histogram(colour="black", fill="yellow") +
scale_x_continuous(limits=c(-4, 4), name="Simulated effects of Interaction term of Cox PH Regression(95% C.I. with red lines)")+
geom_vline(data=Mod.s1, aes(xintercept=quantile(eff,c(.025, .975))),
linetype="dashed", size=1, colour="red")
References
Hlavac, Marek(2014). 2014. Stargazer:LaTex Code and ASCII Text for Well-Formatted Regression and Summary Staistics Tables. http://CRAN.R-project.org/package=stargazer.
Hosmer Jr, Stanley Lemeshow, David W., and Rodney X. Sturdivant. 2013. Applied Logistic Regression. John Wiley & Sons.
Kosuke Imai, Gary King, and Olivia Lau. 2015. Logit: Logistic Regression for Dichotomous Dependent Variables. http://gking.harvard.edu/zelig.