library(readxl)
library(nnet)


# Read in the data
Analysis <- read_excel("~/Library/CloudStorage/OneDrive-Personal/1_Surface_pro_Narasimha/work/AIIMS_Deoghar/Projects/WHO_projects/WHO_COVID_Survey/Analysis/10_06_2023_Fresh_Analysis/Analysis.xlsx")
## New names:
## • `` -> `...3`
# Fit the multinomial logistic regression model
model <- multinom(Severity ~ Tobacco_use + Gender + Occupation + Background + Vaccination + Alcohol_use + Age_years + Co_Morbidity + SES, data = Analysis, maxit = 1000)
## # weights:  39 (24 variable)
## initial  value 1565.522511 
## iter  10 value 664.905660
## iter  20 value 652.835621
## iter  30 value 652.745493
## final  value 652.745096 
## converged
# summarize the results
summary(model)
## Call:
## multinom(formula = Severity ~ Tobacco_use + Gender + Occupation + 
##     Background + Vaccination + Alcohol_use + Age_years + Co_Morbidity + 
##     SES, data = Analysis, maxit = 1000)
## 
## Coefficients:
##          (Intercept) Tobacco_useYes  GenderMale OccupationWork from/at home
## Moderate   -2.177608      0.3384961 -0.09482542                   0.5686107
## Severe     -2.136700     -0.4339554  0.03677957                  -0.3460120
##          OccupationWork in/associated with hospital BackgroundUrban
## Moderate                                  0.2012718     -0.33570598
## Severe                                   -0.4967809     -0.05330268
##          VaccinationNot Completed Alcohol_useYes Age_yearsMiddle age
## Moderate                0.1791730    -0.08575303          -0.5357481
## Severe                  0.3887127     0.14794142          -1.4487024
##          Age_yearsYoung age Co_MorbidityYes SESBPL card holder
## Moderate         -0.5713309        1.082977         -0.5318476
## Severe           -2.0789061        1.207802         -1.2615770
## 
## Std. Errors:
##          (Intercept) Tobacco_useYes GenderMale OccupationWork from/at home
## Moderate   0.6177400      0.2291135  0.2292193                   0.2886111
## Severe     0.8773425      0.3257150  0.2813028                   0.4936754
##          OccupationWork in/associated with hospital BackgroundUrban
## Moderate                                  0.3744682       0.2756775
## Severe                                    0.6259435       0.3920901
##          VaccinationNot Completed Alcohol_useYes Age_yearsMiddle age
## Moderate                0.4989696      0.3220633            0.266678
## Severe                  0.7605862      0.4050153            0.276793
##          Age_yearsYoung age Co_MorbidityYes SESBPL card holder
## Moderate          0.3719060       0.2048159          0.2790671
## Severe            0.5732571       0.2645619          0.4563220
## 
## Residual Deviance: 1305.49 
## AIC: 1353.49
# Assume that your model object is called "model"

# Summarize the results
summary_model <- summary(model)
summary_model
## Call:
## multinom(formula = Severity ~ Tobacco_use + Gender + Occupation + 
##     Background + Vaccination + Alcohol_use + Age_years + Co_Morbidity + 
##     SES, data = Analysis, maxit = 1000)
## 
## Coefficients:
##          (Intercept) Tobacco_useYes  GenderMale OccupationWork from/at home
## Moderate   -2.177608      0.3384961 -0.09482542                   0.5686107
## Severe     -2.136700     -0.4339554  0.03677957                  -0.3460120
##          OccupationWork in/associated with hospital BackgroundUrban
## Moderate                                  0.2012718     -0.33570598
## Severe                                   -0.4967809     -0.05330268
##          VaccinationNot Completed Alcohol_useYes Age_yearsMiddle age
## Moderate                0.1791730    -0.08575303          -0.5357481
## Severe                  0.3887127     0.14794142          -1.4487024
##          Age_yearsYoung age Co_MorbidityYes SESBPL card holder
## Moderate         -0.5713309        1.082977         -0.5318476
## Severe           -2.0789061        1.207802         -1.2615770
## 
## Std. Errors:
##          (Intercept) Tobacco_useYes GenderMale OccupationWork from/at home
## Moderate   0.6177400      0.2291135  0.2292193                   0.2886111
## Severe     0.8773425      0.3257150  0.2813028                   0.4936754
##          OccupationWork in/associated with hospital BackgroundUrban
## Moderate                                  0.3744682       0.2756775
## Severe                                    0.6259435       0.3920901
##          VaccinationNot Completed Alcohol_useYes Age_yearsMiddle age
## Moderate                0.4989696      0.3220633            0.266678
## Severe                  0.7605862      0.4050153            0.276793
##          Age_yearsYoung age Co_MorbidityYes SESBPL card holder
## Moderate          0.3719060       0.2048159          0.2790671
## Severe            0.5732571       0.2645619          0.4563220
## 
## Residual Deviance: 1305.49 
## AIC: 1353.49
# Extract the estimated coefficients
coef <- coef(model)

