HW5
harrissun
2024-05-07
Consider the data (infection.csv) on the times until staphylococcus infection of burn patients.
library(readr)
dt <- read_csv("D:/infection.csv")## Rows: 154 Columns: 18
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (18): Obs, Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8, Z9, Z10, Z11, T1, D1, T2, D2,...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.head(dt)## # A tibble: 6 × 18
## Obs Z1 Z2 Z3 Z4 Z5 Z6 Z7 Z8 Z9 Z10 Z11 T1
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 0 0 15 0 0 1 1 0 0 2 12
## 2 2 0 0 1 20 0 0 1 0 0 0 4 9
## 3 3 0 0 1 15 0 0 0 1 1 0 2 13
## 4 4 0 0 0 20 1 0 1 0 0 0 2 11
## 5 5 0 0 1 70 1 1 1 1 0 0 2 28
## 6 6 0 0 1 20 1 0 1 0 0 0 4 11
## # ℹ 5 more variables: D1 <dbl>, T2 <dbl>, D2 <dbl>, T3 <dbl>, D3 <dbl>Q1
What is the number of patients at each burn site (Z5-Z10)? Which burn site has higher infection rate?
len <- nrow(dt)
z51 <- dt[dt$Z5 == 1,]
z61 <- dt[dt$Z6 == 1,]
z71 <- dt[dt$Z7 == 1,]
z81 <- dt[dt$Z8 == 1,]
z91 <- dt[dt$Z9 == 1,]
z101 <- dt[dt$Z10 == 1,]
cat(paste(cat('the # of Z5:',sum(dt['Z5']),', Infection rate:',sum(z51['D3'])/sum(dt['Z5']),'\n'),
cat('the # of Z6:',sum(dt['Z6']),', Infection rate:',sum(z61['D3'])/sum(dt['Z6']),'\n'),
cat('the # of Z7:',sum(dt['Z7']),', Infection rate:',sum(z71['D3'])/sum(dt['Z7']),'\n'),
cat('the # of Z8:',sum(dt['Z8']),', Infection rate:',sum(z81['D3'])/sum(dt['Z8']),'\n'),
cat('the # of Z9:',sum(dt['Z9']),', Infection rate:',sum(z91['D3'])/sum(dt['Z9']),'\n'),
cat('the # of Z10:',sum(dt['Z10']),', Infection rate:',sum(z101['D3'])/sum(dt['Z10']),'\n'),
cat('=> Z6 has the highest infection rate.')))## the # of Z5: 70 , Infection rate: 0.3428571
## the # of Z6: 35 , Infection rate: 0.4
## the # of Z7: 130 , Infection rate: 0.3230769
## the # of Z8: 63 , Infection rate: 0.3174603
## the # of Z9: 47 , Infection rate: 0.2765957
## the # of Z10: 45 , Infection rate: 0.3777778
## => Z6 has the highest infection rate.Q2
Please show the KM curves for different type of burn Z11.
library(survival)
library(survminer)## 載入需要的套件：ggplot2## 載入需要的套件：ggpubr##
## 載入套件：'survminer'## 下列物件被遮斷自 'package:survival':
##
## myelomalibrary(ggplot2)
library(ggpubr)
dt$Z11 <- factor(dt$Z11)
ggsurvplot(survfit(Surv(T3,D3) ~ Z11 ,data=dt),xlab='Times',ylab='Survival probability')
Q3
Please test if using different treatments (Z1) results in different survival probabilities. What is the median survival time when using routine bathing?
survdiff(Surv(T3,D3)~Z1, data=dt)## Call:
## survdiff(formula = Surv(T3, D3) ~ Z1, data = dt)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Z1=0 70 28 21.4 2.07 3.79
## Z1=1 84 20 26.6 1.66 3.79
##
## Chisq= 3.8 on 1 degrees of freedom, p= 0.05Under \(\alpha\) = 0.05, it cannot reject H0. We cannot conclude that there is significant survival difference between routine bathing and Body cleansing.
bath <- dt[dt$Z1 == 0,]
cat('The median of survival time for routine bathing :',median(bath$T3))## The median of survival time for routine bathing : 17Q4
Please use the Cox regression model to find the relationship between the time to straphylocous aureaus infection and covariates Z1 and Z11. Note that Z11 (Type of burn) is a categorical variable.
dt$Z1 <- factor(dt$Z1)
dt$Z11 <- factor(dt$Z11)
cox_model <- coxph(Surv(T3, D3) ~ Z1 + Z11, data = dt)
summary(cox_model)## Call:
## coxph(formula = Surv(T3, D3) ~ Z1 + Z11, data = dt)
##
## n= 154, number of events= 48
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Z11 -0.5958 0.5511 0.2968 -2.008 0.0447 *
## Z112 1.1328 3.1044 1.0828 1.046 0.2955
## Z113 2.2660 9.6407 1.0837 2.091 0.0365 *
## Z114 0.9888 2.6879 1.0160 0.973 0.3305
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Z11 0.5511 1.8146 0.3081 0.9859
## Z112 3.1044 0.3221 0.3718 25.9197
## Z113 9.6407 0.1037 1.1526 80.6352
## Z114 2.6879 0.3720 0.3669 19.6898
##
## Concordance= 0.641 (se = 0.042 )
## Likelihood ratio test= 11.54 on 4 df, p=0.02
## Wald test = 12.55 on 4 df, p=0.01
## Score (logrank) test = 13.87 on 4 df, p=0.008Comparing ratio ==> If there is no relationship between two indicators, the ratio will close to 1. Above this table, we can see that the CIs of Z11=2 and Z11=4 include 1 and the others don’t, so we can conclude that there is a significant effect when one of Z1=1 or Z11=3 change. # Q5
Check if the proportion assumption is satisfied.
cox_zph <- cox.zph(cox_model)
cox_zph## chisq df p
## Z1 0.44 1 0.507
## Z11 8.14 3 0.043
## GLOBAL 9.09 4 0.059There is the only assumption for Cox regression is the ratio of hazard rate for random two people is a constant.
plot(cox_zph)
Q6
Based on the above question, what’s your conclusion (hazard ratios and confidence intervals)? Please provide some suggestions to the doctors. Above this table, the P-value of indicator Z1 is upper than 0.05 and the indicator Z11’s one isn’t, so Z1 satisfies the assumption and Z11 doesn’t. But the global model slightly not reject H0, so this model is still valid. Theoretically, this model pass the assumption. However, because of the bad performance of the indicator Z11, the validation will decrease. If we wanna continue this project, above this table, we should filter outliers from the indicator Z11 first.
Q7
Please provide a nomogram for users.
library(rms)## Warning: 套件 'rms' 是用 R 版本 4.2.3 來建造的## 載入需要的套件：Hmisc## Warning: 套件 'Hmisc' 是用 R 版本 4.2.3 來建造的##
## 載入套件：'Hmisc'## 下列物件被遮斷自 'package:base':
##
## format.pval, unitsdd <- datadist(dt)
options(datadist="dd")
surv <- Surv(dt$T3,dt$D3)
fp <- psm(surv ~ Z1 + Z11, data = dt, dist='lognormal')
med <- Quantile(fp)
survival <- Survival(fp)
nomM <- nomogram(fp, fun=function(x) med(lp=x),
funlabel="Median Survival Time")
plot(nomM)