biometrician / abe

An R package for Augmented Backward Elimination
GNU General Public License v3.0
3 stars 0 forks source link

#47: problems with coxph objects computed by ´abe´ #52

Open biometrician opened 1 year ago

biometrician commented 1 year ago

Hi Rok,

Working with the package I found another issue, which would be really nice, if we could fix it for the paper.

When working with Cox models typically I would like to check the proportional hazards assumption using the cox.zph function. For some model objects computed with the ´abe´ function this works, but unfortunately not for all. Since we state that these objects are of class ´coxph´ it would be nice if standard post-precessing functions like cox.zph also work.

Here is some code with different combinations (only continuous variables, adding factors, adding robust SEs). Concerning the robust SEs, in my examples with abe these work just fine with this exception.

pacman::p_load(abe, shrink, survival)

data(deepvein)

levels(deepvein$sex) <- paste(".", levels(deepvein$sex), sep = "")
levels(deepvein$loc) <- paste(".", levels(deepvein$loc), sep = "")
levels(deepvein$fiimut) <- paste(".", levels(deepvein$fiimut), sep = "")
levels(deepvein$fvleid) <- paste(".", levels(deepvein$fvleid), sep = "")

deepvein$sex <- relevel(deepvein$sex, ref = ".female")
deepvein$fiimut <- relevel(deepvein$fiimut, ref = ".absent")
deepvein$fvleid <- relevel(deepvein$fvleid, ref = ".absent")
deepvein$loc <- relevel(deepvein$loc, ref = ".PE")

#only continuous variables in the model
global_model1 <- coxph(Surv(time, status) ~ age + bmi + log2ddim + durther, 
                  data = deepvein, x = TRUE, y = TRUE)
abe_model <- abe(fit = global_model1,
                        data = deepvein,
                        criterion = "AIC", 
                        tau = 0.05, verbose = FALSE)
# works
cox.zph(abe_model, transform = "identity")

#add factor variables and use type.factor = "factor"
global_model2 <- coxph(Surv(time, status) ~ sex + age + bmi + log2ddim + durther + loc + 
                        fvleid + fiimut, 
                  data = deepvein, x = TRUE, y = TRUE)
abe_factor_model <- abe(fit = global_model2,
                        data = deepvein,
                        criterion = "AIC", 
                        tau = 0.05,
                        type.factor = "factor",
                        exact = TRUE, verbose = FALSE)
# DOES NOT WORK
cox.zph(abe_factor_model, transform = "identity")

#add factor variables and use type.factor = "individual"
deepvein$loc_proximal <- c(deepvein$loc == ".proximal") * 1
deepvein$loc_distal <- c(deepvein$loc == ".distal") * 1

global_dummy_model <- coxph(Surv(time, status) ~ sex + age + bmi + log2ddim + 
    durther + loc_proximal + loc_distal + fvleid + fiimut, 
                            data = deepvein, 
                            x = TRUE, 
                            y = TRUE)

abe_individual_model <- abe(global_dummy_model,
                            data = deepvein,
                            criterion ="AIC", 
                            tau = 0.05, type.factor = "individual", verbose = FALSE)
# This works!
cox.zph(abe_individual_model, transform = "identity")

#add factor variables and use type.factor = "individual" and request robust SE
global_dummy_model2 <- coxph(Surv(time, status) ~ cluster(pnr) + sex + age + bmi + log2ddim + 
    durther + loc_proximal + loc_distal + fvleid + fiimut, 
                            data = deepvein, 
                            x = TRUE, 
                            y = TRUE)

abe_individual_model2 <- abe(global_dummy_model2,
                            data = deepvein,
                            criterion ="AIC", 
                            tau = 0.05, type.factor = "individual", verbose = FALSE)
# DOES NOT WORK, but different error
cox.zph(abe_individual_model2, transform = "identity")

Do you think it make sense simply to refit the models when factors are included in the model to acieve a "clean" fit?

Thanks a lot! Daniela

rokblagus commented 1 year ago

