This project examines relevant variables for potential collinearity when predicting positive and negative sentiment using clinical notes. We specifically examine race and gender. Other covariates include Charlson scores, age, insurance status, physician specialty, and ethnicity.

We also deploy Logistic Regression, Naive Bayes, Random Forest, and Artificial Neural Networks to predict positive or negative sentiment.

sent_04_23 <- read.csv("/Users/KKW4/Desktop/AIM-AHEAD/Research projects/Sentiment analysis/final_sent_0423.csv")


library(corrplot)
## corrplot 0.92 loaded
library(jtools)
library(mediation)
## Loading required package: MASS
## Loading required package: Matrix
## Loading required package: mvtnorm
## Loading required package: sandwich
## mediation: Causal Mediation Analysis
## Version: 4.5.0
library(e1071)
library(naivebayes)
## naivebayes 0.9.7 loaded
library(randomForest)
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
library(caret)
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:randomForest':
## 
##     margin
## Loading required package: lattice
library(keras)
## 
## Attaching package: 'keras'
## The following object is masked from 'package:jtools':
## 
##     get_weights
library(mlbench)
library(neuralnet)
library(tidyverse)
## Registered S3 methods overwritten by 'broom':
##   method            from  
##   tidy.glht         jtools
##   tidy.summary.glht jtools
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ✔ purrr   0.3.5      
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::combine()  masks randomForest::combine()
## ✖ dplyr::compute()  masks neuralnet::compute()
## ✖ tidyr::expand()   masks Matrix::expand()
## ✖ dplyr::filter()   masks stats::filter()
## ✖ dplyr::lag()      masks stats::lag()
## ✖ purrr::lift()     masks caret::lift()
## ✖ ggplot2::margin() masks randomForest::margin()
## ✖ tidyr::pack()     masks Matrix::pack()
## ✖ dplyr::select()   masks MASS::select()
## ✖ tidyr::unpack()   masks Matrix::unpack()
library(devtools)
## Loading required package: usethis

This model includes all variables. Set is the setting or physician specialty. We used packge jtools to produce clear tables. Charlson was removed given the number of NAs or missings it produced in the regression model. These are all now factor variables and previously were not. That was an error.

#with all variables included - age is a categorical variable here
hin_mod <- glm(as.factor(sent_bin) ~ as.factor(eth) + as.factor(sex) + as.factor(race) + as.factor(age_group) + as.factor(ins.y) + as.factor(set), data = sent_04_23, family = "binomial")

summ(hin_mod)
## MODEL INFO:
## Observations: 1155 (662 missing obs. deleted)
## Dependent Variable: as.factor(sent_bin)
## Type: Generalized linear model
##   Family: binomial 
##   Link function: logit 
## 
## MODEL FIT:
## χ²(15) = 45.24, p = 0.00
## Pseudo-R² (Cragg-Uhler) = 0.05
## Pseudo-R² (McFadden) = 0.03
## AIC = 1529.36, BIC = 1610.19 
## 
## Standard errors: MLE
## ----------------------------------------------------------
##                                Est.   S.E.   z val.      p
## --------------------------- ------- ------ -------- ------
## (Intercept)                   -1.13   0.45    -2.52   0.01
## as.factor(eth)1               -0.18   0.34    -0.52   0.61
## as.factor(eth)2               -1.78   0.78    -2.28   0.02
## as.factor(sex)1               -0.25   0.13    -1.96   0.05
## as.factor(race)2               0.03   0.13     0.26   0.80
## as.factor(race)3              -0.41   0.45    -0.92   0.36
## as.factor(race)4               1.16   0.84     1.38   0.17
## as.factor(race)5              -0.73   0.52    -1.40   0.16
## as.factor(age_group)2          0.74   0.35     2.14   0.03
## as.factor(age_group)3          0.40   0.35     1.12   0.26
## as.factor(age_group)4          0.68   0.38     1.79   0.07
## as.factor(ins.y)2              0.48   0.29     1.67   0.09
## as.factor(ins.y)3              0.56   0.30     1.91   0.06
## as.factor(ins.y)4             -0.43   0.51    -0.84   0.40
## as.factor(ins.y)5              0.48   0.60     0.80   0.42
## as.factor(set)1               -0.23   0.13    -1.80   0.07
## ----------------------------------------------------------

