DEM 7223 - Event History Analysis - Parametric Hazard Models
Import the Dataset
Code
tracts.texas.pharmacy <- import("tracts.texas.pharmacy.rds") %>%
mutate(phy_type=ifelse(phy_type %in%c("Community Indep","Community Independent"), "Community Indep",
ifelse(phy_type=="", "Unknown",
ifelse(phy_type=="Comm Multi 5","Community Multi",
ifelse(phy_type%in%c("Hospital (Independent)"),"Indep Hospital",
ifelse(phy_type=="Hospital (Multiple/Chain)", "Multi Hospital",
ifelse(phy_type=="Ambulatory Surgical Center","Amb Surg Center",phy_type))))))) %>%
mutate(days=(as.numeric(status_date - lic_orig_date)), rurality=as.factor(rurality), phy_type=as.factor(phy_type)) Age at Pharmacy Closure in Texas
This analysis uses pharmacies’ age-at-closure as the outcome variable. I calculated age at pharmacy closure as the difference between the Pharmacy date of closure (Status_date) and their first date of operation(lic_orig_date). However, in the case of active pharmacies, I calculated their age using the censoring date of pharmacy closure follow-up. In this case, the latest follow-up date was December 31, 2021
Code
tracts.texas.pharmacy <- tracts.texas.pharmacy %>%
mutate(death_age = days,
d.event = ifelse(status == "Closed", 1, 0)) %>%
filter(complete.cases(.))
head(tracts.texas.pharmacy)Code
fit.haz.km <- kphaz.fit(tracts.texas.pharmacy$death_age, tracts.texas.pharmacy$d.event,
method="product-limit")
#smoothed hazard using a Kernal-density method
fit.haz.sm <- muhaz(tracts.texas.pharmacy$death_age, tracts.texas.pharmacy$d.event)
kphaz.plot(fit.haz.km,
main="Plot of the harzard of a Pharmacy closing")
lines(fit.haz.sm, col=2, lwd=3)
legend("topleft",
legend=c("KM Harzard", "Smoothed Hazard"),
col=c(1,2),
lty=c(1,1))
Since i am mainly interested in knowing if the hazard rate of pharmacy closure will increase or decrease giving the types of pharmacy and the location of the pharmacy(metro or non-metro), i will fit a proportional hazard model. For these models, the covariates are the types of pharmacy and residential characteristics of the pharmacy(metro or non metro)
Fitting the Exponential distribution model
Code
fit.exp <- phreg(Surv(death_age,d.event)~ phy_type+rurality+phy_type*rurality,
data=tracts.texas.pharmacy,
dist="weibull",
shape = 1)
summary(fit.exp)Single term deletions
Model:
Surv(death_age, d.event) ~ phy_type + rurality + phy_type * rurality
Df AIC LRT Pr(>Chi)
<none> 29758
phy_type:rurality 7 29762 18.8 0.0089 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Covariate Mean Coef Rel.Risk S.E. Wald p
phy_type
Amb Surg Center 0.057 0 1 (reference)
Community Indep 0.255 1.727 5.625 0.170 0.0000
Community Multi 0.503 -0.281 0.755 0.179 0.1168
Indep Hospital 0.054 0.417 1.517 0.233 0.0732
Multi Hospital 0.035 -0.604 0.547 0.345 0.0797
Other 0.057 1.316 3.729 0.189 0.0000
Public Health 0.036 0.258 1.294 0.267 0.3339
Unknown 0.003 0.106 1.112 0.726 0.8835
rurality
Metro 0.854 0 1 (reference)
Non-Metro 0.146 0.235 1.265 0.601 0.6956
phy_type:rurality
Community Indep:Non-Me -1.382 0.251 0.618 0.0253
Community Multi:Non-Me -0.443 0.642 0.637 0.4868
Indep Hospital:Non-Met -1.515 0.220 0.766 0.0480
Multi Hospital:Non-Met -29.703 0.000 1749763.643 1.0000
Other:Non-Metro -1.103 0.332 0.715 0.1230
Public Health:Non-Metr -1.250 0.287 0.859 0.1457
Unknown:Non-Metro 0.659 1.932 1.364 0.6292
Events 1370
Total time at risk 39369349
Max. log. likelihood -14863
LR test statistic 1143.08
Degrees of freedom 15
Overall p-value 0Code
plot(fit.exp)
Code
AIC(fit.exp)[1] 29757.58Fitting the Weibull Model
Code
fit.weibull <- phreg(Surv(death_age,d.event)~ phy_type+rurality+phy_type*rurality,
data=tracts.texas.pharmacy,
dist="weibull",
)
summary(fit.weibull)Single term deletions
Model:
Surv(death_age, d.event) ~ phy_type + rurality + phy_type * rurality
Df AIC LRT Pr(>Chi)
<none> 29744
phy_type:rurality 7 29748 17.9 0.012 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Covariate Mean Coef Rel.Risk S.E. Wald p
phy_type
Amb Surg Center 0.057 0 1 (reference)
Community Indep 0.255 1.713 5.548 0.170 0.0000
Community Multi 0.503 -0.261 0.770 0.179 0.1460
Indep Hospital 0.054 0.423 1.526 0.233 0.0690
Multi Hospital 0.035 -0.585 0.557 0.345 0.0898
Other 0.057 1.318 3.736 0.189 0.0000
Public Health 0.036 0.291 1.337 0.267 0.2765
Unknown 0.003 0.129 1.137 0.727 0.8593
rurality
Metro 0.854 0 1 (reference)
Non-Metro 0.146 0.260 1.297 0.601 0.6652
phy_type:rurality
Community Indep:Non-Me -1.370 0.254 0.618 0.0265
Community Multi:Non-Me -0.466 0.628 0.637 0.4645
Indep Hospital:Non-Met -1.506 0.222 0.766 0.0494
Multi Hospital:Non-Met -30.115 0.000 2137719.468 1.0000
Other:Non-Metro -1.103 0.332 0.715 0.1231
Public Health:Non-Metr -1.295 0.274 0.859 0.1316
Unknown:Non-Metro 0.639 1.895 1.364 0.6393
Events 1370
Total time at risk 39369349
Max. log. likelihood -14855
LR test statistic 1082.57
Degrees of freedom 15
Overall p-value 0Code
plot(fit.weibull)
Code
AIC(fit.weibull)[1] 29744.46Graphical checks on the model fit
Code
emp<-coxreg(Surv(death_age, d.event)~ phy_type + rurality + phy_type*rurality,
data=tracts.texas.pharmacy)
check.dist(sp=emp,pp=fit.exp, main = "Empirical vs. Exponential")
Code
check.dist(sp=emp,pp=fit.weibull, main = "Empirical vs. Weibull")
Using the AIC for the model fits and the graphical representations of the fit, i concluded that the Weibull Model(AIC=29744.46) had a best fit for the data. Community independent pharmacies are at an increase risk of closure relative to other types of pharmacies.