# Is this a correct way to confirm Accelerated Failure Time (AFT) method's suitability for analysis? An example

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)

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


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"))