Odds ratios for model 1

exp(cbind(OR = coef(hin_mod), confint(hin_mod)))
## Waiting for profiling to be done...
##                              OR      2.5 %     97.5 %
## (Intercept)           0.3214687 0.12972229  0.7605121
## as.factor(eth)1       0.8390067 0.42382270  1.6184089
## as.factor(eth)2       0.1678659 0.02540262  0.6364133
## as.factor(sex)1       0.7790934 0.60704151  1.0000423
## as.factor(race)2      1.0346355 0.79977547  1.3394809
## as.factor(race)3      0.6615598 0.26124355  1.5430609
## as.factor(race)4      3.1895157 0.62150436 18.6252228
## as.factor(race)5      0.4835674 0.16606506  1.2981347
## as.factor(age_group)2 2.1049895 1.08803657  4.3078544
## as.factor(age_group)3 1.4875850 0.75920944  3.0761767
## as.factor(age_group)4 1.9716999 0.95424831  4.2704518
## as.factor(ins.y)2     1.6227218 0.93193976  2.9094768
## as.factor(ins.y)3     1.7565041 0.99549144  3.1852397
## as.factor(ins.y)4     0.6520059 0.22695007  1.7188203
## as.factor(ins.y)5     1.6218001 0.47565202  5.2252104
## as.factor(set)1       0.7938745 0.61697297  1.0213846

All variables are included in this model and age is numeric. We run interactions with age.

#with all variables included - age is a numeric variable here and interacts with race
hin_mod_1a <- glm(as.factor(sent_bin) ~ as.factor(eth) + as.factor(sex)*Patient_age + as.factor(race)*Patient_age +  as.factor(ins.y) + as.factor(set), data = sent_04_23, family = "binomial")

summ(hin_mod_1a)
## MODEL INFO:
## Observations: 1291 (526 missing obs. deleted)
## Dependent Variable: as.factor(sent_bin)
## Type: Generalized linear model
##   Family: binomial 
##   Link function: logit 
## 
## MODEL FIT:
## χ²(18) = 55.25, p = 0.00
## Pseudo-R² (Cragg-Uhler) = 0.06
## Pseudo-R² (McFadden) = 0.03
## AIC = 1716.00, BIC = 1814.10 
## 
## Standard errors: MLE
## -----------------------------------------------------------------
##                                       Est.   S.E.   z val.      p
## ---------------------------------- ------- ------ -------- ------
## (Intercept)                           0.83   0.69     1.20   0.23
## as.factor(eth)1                      -0.44   0.33    -1.31   0.19
## as.factor(eth)2                      -1.16   0.66    -1.76   0.08
## as.factor(sex)1                      -0.85   0.62    -1.37   0.17
## Patient_age                          -0.02   0.01    -2.29   0.02
## as.factor(race)2                     -0.28   0.59    -0.48   0.63
## as.factor(race)3                     -1.04   2.68    -0.39   0.70
## as.factor(race)4                      3.28   2.64     1.24   0.21
## as.factor(race)5                     -5.09   2.22    -2.30   0.02
## as.factor(ins.y)2                     0.67   0.28     2.38   0.02
## as.factor(ins.y)3                     0.54   0.29     1.87   0.06
## as.factor(ins.y)4                    -0.22   0.48    -0.45   0.65
## as.factor(ins.y)5                     0.70   0.59     1.18   0.24
## as.factor(set)1                      -0.46   0.12    -3.84   0.00
## as.factor(sex)1:Patient_age           0.01   0.01     1.31   0.19
## Patient_age:as.factor(race)2          0.00   0.01     0.35   0.73
## Patient_age:as.factor(race)3          0.01   0.04     0.30   0.77
## Patient_age:as.factor(race)4         -0.05   0.05    -0.99   0.32
## Patient_age:as.factor(race)5          0.07   0.04     2.08   0.04
## -----------------------------------------------------------------

