Inquiries regarding this code contact gabriel.carrasco@upch.pe

rm(list=ls())
library(tidyverse); library(epiDisplay); library(tidyr); library(nnet); library(questionr)
library(AER); library(afex); library(car); library(lmtest); library(haven)

Example dep: 3cat, indep: 2 cat

Go to SAS output, pp 115-119

rm(list=ls())
data <- data.frame(outcome=c(0,0,1,1,2,2),
                   freckling=c(0,1,0,1,0,1),
                   weights=c(402,114,58,18,155,129))
data$outcome <- factor(data$outcome,
                    levels = c(0,1,2),
                    labels = c("controls", "skin_les", "melanoma"))
data$freckling <- factor(data$freckling,
                       levels = c(0,1),
                       labels = c("No", "Yes"))

wtd.table(data$outcome, data$freckling, weights = data$weights) # pp7 Notes
##           No Yes
## controls 402 114
## skin_les  58  18
## melanoma 155 129
wtd.table(data$outcome, weights = data$weights)
## controls skin_les melanoma 
##      516       76      284
test <- multinom(outcome ~ freckling, weights = weights, data = data)
## # weights:  9 (4 variable)
## initial  value 962.384365 
## iter  10 value 754.958109
## final  value 754.958093 
## converged
summary(test) # Residual Deviance = -2LogL in SAS
## Call:
## multinom(formula = outcome ~ freckling, data = data, weights = weights)
## 
## Coefficients:
##          (Intercept) frecklingYes
## skin_les  -1.9356147   0.08876193
## melanoma  -0.9531096   1.07641979
## 
## Std. Errors:
##          (Intercept) frecklingYes
## skin_les   0.1404367    0.2899875
## melanoma   0.0945508    0.1595695
## 
## Residual Deviance: 1509.916 
## AIC: 1517.916
coeftest(test)
## 
## z test of coefficients:
## 
##                        Estimate Std. Error  z value  Pr(>|z|)    
## skin_les:(Intercept)  -1.935615   0.140437 -13.7828 < 2.2e-16 ***
## skin_les:frecklingYes  0.088762   0.289988   0.3061    0.7595    
## melanoma:(Intercept)  -0.953110   0.094551 -10.0804 < 2.2e-16 ***
## melanoma:frecklingYes  1.076420   0.159570   6.7458 1.522e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef(test)) # OR
##          (Intercept) frecklingYes
## skin_les   0.1443355     1.092820
## melanoma   0.3855403     2.934156
exp(confint(test)) # 95% CI - OR
## , , skin_les
## 
##                  2.5 %    97.5 %
## (Intercept)  0.1096057 0.1900698
## frecklingYes 0.6190290 1.9292417
## 
## , , melanoma
## 
##                  2.5 %    97.5 %
## (Intercept)  0.3203229 0.4640359
## frecklingYes 2.1461366 4.0115202
set_sum_contrasts() # use sum coding, necessary to make type III LR tests valid
Anova(test,type="III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: outcome
##           LR Chisq Df Pr(>Chisq)    
## freckling   47.654  2  4.488e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Example dep: 3cat, indep: 4 cat

Go to SAS output, pp 120-125

rm(list=ls())
data <- data.frame(outcome=c(0,0,0,0,1,1,1,1,2,2,2,2),
                   tanning=c(0,1,2,3,0,1,2,3,0,1,2,3),
                   weights=c(42,105,222,154,11,28,22,14,34,90,101,48))
data$outcome <- factor(data$outcome,
                       levels = c(0,1,2),
                       labels = c("controls", "skin_les", "melanoma"))
data$tanning <- factor(data$tanning,
                       levels = c(0,1,2,3),
                       labels = c("None", "Light", "Average","Dark"))

wtd.table(data$outcome, data$tanning, weights = data$weights) # pp11 Notes
##          None Light Average Dark
## controls   42   105     222  154
## skin_les   11    28      22   14
## melanoma   34    90     101   48
wtd.table(data$outcome, weights = data$weights)
## controls skin_les melanoma 
##      523       75      273
test <- multinom(outcome ~ tanning, weights = weights, data = data)
## # weights:  15 (8 variable)
## initial  value 956.891303 
## iter  10 value 750.774909
## final  value 749.230971 
## converged
summary(test) # Residual Deviance = -2LogL in SAS
## Call:
## multinom(formula = outcome ~ tanning, data = data, weights = weights)
## 
## Coefficients:
##          (Intercept)  tanning1  tanning2   tanning3
## skin_les  -1.8427601 0.5029940 0.5210152 -0.4688781
## melanoma  -0.5796912 0.3683838 0.4255423 -0.2078637
## 
## Std. Errors:
##          (Intercept)  tanning1  tanning2  tanning3
## skin_les  0.13412607 0.2744974 0.2015160 0.2072906
## melanoma  0.08499582 0.1839428 0.1324453 0.1201126
## 
## Residual Deviance: 1498.462 
## AIC: 1514.462
coeftest(test)
## 
## z test of coefficients:
## 
##                       Estimate Std. Error  z value  Pr(>|z|)    
## skin_les:(Intercept) -1.842760   0.134126 -13.7390 < 2.2e-16 ***
## skin_les:tanning1     0.502994   0.274497   1.8324  0.066889 .  
## skin_les:tanning2     0.521015   0.201516   2.5855  0.009724 ** 
## skin_les:tanning3    -0.468878   0.207291  -2.2619  0.023701 *  
## melanoma:(Intercept) -0.579691   0.084996  -6.8202 9.089e-12 ***
## melanoma:tanning1     0.368384   0.183943   2.0027  0.045209 *  
## melanoma:tanning2     0.425542   0.132445   3.2130  0.001314 ** 
## melanoma:tanning3    -0.207864   0.120113  -1.7306  0.083528 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef(test)) # OR
##          (Intercept) tanning1 tanning2  tanning3
## skin_les   0.1583797 1.653665 1.683736 0.6257038
## melanoma   0.5600713 1.445397 1.530420 0.8123177
exp(confint(test)) # 95% CI - OR
## , , skin_les
## 
##                 2.5 %    97.5 %
## (Intercept) 0.1217674 0.2060003
## tanning1    0.9655947 2.8320449
## tanning2    1.1343401 2.4992217
## tanning3    0.4167952 0.9393229
## 
## , , melanoma
## 
##                 2.5 %    97.5 %
## (Intercept) 0.4741270 0.6615945
## tanning1    1.0078935 2.0728098
## tanning2    1.1805190 1.9840303
## tanning3    0.6419274 1.0279357
set_sum_contrasts() # use sum coding, necessary to make type III LR tests valid
Anova(test,type="III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: outcome
##         LR Chisq Df Pr(>Chisq)    
## tanning   36.337  6   2.37e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Example dep: 3cat, indep: continuous

Go to SAS output, pp 120-125

head(dat)
## # A tibble: 6 x 4
##   SHOCKTYP    Prin1  Prin2    UO
##   <fct>       <dbl>  <dbl> <dbl>
## 1 NONSHOCK  1.17     0.527   110
## 2 NONSHOCK  0.658    3.16     40
## 3 NONSHOCK -0.00619  2.41      0
## 4 NONSHOCK  2.51     0.227     0
## 5 NONSHOCK  1.84     0.710    42
## 6 NONSHOCK -1.50    -0.537     0
tab1(dat$SHOCKTYP)

## dat$SHOCKTYP : 
##          Frequency Percent Cum. percent
## NONSHOCK        34    47.9         47.9
## HYPOVOL         17    23.9         71.8
## CARDIO          20    28.2        100.0
##   Total         71   100.0        100.0
test <- multinom(SHOCKTYP ~ ., data = dat)
## # weights:  15 (8 variable)
## initial  value 78.001472 
## iter  10 value 53.853016
## final  value 53.189808 
## converged
summary(test) # Residual Deviance = -2LogL in SAS
## Call:
## multinom(formula = SHOCKTYP ~ ., data = dat)
## 
## Coefficients:
##         (Intercept)       Prin1      Prin2           UO
## HYPOVOL  -0.4426198 -0.09010694 -0.6198120 -0.002351107
## CARDIO   -0.7608403  0.84711096 -0.7192353 -0.013319075
## 
## Std. Errors:
##         (Intercept)     Prin1     Prin2         UO
## HYPOVOL    0.380034 0.1910667 0.2423933 0.00348166
## CARDIO     0.485351 0.2521337 0.2682082 0.01025010
## 
## Residual Deviance: 106.3796 
## AIC: 122.3796
coeftest(test) # p-values
## 
## z test of coefficients:
## 
##                       Estimate Std. Error z value  Pr(>|z|)    
## HYPOVOL:(Intercept) -0.4426198  0.3800340 -1.1647 0.2441466    
## HYPOVOL:Prin1       -0.0901069  0.1910667 -0.4716 0.6372127    
## HYPOVOL:Prin2       -0.6198120  0.2423933 -2.5571 0.0105564 *  
## HYPOVOL:UO          -0.0023511  0.0034817 -0.6753 0.4994957    
## CARDIO:(Intercept)  -0.7608403  0.4853510 -1.5676 0.1169725    
## CARDIO:Prin1         0.8471110  0.2521337  3.3598 0.0007801 ***
## CARDIO:Prin2        -0.7192353  0.2682082 -2.6816 0.0073264 ** 
## CARDIO:UO           -0.0133191  0.0102501 -1.2994 0.1938035    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef(test)) # OR
##         (Intercept)     Prin1     Prin2        UO
## HYPOVOL   0.6423514 0.9138335 0.5380456 0.9976517
## CARDIO    0.4672736 2.3328973 0.4871246 0.9867692
exp(confint(test)) # 95% CI - OR
## , , HYPOVOL
## 
##                 2.5 %   97.5 %
## (Intercept) 0.3049911 1.352876
## Prin1       0.6283922 1.328934
## Prin2       0.3345750 0.865256
## UO          0.9908669 1.004483
## 
## , , CARDIO
## 
##                 2.5 %    97.5 %
## (Intercept) 0.1804844 1.2097699
## Prin1       1.4232429 3.8239499
## Prin2       0.2879658 0.8240229
## UO          0.9671430 1.0067937
set_sum_contrasts() # use sum coding, necessary to make type III LR tests valid
Anova(test,type="III") # p-value for general variable across categories in the dependent variable
## Analysis of Deviance Table (Type III tests)
## 
## Response: SHOCKTYP
##       LR Chisq Df Pr(>Chisq)    
## Prin1  21.5865  2  2.054e-05 ***
## Prin2  13.1803  2   0.001374 ** 
## UO      2.8598  2   0.239331    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1