# Exponentiate the coefficients to obtain the odds ratios
odds_ratios <- exp(coef)
odds_ratios_ci <- exp(confint(model, level = 0.95))

# Print the odds ratios
odds_ratios
##          (Intercept) Tobacco_useYes GenderMale OccupationWork from/at home
## Moderate   0.1133122      1.4028363  0.9095317                    1.765812
## Severe     0.1180438      0.6479412  1.0374643                    0.707504
##          OccupationWork in/associated with hospital BackgroundUrban
## Moderate                                  1.2229572       0.7148332
## Severe                                    0.6084863       0.9480930
##          VaccinationNot Completed Alcohol_useYes Age_yearsMiddle age
## Moderate                 1.196228      0.9178209           0.5852313
## Severe                   1.475081      1.1594450           0.2348749
##          Age_yearsYoung age Co_MorbidityYes SESBPL card holder
## Moderate          0.5647733        2.953458          0.5875185
## Severe            0.1250670        3.346121          0.2832070
odds_ratios_ci
## , , Moderate
## 
##                                                 2.5 %    97.5 %
## (Intercept)                                0.03376413 0.3802752
## Tobacco_useYes                             0.89533435 2.1980053
## GenderMale                                 0.58037137 1.4253769
## OccupationWork from/at home                1.00294774 3.1089281
## OccupationWork in/associated with hospital 0.58703406 2.5477640
## BackgroundUrban                            0.41643536 1.2270490
## VaccinationNot Completed                   0.44987321 3.1808090
## Alcohol_useYes                             0.48822176 1.7254355
## Age_yearsMiddle age                        0.34700112 0.9870162
## Age_yearsYoung age                         0.27246279 1.1706879
## Co_MorbidityYes                            1.97692941 4.4123540
## SESBPL card holder                         0.34000021 1.0152286
## 
## , , Severe
## 
##                                                 2.5 %    97.5 %
## (Intercept)                                0.02114695 0.6589286
## Tobacco_useYes                             0.34220508 1.2268308
## GenderMale                                 0.59776125 1.8006055
## OccupationWork from/at home                0.26885098 1.8618565
## OccupationWork in/associated with hospital 0.17842131 2.0751755
## BackgroundUrban                            0.43964608 2.0445544
## VaccinationNot Completed                   0.33220309 6.5497980
## Alcohol_useYes                             0.52420424 2.5644826
## Age_yearsMiddle age                        0.13653056 0.4040575
## Age_yearsYoung age                         0.04066167 0.3846802
## Co_MorbidityYes                            1.99226080 5.6200108
## SESBPL card holder                         0.11579293 0.6926695
z <- summary_model$coefficients/summary_model$standard.errors
p <- (1 - pnorm(abs(z), 0, 1))*2 # we are using two-tailed z test
p
##           (Intercept) Tobacco_useYes GenderMale OccupationWork from/at home
## Moderate 0.0004232898      0.1395639  0.6791022                  0.04881975
## Severe   0.0148744013      0.1827563  0.8959752                  0.48337193
##          OccupationWork in/associated with hospital BackgroundUrban
## Moderate                                  0.5909312       0.2233195
## Severe                                    0.4273984       0.8918648
##          VaccinationNot Completed Alcohol_useYes Age_yearsMiddle age
## Moderate                0.7195307      0.7900379        4.454030e-02
## Severe                  0.6093021      0.7149071        1.659851e-07
##          Age_yearsYoung age Co_MorbidityYes SESBPL card holder
## Moderate       0.1244834890    1.239576e-07        0.056675448
## Severe         0.0002873096    4.988032e-06        0.005698147
library(knitr)
Pclass2 <- rbind(summary_model$coefficients[1,],summary_model$standard.errors[1,],z[1,],p[1,])
rownames(Pclass2) <- c("Coefficient","Std. Errors","z stat","p value")
knitr::kable(Pclass2)
(Intercept) Tobacco_useYes GenderMale OccupationWork from/at home OccupationWork in/associated with hospital BackgroundUrban VaccinationNot Completed Alcohol_useYes Age_yearsMiddle age Age_yearsYoung age Co_MorbidityYes SESBPL card holder
Coefficient -2.1776081 0.3384961 -0.0948254 0.5686107 0.2012718 -0.3357060 0.1791730 -0.0857530 -0.5357481 -0.5713309 1.0829766 -0.5318476
Std. Errors 0.6177400 0.2291135 0.2292193 0.2886111 0.3744682 0.2756775 0.4989696 0.3220633 0.2666780 0.3719060 0.2048159 0.2790671
z stat -3.5251207 1.4774168 -0.4136886 1.9701625 0.5374871 -1.2177488 0.3590860 -0.2662614 -2.0089702 -1.5362238 5.2875616 -1.9058055
p value 0.0004233 0.1395639 0.6791022 0.0488198 0.5909312 0.2233195 0.7195307 0.7900379 0.0445403 0.1244835 0.0000001 0.0566754
Pclass3 <- rbind(summary_model$coefficients[2,],summary_model$standard.errors[2,],z[2,],p[2,])
rownames(Pclass3) <- c("Coefficient","Std. Errors","z stat","p value")
knitr::kable(Pclass3)
(Intercept) Tobacco_useYes GenderMale OccupationWork from/at home OccupationWork in/associated with hospital BackgroundUrban VaccinationNot Completed Alcohol_useYes Age_yearsMiddle age Age_yearsYoung age Co_MorbidityYes SESBPL card holder
Coefficient -2.1366999 -0.4339554 0.0367796 -0.3460120 -0.4967809 -0.0533027 0.3887127 0.1479414 -1.4487024 -2.0789061 1.2078018 -1.2615770
Std. Errors 0.8773425 0.3257150 0.2813028 0.4936754 0.6259435 0.3920901 0.7605862 0.4050153 0.2767930 0.5732571 0.2645619 0.4563220
z stat -2.4354227 -1.3323163 0.1307473 -0.7008896 -0.7936514 -0.1359450 0.5110699 0.3652737 -5.2338835 -3.6264813 4.5652905 -2.7646639
p value 0.0148744 0.1827563 0.8959752 0.4833719 0.4273984 0.8918648 0.6093021 0.7149071 0.0000002 0.0002873 0.0000050 0.0056981
# Fit the multinomial logistic regression model
model_complications <- multinom(Any_complication ~ Tobacco_use + Gender + Occupation + Background + Vaccination + Alcohol_use + Age_years + Co_Morbidity + SES, data = Analysis, maxit = 1000)
## # weights:  13 (12 variable)
## initial  value 987.734732 
## iter  10 value 773.421388
## iter  20 value 752.850775
## iter  20 value 752.850775
## iter  20 value 752.850775
## final  value 752.850775 
## converged
# summarize the results
summary(model_complications)
## Call:
## multinom(formula = Any_complication ~ Tobacco_use + Gender + 
##     Occupation + Background + Vaccination + Alcohol_use + Age_years + 
##     Co_Morbidity + SES, data = Analysis, maxit = 1000)
## 
## Coefficients:
##                                                  Values Std. Err.
## (Intercept)                                -2.762659089 0.4741675
## Tobacco_useYes                              0.458443043 0.1560563
## GenderMale                                 -0.004634258 0.1514233
## OccupationWork from/at home                 0.487846406 0.2163415
## OccupationWork in/associated with hospital  0.651642875 0.2394029
## BackgroundUrban                             0.204260205 0.1948345
## VaccinationNot Completed                    1.175754775 0.3937252
## Alcohol_useYes                              0.431592031 0.2007126
## Age_yearsMiddle age                        -0.036655639 0.1952475
## Age_yearsYoung age                         -0.384044784 0.2606596
## Co_MorbidityYes                             1.286184811 0.1315558
## SESBPL card holder                         -0.703859410 0.1842346
## 
## Residual Deviance: 1505.702 
## AIC: 1529.702
# Assume that your model object is called "model"

