ns-rse
commented
1 month ago I've been doing some more investigations and have found that there are potential issues with missing data for some of the variables that may be causing problems.
I've used the following script based on logistic regression models that @mdp21oe had written
NB - In order to be comparable with the recipe/workflow that we are investigating I have changed the data being used to train instead of the full raw (but tidy) dataset df. This later (i.e. df) includes 214 individuals who are missing final_pathology and can not be included in the analysis anyway since there is no assessment of their thyroid cancer status. Further we want to work out why the various Tidymodelling tuning doesn't work and that uses the train subset so we need to compare like for like.
Setup
Load libraries and data, remove those with missing final_pathology and split into train and test.
library(dplyr)
library(Hmisc)
library(knitr)
library(readr)
library(rmarkdown)
## Libraries for Tidymodelling
library(dials)
library(kernlab)
library(knitr)
library(tidymodels)
library(tidyverse)
library(vip)
## Load and retain only those with final_pathology
df <- readRDS(paste(r_dir, "clean.rds", sep = "/"))
df_complete <- df |>
dplyr::select(
age_at_scan,
gender,
ethnicity,
incidental_nodule,
palpable_nodule,
rapid_enlargement,
compressive_symptoms,
hypertension,
vocal_cord_paresis,
graves_disease,
hashimotos_thyroiditis,
family_history_thyroid_cancer,
exposure_radiation,
albumin,
tsh_value,
lymphocytes,
monocyte,
bta_u_classification,
solitary_nodule,
size_nodule_mm,
cervical_lymphadenopathy,
thy_classification,
final_pathology) |>
dplyr::filter(!is.na(final_pathology))
## Split data into train/test
split <- rsample::initial_split(df_complete, prop = 0.75)
train <- rsample::training(split)
test <- rsample::testing(split)
## Setup cross-validation with k-folds (n = 10)
cv_folds <- rsample::vfold_cv(train, v = 10, repeats = 10)Run glm() on subsets of variables
Clinical variables
## Clinical variables
glm_clin <- glm(final_pathology ~ age_at_scan + gender + incidental_nodule + palpable_nodule +
rapid_enlargement + compressive_symptoms + hashimotos_thyroiditis + family_history_thyroid_cancer + exposure_radiation +
tsh_value + size_nodule_mm + solitary_nodule + cervical_lymphadenopathy,
data = train,
family = binomial(link = "logit"))
summary(glm_clin)
Call:
glm(formula = final_pathology ~ age_at_scan + gender + incidental_nodule +
palpable_nodule + rapid_enlargement + compressive_symptoms +
hashimotos_thyroiditis + family_history_thyroid_cancer +
exposure_radiation + tsh_value + size_nodule_mm + solitary_nodule +
cervical_lymphadenopathy, family = binomial(link = "logit"),
data = train)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.004e+00 7.626e-01 -2.627 0.008608 **
age_at_scan -2.977e-02 1.171e-02 -2.542 0.011026 *
genderM 7.614e-01 4.115e-01 1.850 0.064303 .
incidental_noduleYes -1.308e-01 4.969e-01 -0.263 0.792384
palpable_noduleYes -1.011e-01 5.328e-01 -0.190 0.849506
rapid_enlargementYes 5.313e-03 9.793e-01 0.005 0.995671
compressive_symptomsYes 1.170e-01 5.556e-01 0.211 0.833266
hashimotos_thyroiditisYes -1.459e+01 2.400e+03 -0.006 0.995149
family_history_thyroid_cancerYes 1.701e+01 2.400e+03 0.007 0.994343
exposure_radiationYes -1.489e+01 1.671e+03 -0.009 0.992889
tsh_value 1.336e-01 7.379e-02 1.811 0.070200 .
size_nodule_mm 4.478e-02 1.174e-02 3.816 0.000135 ***
solitary_noduleYes 6.990e-01 3.594e-01 1.945 0.051776 .
cervical_lymphadenopathyYes 3.017e+00 9.198e-01 3.280 0.001038 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 286.13 on 325 degrees of freedom
Residual deviance: 223.54 on 312 degrees of freedom
(536 observations deleted due to missingness)
AIC: 251.54
Number of Fisher Scoring iterations: 15
glm_clin
Call: glm(formula = final_pathology ~ age_at_scan + gender + incidental_nodule +
palpable_nodule + rapid_enlargement + compressive_symptoms +
hashimotos_thyroiditis + family_history_thyroid_cancer +
exposure_radiation + tsh_value + size_nodule_mm + solitary_nodule +
cervical_lymphadenopathy + incidental_nodule + tsh_value +
size_nodule_mm + solitary_nodule + bta_u_classification,
family = binomial(link = "logit"), data = train)
Coefficients:
(Intercept) age_at_scan genderM
-2.98963 -0.02312 0.61605
incidental_noduleYes palpable_noduleYes rapid_enlargementYes
-0.43710 -0.44385 -0.18599
compressive_symptomsYes hashimotos_thyroiditisYes family_history_thyroid_cancerYes
0.22305 -16.28615 16.10881
exposure_radiationYes tsh_value size_nodule_mm
-17.31802 0.15792 0.03851
solitary_noduleYes cervical_lymphadenopathyYes bta_u_classificationU3
0.77457 2.03330 1.94868
bta_u_classificationU4 bta_u_classificationU5
3.12813 17.74532
Degrees of Freedom: 312 Total (i.e. Null); 296 Residual
(549 observations deleted due to missingness) <<<<<<<<< IMPORTANT!!!
Null Deviance: 275
Residual Deviance: 180.3 AIC: 214.3549 observations dropped due to missing data.
Biochemical markers
glm_biochem <- glm(final_pathology ~ age_at_scan + gender + tsh_value + albumin + lymphocytes + monocyte,
data = train,
family = binomial(link = "logit")
)
summary(glm_biochem)
Call:
glm(formula = final_pathology ~ age_at_scan + gender + tsh_value +
albumin + lymphocytes + monocyte, family = binomial(link = "logit"),
data = train)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.73160 2.67009 -0.649 0.516652
age_at_scan -0.03118 0.01122 -2.780 0.005441 **
genderM 1.34973 0.39931 3.380 0.000725 ***
tsh_value 0.11219 0.07189 1.561 0.118623
albumin 0.04716 0.05323 0.886 0.375708
lymphocytes 0.27508 0.26067 1.055 0.291292
monocyte -3.82506 1.27730 -2.995 0.002748 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 264.17 on 347 degrees of freedom
Residual deviance: 227.00 on 341 degrees of freedom
(514 observations deleted due to missingness)
AIC: 241
Number of Fisher Scoring iterations: 6
glm_biochem
Call: glm(formula = final_pathology ~ age_at_scan + gender + tsh_value +
albumin + lymphocytes + monocyte, family = binomial(link = "logit"),
data = train)
Coefficients:
(Intercept) age_at_scan genderM tsh_value albumin lymphocytes monocyte
-1.73160 -0.03118 1.34973 0.11219 0.04716 0.27508 -3.82506
Degrees of Freedom: 347 Total (i.e. Null); 341 Residual
(514 observations deleted due to missingness) <<<<<<<<< IMPORTANT!!!
Null Deviance: 264.2
Residual Deviance: 227 AIC: 241514 observations deleted due to missing data.
Ultrasound features
glm_ultrasound <- glm(
final_pathology ~ age_at_scan + gender + incidental_nodule + tsh_value + size_nodule_mm +
solitary_nodule + cervical_lymphadenopathy + bta_u_classification, # + thy_classification,
data = train,
family = binomial(link = "logit")
)
summary(glm_ultrasound)
Call:
glm(formula = final_pathology ~ age_at_scan + gender + incidental_nodule +
tsh_value + size_nodule_mm + solitary_nodule + cervical_lymphadenopathy +
bta_u_classification, family = binomial(link = "logit"),
data = train)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.92562 0.67664 -4.324 1.53e-05 ***
age_at_scan -0.02624 0.01140 -2.302 0.02135 *
genderM 0.62732 0.40513 1.548 0.12152
incidental_noduleYes -0.35053 0.38496 -0.911 0.36252
tsh_value 0.14586 0.07193 2.028 0.04257 *
size_nodule_mm 0.03141 0.01068 2.940 0.00328 **
solitary_noduleYes 0.32866 0.36267 0.906 0.36482
cervical_lymphadenopathyYes 0.64249 0.88127 0.729 0.46597
bta_u_classificationU3 1.98761 0.38883 5.112 3.19e-07 ***
bta_u_classificationU4 3.15909 0.69606 4.539 5.67e-06 ***
bta_u_classificationU5 18.13639 1164.29414 0.016 0.98757
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 329.11 on 404 degrees of freedom
Residual deviance: 225.02 on 394 degrees of freedom
(457 observations deleted due to missingness)
AIC: 247.02
Number of Fisher Scoring iterations: 15
glm_ultrasound
Call: glm(formula = final_pathology ~ age_at_scan + gender + incidental_nodule +
tsh_value + size_nodule_mm + solitary_nodule + cervical_lymphadenopathy +
bta_u_classification, family = binomial(link = "logit"),
data = train)
Coefficients:
(Intercept) age_at_scan genderM incidental_noduleYes
-2.92562 -0.02624 0.62732 -0.35053
tsh_value size_nodule_mm solitary_noduleYes cervical_lymphadenopathyYes
0.14586 0.03141 0.32866 0.64249
bta_u_classificationU3 bta_u_classificationU4 bta_u_classificationU5
1.98761 3.15909 18.13639
Degrees of Freedom: 404 Total (i.e. Null); 394 Residual
(457 observations deleted due to missingness) <<<<<<<<< IMPORTANT!!!
Null Deviance: 329.1
Residual Deviance: 225 AIC: 247457 observations dropped due to missing data.
Saturated model
In the Tidymodeling we are attempting to use a saturated model with all variables because as can be seen from above the coefficients vary wildly depending on the other co-variates. If we look at the values for age_at_scan and gender which are included in all and they are bouncing around all over the place. What we want is a fully adjusted estimate for each co-variate taking into account all the other variables
glm_all <- glm(
final_pathology ~ age_at_scan + gender + incidental_nodule + palpable_nodule +
rapid_enlargement + compressive_symptoms + hashimotos_thyroiditis + family_history_thyroid_cancer + exposure_radiation +
tsh_value +
size_nodule_mm +
solitary_nodule +
cervical_lymphadenopathy +
albumin +
lymphocytes +
monocyte +
incidental_nodule +
tsh_value +
size_nodule_mm +
solitary_nodule +
bta_u_classification,
data = train,
family = binomial(link = "logit")
)
glm_all <- glm(
+ final_pathology ~ age_at_scan +
+ gender +
+ incidental_nodule +
+ palpable_nodule +
+ rapid_enlargement +
+ compressive_symptoms +
+ hashimotos_thyroiditis +
+ family_history_thyroid_cancer +
+ exposure_radiation +
+ tsh_value +
+ size_nodule_mm +
+ solitary_nodule +
+ cervical_lymphadenopathy +
+ albumin +
+ lymphocytes +
+ monocyte +
+ incidental_nodule +
+ tsh_value +
+ size_nodule_mm +
+ solitary_nodule +
+ bta_u_classification,
+ data = train,
+ family = binomial(link = "logit")
+ )
+ + + + + + + + + + + + + + + + + + + + + + + + Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levelsOh no our model didn't fit we've got a weird error cropping up here?!?!?
Looking at the pattern of missing data (see above graph and tables in current rendering of index.qmd) we see the variables with the most missing data are...
| Variable | N | % |
|---|---|---|
thy_classification |
915 | 79.6 |
albumin |
515 | 44.8 |
tsh_value |
413 | 35.9 |
monocyte |
363 | 31.6 |
lymphocytes |
359 | 31.2 |
size_nodule_mm |
319 | 27.7 |
family_history_thyroid_cancer |
281 | 24.4 |
We have already excluded thy_classification in the above but there are some others with very high levels of missing data, namely the biochemical markers obtained from blood assays albumin, lymphocytes, monocyte and tsh_value.
If individuals are excluded based on these and it leaves a individuals who for one of the other binary variables are invariant, e.g. exposure_radiation where only 6 individuals are Yes and the rest are No then the glm() model will complain with the above message that contrasts can be applied only to factors with 2 or more levels because we have to have variation in our predictor for it to be used in the model (if the subset that remains after dropping out those with missing across all factors are all exposure_radiation = No then we will see this sort of error, although it's worth noting I've not investigated whether this is the problem variable).
But the biochemical markers are also problematic because they have so much missing data. If we exclude these biochemical markers and the exposure_radiation too for good measure as there are so few people with Yes we can fit a saturated model.
## Saturated model
glm_all <- glm(
final_pathology ~ age_at_scan +
gender +
incidental_nodule +
palpable_nodule +
rapid_enlargement +
compressive_symptoms +
hashimotos_thyroiditis +
family_history_thyroid_cancer +
## exposure_radiation + <<< We can easily exclude variables by laying out code in
## tsh_value + <<< this manner and adding comments to start of lines ;-)
size_nodule_mm +
solitary_nodule +
cervical_lymphadenopathy +
## albumin +
## lymphocytes +
## monocyte +
incidental_nodule +
size_nodule_mm +
solitary_nodule +
bta_u_classification,
data = train,
family = binomial(link = "logit")
)
summary(glm_all)
Call:
glm(formula = final_pathology ~ age_at_scan + gender + incidental_nodule +
palpable_nodule + rapid_enlargement + compressive_symptoms +
hashimotos_thyroiditis + family_history_thyroid_cancer +
size_nodule_mm + solitary_nodule + cervical_lymphadenopathy +
incidental_nodule + size_nodule_mm + solitary_nodule + bta_u_classification,
family = binomial(link = "logit"), data = train)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.75668 0.75789 -3.637 0.000276 ***
age_at_scan -0.02455 0.01143 -2.147 0.031782 *
genderM 0.49636 0.40278 1.232 0.217822
incidental_noduleYes -0.61578 0.49288 -1.249 0.211540
palpable_noduleYes -0.18385 0.51993 -0.354 0.723637
rapid_enlargementYes -0.27294 1.08007 -0.253 0.800497
compressive_symptomsYes -0.12914 0.57947 -0.223 0.823640
hashimotos_thyroiditisYes -15.60853 2399.54500 -0.007 0.994810
family_history_thyroid_cancerYes 0.82812 2.04767 0.404 0.685904
size_nodule_mm 0.03477 0.01109 3.136 0.001711 **
solitary_noduleYes 0.67604 0.36639 1.845 0.065018 .
cervical_lymphadenopathyYes 1.32420 0.98795 1.340 0.180132
bta_u_classificationU3 2.16334 0.39294 5.506 3.68e-08 ***
bta_u_classificationU4 3.41627 0.75777 4.508 6.53e-06 ***
bta_u_classificationU5 18.92801 954.07210 0.020 0.984172
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 338.70 on 437 degrees of freedom
Residual deviance: 219.56 on 423 degrees of freedom
(424 observations deleted due to missingness)
AIC: 249.56
Number of Fisher Scoring iterations: 15
glm_all
>
Call: glm(formula = final_pathology ~ age_at_scan + gender + incidental_nodule +
palpable_nodule + rapid_enlargement + compressive_symptoms +
hashimotos_thyroiditis + family_history_thyroid_cancer +
size_nodule_mm + solitary_nodule + cervical_lymphadenopathy +
incidental_nodule + size_nodule_mm + solitary_nodule + bta_u_classification,
family = binomial(link = "logit"), data = train)
Coefficients:
(Intercept) age_at_scan genderM
-2.75668 -0.02455 0.49636
incidental_noduleYes palpable_noduleYes rapid_enlargementYes
-0.61578 -0.18385 -0.27294
compressive_symptomsYes hashimotos_thyroiditisYes family_history_thyroid_cancerYes
-0.12914 -15.60853 0.82812
size_nodule_mm solitary_noduleYes cervical_lymphadenopathyYes
0.03477 0.67604 1.32420
bta_u_classificationU3 bta_u_classificationU4 bta_u_classificationU5
2.16334 3.41627 18.92801
Degrees of Freedom: 437 Total (i.e. Null); 423 Residual
(424 observations deleted due to missingness)
Null Deviance: 338.7
Residual Deviance: 219.6 AIC: 249.6424 observations are still dropped due to missing data though!!!
This hasn't resolved all of the problems but is something to bear in mind.
I've asked a question on the Data Science Learning Slack channel (specifically in the #help-r-model channel, but the first link in this paragraph will take you there, although you will need to sign up to the Slack channel first @mdp21oe).
In the current workflow where we select all variables to be included we find that on fitting models (
glm()/glmnet()/ etc.) onlyage_at_scanandgenderare retained as variables.Why is this?
Simple model
Start by taking a simplified model with less variables and test it with the Tidymodels approach to running logistic regression
We see that 4 observations are dropped due to "missingness". Lets add a few more variables in to our predictors.
Weird, neither
incidental_nodulenorpalpable_nodulehave coefficients returned, lets checkpalpable_noduleon its own.There is an effect there for
palpable_nodulealthough 43 observations are dropped (NB no assessment of significance is made at this point, confidence intervals may be very wide around that point estimate).Missing data patterns
Lets check the patterns of missing data using visdat
This is a useful plot, it shows the pattern of missing data across the selected variables and of particular note is that
thyroid_classificationis missing for almost everyone, we will have to do one of two things, either exclude this variable or derive a category for those with missing. I'm not sure but suspect this is only available for those who have had biopsies and/or surgery.