4
$\begingroup$

I need to conduct an adjusted survival analysis, however, Cox PH assumption was not met and data stratification isn't a good solution as one of my study groups already has a small number of patients. Thus, I tried to fit my first AFT model, however, I am not sure if everything is fine. Is the following approach OK for answering my study question?

Study question: to compared adjusted survival of two patient groups, A (n=11000) and B (n=50)

Follow up time: 460 days

Variables used for adjusting: age (50-100), sex (m/f), comorbidity(integer, 0-12)

The approach I used (my current best understanding of AFT)

1. I run AFT models with different distributions (here I only show two models to reduce the length of this post)

library(flexsurv)
fit_genf = flexsurvreg(Surv(time, status) ~ group, dist="genf", data = data)
fit_llogis = flexsurvreg(Surv(time, status) ~ group, dist="llogis", data = data)

2. Then I plot the models to assess their fit with the data and choose the best

plot(fit_genf)

Plot: Generalised F distribution AFT model

plot(fit_llogis)

Plot: Logistic scale distribution AFT model

3. Then I check their log-likelihoods and choose the best

fit_genf$loglik

-29540.38

fit_llogis$loglik

-29748.55

4. Then I check their AIC values and choose the best

fit_genf$AIC

59090.77

fit_llogis$AIC

59503.1

5. RESULT: Generalised F distribution AFT modelling (fit_genf) is an acceptable method for my analysis since this had the best fit, better log-likelihood and smaller AIC. Or is something wrong in my approach (large AIC number, log-likelihood etc)?

Full fit_genf model summary:

Call:
flexsurvreg(formula = Surv(time, status) ~ group, data = data, 
    dist = "genf")

Estimates: 
        data mean  est       L95%      U95%      se        exp(est)
mu            NA    4.82150   4.29511   5.34788   0.26857        NA
sigma         NA    2.37833   1.89590   2.98351   0.27509        NA
Q             NA   -2.79284  -3.71926  -1.86641   0.47267        NA
P             NA    2.22999   1.05947   4.69368   0.84675        NA
groupB   0.00433    0.62147  -0.29332   1.53627   0.46674   1.86167
        L95%      U95%    
mu            NA        NA
sigma         NA        NA
Q             NA        NA
P             NA        NA
groupB   0.74578   4.64721

N = 11541,  Events: 3851,  Censored: 7690
Total time at risk: 4024558
Log-likelihood = -29540.38, df = 5
AIC = 59090.77

ADDED AFTER CARLO'S RESPONSE

6. COMPLETE MODEL, named "final"

final = flexsurvreg(formula = Surv(time, status) ~ group + sex + 
        age + comorbidity, data = data, dist = "gompertz")
final


Call:
flexsurvreg(formula = Surv(time, status) ~ group + sex + 
    age + comorbidity, data = data, dist = "gompertz")

Estimates: 
                   data mean   est         L95%        U95%        se          exp(est)    L95%        U95%      
shape                      NA  -0.1679218  -0.1766458  -0.1591978   0.0044511          NA          NA          NA
rate                       NA   0.0005775   0.0004209   0.0007923   0.0000932          NA          NA          NA
groupB              0.0043324  -0.2059110  -0.7134321   0.3016101   0.2589441   0.8139055   0.4899597   1.3520339
sexFemale           0.7172689  -0.4684015  -0.5421893  -0.3946137   0.0376475   0.6260021   0.5814739   0.6739403
age                78.7702972   0.0615099   0.0576167   0.0654030   0.0019863   1.0634410   1.0593089   1.0675892
comorbidity         1.6603414   0.1387881   0.1209541   0.1566220   0.0090991   1.1488806   1.1285732   1.1695534

N = 11541,  Events: 3851,  Censored: 7690
Total time at risk: 132299.7
Log-likelihood = -15802.47, df = 6
AIC = 31616.93

7. HOW TO CHECK COMPLETE MODEL'S FIT?

I tried to do this as follows using different values for replacing the question marks, but got this error: Error in match.arg(type) : 'arg' must be NULL or a character vector

KaplanMeier = survfit(Surv(time, status) ~ group, data = data)

plot(KaplanMeier)
lines(predict(final, newdata=list(group = "A",
                                 sex = ?, 
                                 age = ?, 
                                 comorbidity_score = ?,
                                 type="quantile",
                                 p=seq(.01,.99,by=.01)),
                                 seq(.99,.01,by=-.01),
                                 col="blue"))
lines(predict(final, newdata=list(group = "B",
                                 sex = ?, 
                                 age = ?, 
                                 comorbidity_score = ?,
                                 type="quantile",
                                 p=seq(.01,.99,by=.01)),
                                 seq(.99,.01,by=-.01),
                                 col="blue"))
$\endgroup$

1 Answer 1

2
$\begingroup$

I know that it's been a while since this question was first posted, but I came across this resource that may help you. Take a look at the slide that's third from the bottom: http://www.ams.sunysb.edu/~zhu/ams588/Lecture_5_AFT.pdf

It seems that you can fit a couple of models, and given that, for example, the exponential distribution is a subset of the Gamma distribution, you can use the log-likelihood values to determine which ones are better fit -- but only for nested models. Smaller log-likelihood values is indicative of a better fit (i.e., closer to 0). However, the slides also state that you assume that the Gamma distribution is reasonable to start with.

There's also this paper here as well: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2376927/pdf/89-6601120a.pdf

Hope this helps someone else who comes across a similar question.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.