# Summarize the results
summary_model_complications <- summary(model_complications)
summary_model_complications
## Call:
## multinom(formula = Any_complication ~ Tobacco_use + Gender + 
##     Occupation + Background + Vaccination + Alcohol_use + Age_years + 
##     Co_Morbidity + SES, data = Analysis, maxit = 1000)
## 
## Coefficients:
##                                                  Values Std. Err.
## (Intercept)                                -2.762659089 0.4741675
## Tobacco_useYes                              0.458443043 0.1560563
## GenderMale                                 -0.004634258 0.1514233
## OccupationWork from/at home                 0.487846406 0.2163415
## OccupationWork in/associated with hospital  0.651642875 0.2394029
## BackgroundUrban                             0.204260205 0.1948345
## VaccinationNot Completed                    1.175754775 0.3937252
## Alcohol_useYes                              0.431592031 0.2007126
## Age_yearsMiddle age                        -0.036655639 0.1952475
## Age_yearsYoung age                         -0.384044784 0.2606596
## Co_MorbidityYes                             1.286184811 0.1315558
## SESBPL card holder                         -0.703859410 0.1842346
## 
## Residual Deviance: 1505.702 
## AIC: 1529.702
# Extract the estimated coefficients
coef_complications <- coef(model_complications)

# Exponentiate the coefficients to obtain the odds ratios
odds_ratios_complications <- exp(coef_complications)
odds_ratios_ci_complications <- exp(confint(model_complications, level = 0.95))

