R Under development (unstable) (2023-08-12 r84939 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > 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)) [1] TRUE > > # 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) [1] TRUE > > # 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)) [1] TRUE > aeq(resid(kfit1, type='score'), resid(kfitx, type='score')) [1] TRUE > aeq(resid(kfit1, type='schoe'), resid(kfitx, type='schoe')) [1] TRUE > > # 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')) [1] "Mean relative difference: 0.5216263" > > # Force equality > zed <- kfitx > zed$var <- kfit1$var > aeq(resid(kfit1, type='dfbeta'), resid(zed, type='dfbeta')) [1] TRUE > > # The score residuals are equal, however. > > temp1 <- resid(kfit1, type='score') > temp2 <- resid(kfitx, type='score') > aeq(temp1, temp2) [1] TRUE > > # > # Now for some tests of predicted values > # > aeq(predict(kfit1, type='expected'), predict(kfitx, type='expected')) [1] TRUE > aeq(predict(kfit1, type='lp'), predict(kfitx, type='lp')) [1] TRUE > > temp1 <- predict(kfit1, type='terms', se.fit=T) > temp2 <- predict(kfitx, type='terms', se.fit=T) > aeq(temp1$fit[,1:2], temp2$fit) [1] TRUE > # 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)) [1] TRUE > mean(temp1$se.fit[,1:2]/ temp2$se.fit) [1] 1.433017 > aeq(as.vector(temp1$se.fit[,3])^2, + as.vector(kfit1$fvar[match(kidney$id, sort(unique(kidney$id)))])) [1] TRUE > > print(temp1) $fit age sex frailty(id, dist = "gauss") 1 -0.073981042 1.039553 0.59814468 2 -0.073981042 1.039553 0.59814468 3 0.020278123 -0.371269 0.38512389 4 0.020278123 -0.371269 0.38512389 5 -0.055129209 1.039553 0.20210998 6 -0.055129209 1.039553 0.20210998 7 -0.059842167 -0.371269 -0.55932015 8 -0.055129209 -0.371269 -0.55932015 9 -0.158814289 1.039553 0.28558387 10 -0.158814289 1.039553 0.28558387 11 -0.130536540 -0.371269 0.06628942 12 -0.125823582 -0.371269 0.06628942 13 0.034416998 1.039553 0.80505119 14 0.034416998 1.039553 0.80505119 15 0.053268830 -0.371269 -0.43834241 16 0.057981789 -0.371269 -0.43834241 17 0.119250245 -0.371269 -0.05631649 18 0.119250245 -0.371269 -0.05631649 19 0.034416998 1.039553 -0.49980572 20 0.039129956 1.039553 -0.49980572 21 0.001426290 -0.371269 -0.13028264 22 0.001426290 -0.371269 -0.13028264 23 -0.045703292 -0.371269 0.06377401 24 -0.045703292 -0.371269 0.06377401 25 -0.040990334 -0.371269 0.38815296 26 -0.040990334 -0.371269 0.38815296 27 -0.007999626 -0.371269 -0.47650510 28 -0.007999626 -0.371269 -0.47650510 29 -0.125823582 -0.371269 -0.66986830 30 -0.125823582 -0.371269 -0.66986830 31 0.076833621 1.039553 0.19359678 32 0.076833621 1.039553 0.19359678 33 0.076833621 -0.371269 -0.16483200 34 0.076833621 -0.371269 -0.16483200 35 -0.003286668 -0.371269 -0.15794998 36 0.001426290 -0.371269 -0.15794998 37 0.043842914 -0.371269 -0.46236014 38 0.043842914 -0.371269 -0.46236014 39 0.001426290 -0.371269 0.12603308 40 0.001426290 -0.371269 0.12603308 41 0.010852206 1.039553 -1.74303142 42 0.015565165 1.039553 -1.74303142 43 -0.064555125 -0.371269 -0.45211210 44 -0.064555125 -0.371269 -0.45211210 45 0.086259538 -0.371269 0.51574106 46 0.090972496 -0.371269 0.51574106 47 -0.007999626 -0.371269 0.09475123 48 -0.003286668 -0.371269 0.09475123 49 -0.003286668 1.039553 0.05790354 50 -0.003286668 1.039553 0.05790354 51 0.062694747 -0.371269 -0.37933234 52 0.067407705 -0.371269 -0.37933234 53 -0.158814289 -0.371269 0.11248891 54 -0.158814289 -0.371269 0.11248891 55 0.039129956 -0.371269 0.54791210 56 0.039129956 -0.371269 0.54791210 57 0.043842914 1.039553 0.45873482 58 0.043842914 1.039553 0.45873482 59 0.048555872 -0.371269 0.35639797 60 0.048555872 -0.371269 0.35639797 61 0.057981789 -0.371269 0.48803342 62 0.057981789 -0.371269 0.48803342 63 0.029704039 -0.371269 0.25597325 64 0.034416998 -0.371269 0.25597325 65 0.062694747 -0.371269 0.23054948 66 0.062694747 -0.371269 0.23054948 67 0.001426290 -0.371269 -0.13680005 68 0.006139248 -0.371269 -0.13680005 69 -0.102258791 -0.371269 0.51977995 70 -0.102258791 -0.371269 0.51977995 71 -0.007999626 -0.371269 -0.23878154 72 -0.007999626 -0.371269 -0.23878154 73 0.039129956 -0.371269 0.17174306 74 0.039129956 -0.371269 0.17174306 75 0.076833621 1.039553 -0.35822829 76 0.076833621 1.039553 -0.35822829 $se.fit age sex frailty(id, dist = "gauss") 1 0.195861829 0.3280279 0.6246430 2 0.195861829 0.3280279 0.6246430 3 0.053685514 0.1171528 0.6954922 4 0.053685514 0.1171528 0.6954922 5 0.145952360 0.3280279 0.5705340 6 0.145952360 0.3280279 0.5705340 7 0.158429727 0.1171528 0.4894541 8 0.145952360 0.1171528 0.4894541 9 0.420454437 0.3280279 0.6071455 10 0.420454437 0.3280279 0.6071455 11 0.345590234 0.1171528 0.5633997 12 0.333112867 0.1171528 0.5633997 13 0.091117615 0.3280279 0.6641707 14 0.091117615 0.3280279 0.6641707 15 0.141027084 0.1171528 0.5101890 16 0.153504451 0.1171528 0.5101890 17 0.315710223 0.1171528 0.5491569 18 0.315710223 0.1171528 0.5491569 19 0.091117615 0.3280279 0.5264083 20 0.103594982 0.3280279 0.5264083 21 0.003776045 0.1171528 0.5180953 22 0.003776045 0.1171528 0.5180953 23 0.120997626 0.1171528 0.6208806 24 0.120997626 0.1171528 0.6208806 25 0.108520259 0.1171528 0.5811421 26 0.108520259 0.1171528 0.5811421 27 0.021178689 0.1171528 0.6247779 28 0.021178689 0.1171528 0.6247779 29 0.333112867 0.1171528 0.5615987 30 0.333112867 0.1171528 0.5615987 31 0.203413919 0.3280279 0.6532405 32 0.203413919 0.3280279 0.6532405 33 0.203413919 0.1171528 0.5247227 34 0.203413919 0.1171528 0.5247227 35 0.008701322 0.1171528 0.5106606 36 0.003776045 0.1171528 0.5106606 37 0.116072349 0.1171528 0.6284328 38 0.116072349 0.1171528 0.6284328 39 0.003776045 0.1171528 0.6320009 40 0.003776045 0.1171528 0.6320009 41 0.028730780 0.3280279 0.5235228 42 0.041208147 0.3280279 0.5235228 43 0.170907094 0.1171528 0.5492095 44 0.170907094 0.1171528 0.5492095 45 0.228368654 0.1171528 0.6058686 46 0.240846021 0.1171528 0.6058686 47 0.021178689 0.1171528 0.6267998 48 0.008701322 0.1171528 0.6267998 49 0.008701322 0.3280279 0.5526664 50 0.008701322 0.3280279 0.5526664 51 0.165981818 0.1171528 0.5556706 52 0.178459185 0.1171528 0.5556706 53 0.420454437 0.1171528 0.5849825 54 0.420454437 0.1171528 0.5849825 55 0.103594982 0.1171528 0.6081780 56 0.103594982 0.1171528 0.6081780 57 0.116072349 0.3280279 0.6010279 58 0.116072349 0.3280279 0.6010279 59 0.128549717 0.1171528 0.5762113 60 0.128549717 0.1171528 0.5762113 61 0.153504451 0.1171528 0.5982501 62 0.153504451 0.1171528 0.5982501 63 0.078640248 0.1171528 0.6614053 64 0.091117615 0.1171528 0.6614053 65 0.165981818 0.1171528 0.5609510 66 0.165981818 0.1171528 0.5609510 67 0.003776045 0.1171528 0.5844921 68 0.016253412 0.1171528 0.5844921 69 0.270726031 0.1171528 0.6089631 70 0.270726031 0.1171528 0.6089631 71 0.021178689 0.1171528 0.6795741 72 0.021178689 0.1171528 0.6795741 73 0.103594982 0.1171528 0.6421784 74 0.103594982 0.1171528 0.6421784 75 0.203413919 0.3280279 0.5779661 76 0.203413919 0.3280279 0.5779661 > kfit1 Call: coxph(formula = Surv(time, status) ~ age + sex + frailty(id, dist = "gauss"), data = kidney) coef se(coef) se2 Chisq DF p age 0.00471 0.01248 0.00856 0.14267 1.0 0.7056 sex -1.41082 0.44518 0.31504 10.04319 1.0 0.0015 frailty(id, dist = "gauss 26.54461 14.7 0.0294 Iterations: 6 outer, 39 Newton-Raphson Variance of random effect= 0.569 Degrees of freedom for terms= 0.5 0.5 14.7 Likelihood ratio test=47.5 on 15.7 df, p=5e-05 n= 76, number of events= 58 > kfitx Call: coxph(formula = Surv(time, status) ~ age + sex + offset(temp), data = kidney, init = kfit1$coef, iter = 0) coef exp(coef) se(coef) z p age 0.004713 1.004724 0.008749 0.539 0.59 sex -1.410822 0.243943 0.309164 -4.563 5.03e-06 Likelihood ratio test=0 on 2 df, p=1 n= 76, number of events= 58 > > 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)) [1] TRUE > > # Now fit a model with explicit offset > kfitx <- coxph(Surv(time, status) ~ offset(tempf),kidney, eps=1e-7) > > aeq(resid(kfit1), resid(kfitx)) [1] TRUE > aeq(resid(kfit1, type='deviance'), resid(kfitx, type='deviance')) [1] TRUE > > # > # 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')) [1] TRUE > aeq(predict(kfit1, type='lp'), predict(kfitx, type='lp')) [1] TRUE > > temp1 <- predict(kfit1, type='terms', se.fit=T) > aeq(temp1$fit, kfitx$linear) [1] TRUE > aeq(temp1$se.fit^2, + kfit1$fvar[match(kidney$id, sort(unique(kidney$id)))]) [1] TRUE > > temp1 $fit [1] 0.696003729 0.696003729 0.244575316 0.244575316 0.494175549 [6] 0.494175549 -0.659248798 -0.659248798 0.521423106 0.521423106 [11] -0.114492938 -0.114492938 0.800127481 0.800127481 -0.488101282 [16] -0.488101282 -0.120396647 -0.120396647 0.131121515 0.131121515 [21] -0.214987009 -0.214987009 -0.054872789 -0.054872789 0.184657295 [26] 0.184657295 -0.510007747 -0.510007747 -0.790746805 -0.790746805 [31] 0.324674289 0.324674289 -0.239374060 -0.239374060 -0.264428564 [36] -0.264428564 -0.472698773 -0.472698773 0.006304049 0.006304049 [41] -0.873434085 -0.873434085 -0.530880840 -0.530880840 0.351411783 [46] 0.351411783 -0.037212138 -0.037212138 0.442049266 0.442049266 [51] -0.419206550 -0.419206550 -0.108012854 -0.108012854 0.346332076 [56] 0.346332076 0.659300205 0.659300205 0.197278585 0.197278585 [61] 0.304868889 0.304868889 0.139712997 0.139712997 0.093574024 [66] 0.093574024 -0.209690355 -0.209690355 0.302070834 0.302070834 [71] -0.278962288 -0.278962288 0.068599919 0.068599919 0.078493616 [76] 0.078493616 $se.fit [1] 0.6150025 0.6150025 0.6160184 0.6160184 0.5715622 0.5715622 0.4393615 [8] 0.4393615 0.5761369 0.5761369 0.4834244 0.4834244 0.6421184 0.6421184 [15] 0.4574824 0.4574824 0.4813578 0.4813578 0.5119792 0.5119792 0.4764145 [22] 0.4764145 0.5532477 0.5532477 0.5195437 0.5195437 0.5534327 0.5534327 [29] 0.4775572 0.4775572 0.6364522 0.6364522 0.4708988 0.4708988 0.4670896 [36] 0.4670896 0.5600672 0.5600672 0.5641880 0.5641880 0.4650576 0.4650576 [43] 0.4904715 0.4904715 0.5448430 0.5448430 0.5570120 0.5570120 0.5608187 [50] 0.5608187 0.4996021 0.4996021 0.4831697 0.4831697 0.5452255 0.5452255 [57] 0.6057428 0.6057428 0.5209402 0.5209402 0.5376594 0.5376594 0.5911350 [64] 0.5911350 0.5065368 0.5065368 0.5290283 0.5290283 0.5368433 0.5368433 [71] 0.5996077 0.5996077 0.5762814 0.5762814 0.5782753 0.5782753 > kfit1 Call: coxph(formula = Surv(time, status) ~ frailty(id, dist = "gauss"), data = kidney) coef se(coef) se2 Chisq DF p frailty(id, dist = "gauss 23 13.8 0.057 Iterations: 7 outer, 39 Newton-Raphson Variance of random effect= 0.458 Degrees of freedom for terms= 13.8 Likelihood ratio test=33.4 on 13.8 df, p=0.002 n= 76, number of events= 58 > > > > proc.time() user system elapsed 0.90 0.14 1.03