Odds ratios for the model with patient age as numeric and interactions

exp(cbind(OR = coef(hin_mod_1a), confint(hin_mod_1a)))
## Waiting for profiling to be done...
##                                        OR        2.5 %       97.5 %
## (Intercept)                   2.296447396 0.5936889646 8.943993e+00
## as.factor(eth)1               0.646635023 0.3303708920 1.226995e+00
## as.factor(eth)2               0.313196756 0.0699665589 1.011357e+00
## as.factor(sex)1               0.427100956 0.1257541072 1.436228e+00
## Patient_age                   0.975554474 0.9548810540 9.963549e-01
## as.factor(race)2              0.754160481 0.2356301892 2.410962e+00
## as.factor(race)3              0.352827882 0.0013676301 6.660878e+01
## as.factor(race)4             26.631304464 0.1897955446 1.129080e+04
## as.factor(race)5              0.006174171 0.0000395544 2.888678e-01
## as.factor(ins.y)2             1.956135383 1.1394908522 3.458996e+00
## as.factor(ins.y)3             1.713270557 0.9856524909 3.061643e+00
## as.factor(ins.y)4             0.805783735 0.3048271762 2.012925e+00
## as.factor(ins.y)5             2.004604470 0.6121010778 6.368221e+00
## as.factor(set)1               0.629649941 0.4969798210 7.971138e-01
## as.factor(sex)1:Patient_age   1.013898787 0.9932203181 1.035227e+00
## Patient_age:as.factor(race)2  1.003461653 0.9840816996 1.023246e+00
## Patient_age:as.factor(race)3  1.012866116 0.9292344385 1.105846e+00
## Patient_age:as.factor(race)4  0.955031033 0.8607361275 1.043602e+00
## Patient_age:as.factor(race)5  1.077638306 1.0097984026 1.166720e+00
#Race removed - age is a numeric variable here and interacts with race
# ***** This is the best model to test ******
hin_mod_2 <- glm(as.factor(sent_bin) ~ as.factor(eth) + as.factor(sex)*Patient_age +   as.factor(ins.y) + as.factor(set), data = sent_04_23, family = "binomial")

summ(hin_mod_2)
## MODEL INFO:
## Observations: 1380 (437 missing obs. deleted)
## Dependent Variable: as.factor(sent_bin)
## Type: Generalized linear model
##   Family: binomial 
##   Link function: logit 
## 
## MODEL FIT:
## χ²(10) = 44.70, p = 0.00
## Pseudo-R² (Cragg-Uhler) = 0.04
## Pseudo-R² (McFadden) = 0.02
## AIC = 1829.83, BIC = 1887.36 
## 
## Standard errors: MLE
## ----------------------------------------------------------------
##                                      Est.   S.E.   z val.      p
## --------------------------------- ------- ------ -------- ------
## (Intercept)                          1.02   0.56     1.82   0.07
## as.factor(eth)1                     -0.31   0.20    -1.58   0.11
## as.factor(eth)2                     -1.63   0.62    -2.63   0.01
## as.factor(sex)1                     -1.02   0.58    -1.76   0.08
## Patient_age                         -0.03   0.01    -2.88   0.00
## as.factor(ins.y)2                    0.49   0.26     1.90   0.06
## as.factor(ins.y)3                    0.23   0.27     0.85   0.39
## as.factor(ins.y)4                   -0.19   0.42    -0.44   0.66
## as.factor(ins.y)5                    0.21   0.47     0.44   0.66
## as.factor(set)1                     -0.47   0.11    -4.12   0.00
## as.factor(sex)1:Patient_age          0.02   0.01     1.66   0.10
## ----------------------------------------------------------------

Race was removed in this model and age is categorical with no interactions

#with some mod 2
hin_mod_2a <- glm(as.factor(sent_bin) ~ as.factor(eth) + as.factor(sex) +   as.factor(ins.y) + as.factor(age_group) + as.factor(set), data = sent_04_23, family = "binomial")