# Print the odds ratios
odds_ratios_complications
##                                (Intercept) 
##                                 0.06312369 
##                             Tobacco_useYes 
##                                 1.58160957 
##                                 GenderMale 
##                                 0.99537646 
##                OccupationWork from/at home 
##                                 1.62880466 
## OccupationWork in/associated with hospital 
##                                 1.91869041 
##                            BackgroundUrban 
##                                 1.22661728 
##                   VaccinationNot Completed 
##                                 3.24058793 
##                             Alcohol_useYes 
##                                 1.53970683 
##                        Age_yearsMiddle age 
##                                 0.96400805 
##                         Age_yearsYoung age 
##                                 0.68110092 
##                            Co_MorbidityYes 
##                                 3.61895319 
##                         SESBPL card holder 
##                                 0.49467247
odds_ratios_ci_complications
##                                                 2.5 %    97.5 %
## (Intercept)                                0.02492185 0.1598838
## Tobacco_useYes                             1.16483379 2.1475071
## GenderMale                                 0.73976808 1.3393039
## OccupationWork from/at home                1.06590568 2.4889675
## OccupationWork in/associated with hospital 1.20012035 3.0675031
## BackgroundUrban                            0.83727025 1.7970183
## VaccinationNot Completed                   1.49790516 7.0107310
## Alcohol_useYes                             1.03894174 2.2818384
## Age_yearsMiddle age                        0.65748475 1.4134343
## Age_yearsYoung age                         0.40863691 1.1352339
## Co_MorbidityYes                            2.79642069 4.6834234
## SESBPL card holder                         0.34474420 0.7098041
z_complications <- summary_model_complications$coefficients/summary_model_complications$standard.errors
p_complications <- (1 - pnorm(abs(z_complications), 0, 1))*2 # we are using two-tailed z test
p_complications
##                                (Intercept) 
##                               5.665768e-09 
##                             Tobacco_useYes 
##                               3.306802e-03 
##                                 GenderMale 
##                               9.755848e-01 
##                OccupationWork from/at home 
##                               2.413440e-02 
## OccupationWork in/associated with hospital 
##                               6.489780e-03 
##                            BackgroundUrban 
##                               2.944645e-01 
##                   VaccinationNot Completed 
##                               2.824383e-03 
##                             Alcohol_useYes 
##                               3.153157e-02 
##                        Age_yearsMiddle age 
##                               8.510810e-01 
##                         Age_yearsYoung age 
##                               1.406547e-01 
##                            Co_MorbidityYes 
##                               0.000000e+00 
##                         SESBPL card holder 
##                               1.332075e-04
library(knitr)
Pclass2_complications <- rbind(summary_model_complications$coefficients,summary_model_complications$standard.errors,z_complications,p_complications)
rownames(Pclass2_complications) <- c("Coefficient","Std. Errors","z stat","p value")
knitr::kable(Pclass2_complications)
(Intercept) Tobacco_useYes GenderMale OccupationWork from/at home OccupationWork in/associated with hospital BackgroundUrban VaccinationNot Completed Alcohol_useYes Age_yearsMiddle age Age_yearsYoung age Co_MorbidityYes SESBPL card holder
Coefficient -2.7626591 0.4584430 -0.0046343 0.4878464 0.6516429 0.2042602 1.1757548 0.4315920 -0.0366556 -0.3840448 1.2861848 -0.7038594
Std. Errors 0.4741675 0.1560563 0.1514233 0.2163415 0.2394029 0.1948345 0.3937252 0.2007126 0.1952475 0.2606596 0.1315558 0.1842346
z stat -5.8263354 2.9376781 -0.0306046 2.2549829 2.7219508 1.0483781 2.9862320 2.1502991 -0.1877394 -1.4733574 9.7767280 -3.8204517
p value 0.0000000 0.0033068 0.9755848 0.0241344 0.0064898 0.2944645 0.0028244 0.0315316 0.8510810 0.1406547 0.0000000 0.0001332