# tests/fr_kidney.R In survival: Survival Analysis

```options(na.action=na.exclude) # preserve missings
options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
library(survival)

# From:	McGilchrist and Aisbett, Biometrics 47, 461-66, 1991
# Data on the recurrence times to infection, at the point of insertion of
#  the catheter, for kidney patients using portable dialysis equipment.
#  Catheters may be removed for reasons other than infection, in which case
#  the observation is censored.  Each patient has exactly 2 observations.

# Variables: patient, time, status, age,
#	   sex (1=male, 2=female),
#	   disease type (0=GN, 1=AN, 2=PKD, 3=Other)
#	   author's estimate of the frailty
aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)

# I don't match their answers, and I think that I'm right
kfit <- coxph(Surv(time, status)~ age + sex + disease + frailty(id), kidney)
kfit1<- coxph(Surv(time, status) ~age + sex + disease +
frailty(id, theta=1), kidney, iter=20)
kfit0 <- coxph(Surv(time, status)~ age + sex + disease, kidney)
temp <-  coxph(Surv(time, status) ~age + sex + disease +
frailty(id, theta=1, sparse=F), kidney)

# Check out the EM based score equations
#  temp1 and kfit1 should have essentially the same coefficients
#  temp2 should equal kfit1\$frail
# equality won't be exact because of the different iteration paths
temp1 <- coxph(Surv(time, status) ~ age + sex + disease +
offset(kfit1\$frail[id]), kidney)
rr <- tapply(resid(temp1), kidney\$id, sum)
temp2 <- log(rr/1 +1)
aeq(temp1\$coef, kfit1\$coef, tolerance=.005)
aeq(temp2, kfit1\$frail, tolerance=.005)

kfit
kfit1
kfit0
temp

#
# Now fit the data using REML
#
kfitm1 <- coxph(Surv(time,status) ~ age + sex + disease +
frailty(id, dist='gauss'), kidney)
kfitm2 <- coxph(Surv(time,status) ~ age + sex + disease +
frailty(id, dist='gauss', sparse=F), kidney)
kfitm1
summary(kfitm2)
#
# Fit the kidney data using AIC
#

# gamma, corrected aic
coxph(Surv(time, status) ~ age + sex + frailty(id, method='aic', caic=T),
kidney)

coxph(Surv(time, status) ~ age + sex + frailty(id, dist='t'), kidney)
coxph(Surv(time, status) ~ age + sex + frailty(id, dist='gauss', method='aic',
caic=T), kidney)

# uncorrected aic
coxph(Surv(time, status) ~ age + sex + frailty(id, method='aic', caic=F),
kidney)

coxph(Surv(time, status) ~ age + sex + frailty(id, dist='t', caic=F), kidney)
```

## Try the survival package in your browser

Any scripts or data that you put into this service are public.

survival documentation built on Feb. 16, 2023, 7:34 p.m.