summ(hin_mod_2a)
## MODEL INFO:
## Observations: 1240 (577 missing obs. deleted)
## Dependent Variable: as.factor(sent_bin)
## Type: Generalized linear model
##   Family: binomial 
##   Link function: logit 
## 
## MODEL FIT:
## χ²(11) = 38.16, p = 0.00
## Pseudo-R² (Cragg-Uhler) = 0.04
## Pseudo-R² (McFadden) = 0.02
## AIC = 1641.99, BIC = 1703.46 
## 
## Standard errors: MLE
## ----------------------------------------------------------
##                                Est.   S.E.   z val.      p
## --------------------------- ------- ------ -------- ------
## (Intercept)                   -0.78   0.39    -1.98   0.05
## as.factor(eth)1               -0.20   0.21    -0.94   0.35
## as.factor(eth)2               -2.16   0.74    -2.91   0.00
## as.factor(sex)1               -0.25   0.12    -2.01   0.04
## as.factor(ins.y)2              0.32   0.26     1.23   0.22
## as.factor(ins.y)3              0.22   0.27     0.81   0.42
## as.factor(ins.y)4             -0.35   0.45    -0.78   0.43
## as.factor(ins.y)5              0.18   0.49     0.36   0.72
## as.factor(age_group)2          0.63   0.30     2.10   0.04
## as.factor(age_group)3          0.22   0.31     0.72   0.47
## as.factor(age_group)4          0.45   0.34     1.33   0.18
## as.factor(set)1               -0.23   0.12    -1.87   0.06
## ----------------------------------------------------------
exp(cbind(OR = coef(hin_mod_2a), confint(hin_mod_2a)))
## Waiting for profiling to be done...
##                              OR      2.5 %    97.5 %
## (Intercept)           0.4606102 0.20977898 0.9798854
## as.factor(eth)1       0.8207528 0.54036973 1.2339307
## as.factor(eth)2       0.1147512 0.01828114 0.3932036
## as.factor(sex)1       0.7824357 0.61588635 0.9941461
## as.factor(ins.y)2     1.3838123 0.83065365 2.3520401
## as.factor(ins.y)3     1.2471563 0.73511813 2.1535251
## as.factor(ins.y)4     0.7043600 0.28359006 1.6649588
## as.factor(ins.y)5     1.1935748 0.44068997 3.1035839
## as.factor(age_group)2 1.8802396 1.05932886 3.4633754
## as.factor(age_group)3 1.2464602 0.69350850 2.3199102
## as.factor(age_group)4 1.5642872 0.81995488 3.0710201
## as.factor(set)1       0.7949885 0.62501836 1.0110070

Testing Naive Bayes classification

Creating train and test data sets for training the model and testing - Naive Bayes Classifier with a 80% 20% test train

set.seed(1234)
ind <- sample(2, nrow(sent_04_23), replace = T, prob = c(0.8, 0.2))
train <- sent_04_23[ind == 1,]
test <- sent_04_23[ind == 2,]

NB classification in R

Prediction from the Naive Bayes classifier

Training data

