Nothing
options(na.action=na.exclude) # preserve missings
options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
library(survival)
#
# The residual methods treat a sparse frailty as a fixed offset with
# no variance
#
aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
kfit1 <- coxph(Surv(time, status) ~ age + sex +
frailty(id, dist='gauss'), kidney)
tempf <- predict(kfit1, type='terms')[,3]
temp <- kfit1$frail[match(kidney$id, sort(unique(kidney$id)))]
#all.equal(unclass(tempf), unclass(temp))
all.equal(as.vector(tempf), as.vector(temp))
# Now fit a model with explicit offset
kfitx <- coxph(Surv(time, status) ~ age + sex + offset(tempf),kidney,
eps=1e-7)
# These are not always precisely the same, due to different iteration paths
aeq(kfitx$coef, kfit1$coef)
# This will make them identical
kfitx <- coxph(Surv(time, status) ~ age + sex + offset(temp),kidney,
iter=0, init=kfit1$coef)
aeq(resid(kfit1), resid(kfitx))
aeq(resid(kfit1, type='score'), resid(kfitx, type='score'))
aeq(resid(kfit1, type='schoe'), resid(kfitx, type='schoe'))
# These are not the same, due to a different variance matrix
# The frailty model's variance is about 2x the naive "assume an offset" var
# Expect a value of about 0.5
aeq(resid(kfit1, type='dfbeta'), resid(kfitx, type='dfbeta'))
# Force equality
zed <- kfitx
zed$var <- kfit1$var
aeq(resid(kfit1, type='dfbeta'), resid(zed, type='dfbeta'))
# The score residuals are equal, however.
temp1 <- resid(kfit1, type='score')
temp2 <- resid(kfitx, type='score')
aeq(temp1, temp2)
#
# Now for some tests of predicted values
#
aeq(predict(kfit1, type='expected'), predict(kfitx, type='expected'))
aeq(predict(kfit1, type='lp'), predict(kfitx, type='lp'))
temp1 <- predict(kfit1, type='terms', se.fit=T)
temp2 <- predict(kfitx, type='terms', se.fit=T)
aeq(temp1$fit[,1:2], temp2$fit)
# the next is not equal, all.equal returns a character string in that case
is.character(aeq(temp1$se.fit[,1:2], temp2$se.fit))
mean(temp1$se.fit[,1:2]/ temp2$se.fit)
aeq(as.vector(temp1$se.fit[,3])^2,
as.vector(kfit1$fvar[match(kidney$id, sort(unique(kidney$id)))]))
print(temp1)
kfit1
kfitx
rm(temp1, temp2, kfitx, zed, tempf)
#
# The special case of a single sparse frailty
#
kfit1 <- coxph(Surv(time, status) ~ frailty(id, dist='gauss'), kidney)
tempf <- predict(kfit1, type='terms')
temp <- kfit1$frail[match(kidney$id, sort(unique(kidney$id)))]
all.equal(as.vector(tempf), as.vector(temp))
# Now fit a model with explicit offset
kfitx <- coxph(Surv(time, status) ~ offset(tempf),kidney, eps=1e-7)
aeq(resid(kfit1), resid(kfitx))
aeq(resid(kfit1, type='deviance'), resid(kfitx, type='deviance'))
#
# Some tests of predicted values
#
aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
aeq(predict(kfit1, type='expected'), predict(kfitx, type='expected'))
aeq(predict(kfit1, type='lp'), predict(kfitx, type='lp'))
temp1 <- predict(kfit1, type='terms', se.fit=T)
aeq(temp1$fit, kfitx$linear)
aeq(temp1$se.fit^2,
kfit1$fvar[match(kidney$id, sort(unique(kidney$id)))])
temp1
kfit1
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.