This is all due to us having to use my_updateX functions. I think simply refitting the model should solve all the issues (including us using two different versions of the my_update functions. I tried a very simple fix by evaluating again the final call (its ugly but it kinda works [for the above examples at least]). Can you check if everything still works?

biometrician commented 1 year ago

Hi Rok,

Thanks for the quick fix. The first errors from the example above are now solved. Only the last error still exists. There the error message was different. It is the same error message we had before.

pacman::p_load(abe, shrink, survival)

data(deepvein)
levels(deepvein$sex) <- paste(".", levels(deepvein$sex), sep = "")
levels(deepvein$loc) <- paste(".", levels(deepvein$loc), sep = "")
levels(deepvein$fiimut) <- paste(".", levels(deepvein$fiimut), sep = "")
levels(deepvein$fvleid) <- paste(".", levels(deepvein$fvleid), sep = "")

deepvein$sex <- relevel(deepvein$sex, ref = ".female")
deepvein$fiimut <- relevel(deepvein$fiimut, ref = ".absent")
deepvein$fvleid <- relevel(deepvein$fvleid, ref = ".absent")
deepvein$loc <- relevel(deepvein$loc, ref = ".PE")

#add factor variables and use type.factor = "individual" and request robust SE
global_dummy_model2 <- coxph(Surv(time, status) ~ cluster(pnr) + sex + age + bmi + log2ddim + 
    durther + loc_proximal + loc_distal + fvleid + fiimut, 
                            data = deepvein, 
                            x = TRUE, 
                            y = TRUE)

abe_individual_model2 <- abe(global_dummy_model2,
                            data = deepvein,
                            criterion ="AIC", 
                            tau = 0.05, type.factor = "individual", verbose = FALSE)
# DOES NOT WORK, but different error 
cox.zph(abe_individual_model2, transform = "identity")

When you have time, can you have a look.

Thanks a lot! Best, Daniela

rokblagus commented 1 year ago

Here I already get the error with the coxph function. If I omit variables loc_proximal and loc_distal in the call to coxph, the rest of the code works without any problem. the code also works if I add deepvein$loc_proximal <- c(deepvein$loc == ".proximal") 1 deepvein$loc_distal <- c(deepvein$loc == ".distal") 1 and fit the cox model as specified above.

biometrician commented 1 year ago

Sorry, it seems I deleted too much of my code including the two lines defining the loc_ variables.

Here is the complete code.

pacman::p_load(abe, shrink, survival)

data(deepvein)
levels(deepvein$sex) <- paste(".", levels(deepvein$sex), sep = "")
levels(deepvein$loc) <- paste(".", levels(deepvein$loc), sep = "")
levels(deepvein$fiimut) <- paste(".", levels(deepvein$fiimut), sep = "")
levels(deepvein$fvleid) <- paste(".", levels(deepvein$fvleid), sep = "")

deepvein$sex <- relevel(deepvein$sex, ref = ".female")
deepvein$fiimut <- relevel(deepvein$fiimut, ref = ".absent")
deepvein$fvleid <- relevel(deepvein$fvleid, ref = ".absent")
deepvein$loc <- relevel(deepvein$loc, ref = ".PE")

deepvein$loc_proximal <- c(deepvein$loc == ".proximal") * 1
deepvein$loc_distal <- c(deepvein$loc == ".distal") * 1

#add factor variables and use type.factor = "individual" and request robust SE
global_dummy_model2 <- coxph(Surv(time, status) ~ cluster(pnr) + sex + age + bmi + log2ddim + 
                               durther + loc_proximal + loc_distal + fvleid + fiimut, 
                             data = deepvein, 
                             x = TRUE, 
                             y = TRUE)
# works
cox.zph(global_dummy_model2, transform = "identity")

abe_individual_model2 <- abe(global_dummy_model2,
                             data = deepvein,
                             criterion ="AIC", 
                             tau = 0.05, type.factor = "individual", verbose = FALSE)
# DOES NOT WORK
cox.zph(abe_individual_model2, transform = "identity")

The second call to cox.zph with abe_individual_model2 leads to this error:

Error in model.frame.default(formula = Surv(time, status) ~ sex.male + : 'data' must be a data.frame, environment, or list

Do you have an idea how we can fix this? Thanks, Daniela

rokblagus commented 1 year ago

If I run your code, everything works fine. Are you sure you are using the final version of the function abe? I am using the functions in the branch fixes.