p <- predict(model_train, train, type = 'prob')
## Warning: predict.naive_bayes(): more features in the newdata are provided as
## there are probability tables in the object. Calculation is performed based on
## features to be found in the tables.
head(cbind(p, train))
##           0            1    X    sid sentiment encounter_id.x note_year
## 1 0.9996237 3.762773e-04    1 984695  negative    10023548269      2016
## 2 0.9999448 5.520829e-05  160 164538  negative    10023586843      2016
## 3 0.9999453 5.472220e-05  388 164538  negative    10023587511      2016
## 4 0.9999941 5.944054e-06  616 202237  negative    10023570065      2016
## 6 0.9999764 2.358839e-05 2212 241264  negative    10023554164      2016
## 7 0.9998754 1.246258e-04 2461 905558  negative    10023552295      2016
##   agegroup race sex sent_bin eth age_group Patient_age charlson_score ins.y
## 1    74-81    1   0        0   0         4    77.97399             NA     2
## 2    66-73    2   1        0   0         3    66.97331             NA     2
## 3    66-73    2   1        0   0         3    66.97331             NA     2
## 4    58-65    2   0        0   0         3    58.42574             NA     2
## 6    66-73    2   1        0   0         3    69.01300             NA     2
## 7    58-65    2   0        0   0         3    65.27036             NA     2
##   insurance.y         enc_type                        DEPARTMENT_NAME
## 1    Medicare     Office Visit           EH SC PRIMARY CARE 6TH FLOOR
## 2    Medicare Clinical Support                EH HEALTH WOUND THERAPY
## 3    Medicare Clinical Support                EH HEALTH WOUND THERAPY
## 4    Medicare     Office Visit       EH SM PULMONARY CLINIC 4TH FLOOR
## 6    Medicare     Office Visit ZZZ *INACTIVE* EH W 38TH ADULT PRIMARY
## 7    Medicare     Office Visit               EH SSC UROLOGY 3RD FLOOR
##                SPECIALTY outpx_year outpx_month set
## 1 Geriatric Primary Care       2016          10   1
## 2         Rehabilitation       2016          10   1
## 3         Rehabilitation       2016          10   1
## 4            Pulmonology       2016          10   1
## 6           Primary Care       2016          10   1
## 7                Urology       2016          10   0

The training model is 99.9% accurate with a mis-classification rate of 0.0006821282

p <- predict(model_train, train) 
## Warning: predict.naive_bayes(): more features in the newdata are provided as
## there are probability tables in the object. Calculation is performed based on
## features to be found in the tables.
tab1 <- table(p, train$sent_bin)

1 - sum(diag(tab1)) / sum(tab1)
## [1] 0.0006821282

The test model is 98.6% accurate with a mis-classification rate of 0.01424501

p_test <- predict(model_test, test) 
## Warning: predict.naive_bayes(): more features in the newdata are provided as
## there are probability tables in the object. Calculation is performed based on
## features to be found in the tables.
tab1_test <- table(p_test, test$sent_bin)

1 - sum(diag(tab1_test)) / sum(tab1_test)
## [1] 0.01424501

Testing Random Forest approach in R - 70%, 30% test train

set.seed(222)
ind_RF <- sample(2, nrow(sent_04_23), replace = TRUE, prob = c(0.7, 0.3))
train_RF <- sent_04_23[ind==1,]
test_RF <- sent_04_23[ind==2,]

The error rate is 0%, which seems too good.

rf_train <- randomForest(sent_bin ~., data=train_RF, proximity=TRUE, na.action = na.exclude) 
## Warning in randomForest.default(m, y, ...): The response has five or fewer
## unique values. Are you sure you want to do regression?
print(rf_train)
## 
## Call:
##  randomForest(formula = sent_bin ~ ., data = train_RF, proximity = TRUE,      na.action = na.exclude) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 6
## 
##           Mean of squared residuals: 0.01814366
##                     % Var explained: 92.05

Prediction analytics and Confusion Matrix for train data

Since these data are a bit messy, especially given the frequency drops substantially. The prediction may be off? 100% accuracy?

#p1_RF <- predict(rf_train, train_RF)
#confusionMatrix(p1_RF, train_RF$sent_bin)

Prediction analytics and Confusion Matrix for test data

Since these data are a bit messy, especially given the frequency drops substantially. The prediction may be off? 100% accuracy?

#p2_RF <- predict(rf, test_RF)
#confusionMatrix(p2_RF, test_RF$sent_bin)

Error rate plot from Random Forest Model

plot(rf_train)

Building a neural network using Sentiment Data

We use 12 neurons in the first layer and 7 neurons in the second layer.

sent_04_23 %<>% 
  mutate_if(is.factor, as.numeric)

nn_mm <- model.matrix(~sex+eth+Patient_age+ins.y+set+sent_bin,
               data = sent_04_23)

Creating a formula to build the neural network

col_list <- paste(c(colnames(nn_mm[,-c(1,1380)])),collapse="+")
col_list <- paste(c("sent_bin~",col_list),collapse="")
f <- formula(col_list)

Default Neural Network With One Hidden Node (RPROP+)

With a single hidden layer, ANN predicts sentiment with 99.99% accuracy.

set.seed(8745212)

nmod <- neuralnet(f,data=nn_mm,hidden=1,
                    threshold = 0.01,
                    learningrate.limit = NULL,
                    learningrate.factor =
                         list(minus = 0.5, plus = 1.2),
                    algorithm = "rprop+")

nmod$result.matrix
##                                  [,1]
## error                    1.195490e-06
## reached.threshold        5.168014e-03
## steps                    6.280000e+02
## Intercept.to.1layhid1   -1.475786e+00
## sex.to.1layhid1         -2.078494e-05
## eth.to.1layhid1         -1.398440e-05
## Patient_age.to.1layhid1 -1.559939e-05
## ins.y.to.1layhid1       -3.013378e-04
## set.to.1layhid1         -1.637795e-04
## sent_bin.to.1layhid1     5.062211e+00
## Intercept.to.sent_bin   -2.360234e-01
## 1layhid1.to.sent_bin     1.270315e+00

Changing Resilient Backpropogation algorithm (RPROP+)

The Artificial Neural Network (ANN) with 5 hidden layers predicted sentiment with 99.78% accuracy.

set.seed(7896129)
nn5 <- neuralnet(f,data=nn_mm,hidden=5)

nn5$result.matrix
##                                  [,1]
## error                    3.150236e-04
## reached.threshold        9.558763e-03
## steps                    1.130000e+03
## Intercept.to.1layhid1   -5.574509e+00
## sex.to.1layhid1         -3.828967e+01
## eth.to.1layhid1         -8.844675e-01
## Patient_age.to.1layhid1 -5.378817e-02
## ins.y.to.1layhid1        2.708461e-01
## set.to.1layhid1          1.737011e+00
## sent_bin.to.1layhid1     8.577215e+01
## Intercept.to.1layhid2    5.072177e-01
## sex.to.1layhid2         -3.476411e+00
## eth.to.1layhid2         -1.866817e+00
## Patient_age.to.1layhid2 -1.521402e-01
## ins.y.to.1layhid2        7.873123e-01
## set.to.1layhid2         -4.261276e+01
## sent_bin.to.1layhid2    -7.182104e+00
## Intercept.to.1layhid3   -7.421607e-01
## sex.to.1layhid3         -1.004347e-01
## eth.to.1layhid3         -4.285776e-01
## Patient_age.to.1layhid3 -3.776887e-02
## ins.y.to.1layhid3        7.836825e-02
## set.to.1layhid3          2.882051e-01
## sent_bin.to.1layhid3     1.304398e+01
## Intercept.to.1layhid4    1.017103e+00
## sex.to.1layhid4         -1.784830e-01
## eth.to.1layhid4         -6.501898e-01
## Patient_age.to.1layhid4 -5.633848e-02
## ins.y.to.1layhid4        1.171013e-01
## set.to.1layhid4          4.317848e-01
## sent_bin.to.1layhid4    -7.917356e+00
## Intercept.to.1layhid5   -3.903659e-01
## sex.to.1layhid5          3.091168e+00
## eth.to.1layhid5         -5.695218e-02
## Patient_age.to.1layhid5 -1.686752e-02
## ins.y.to.1layhid5       -1.007441e+01
## set.to.1layhid5         -8.241955e+01
## sent_bin.to.1layhid5    -4.300281e+00
## Intercept.to.sent_bin   -9.182441e-03
## 1layhid1.to.sent_bin     4.316397e-01
## 1layhid2.to.sent_bin     3.409570e-02
## 1layhid3.to.sent_bin     5.775822e-01
## 1layhid4.to.sent_bin    -2.085470e-01
## 1layhid5.to.sent_bin    -2.257257e+00