library(readstata13)

tai = read.dta13("D:\\OneDrive - UMP\\R - lenh 2016\\hc.dta")

## Warning in read.dta13("D:\\OneDrive - UMP\\R - lenh 2016\\hc.dta"): 
##    Factor codes of type double or float detected in variables
## 
##    nhomtuoi, xlloiich, xlhanche, xltacdong,
##    xlnhucau
## 
##    No labels have been assigned.
##    Set option 'nonint.factors = TRUE' to assign labels anyway.

attach(tai)
names(tai)

##  [1] "stt"         "namsinh"     "nganh"       "gioitinh"    "quequan"    
##  [6] "songchung"   "hoctap"      "chucvu"      "kinhte"      "thamgiaclb" 
## [11] "quantamcn"   "khoahocnnlt" "hoinghikh"   "hieubiet"    "nguon"      
## [16] "khoahockhmt" "baigiang"    "nnlt"        "ungdung"     "chuyenkhoa" 
## [21] "truonghoc"   "benhvien"    "ungdungbv"   "l1"          "l2"         
## [26] "l3"          "l4"          "l5"          "l6"          "l7"         
## [31] "l8"          "l9"          "l10"         "l11"         "l12"        
## [36] "l13"         "l14"         "l15"         "l16"         "l17"        
## [41] "h1"          "h2"          "h3"          "h4"          "h5"         
## [46] "h6"          "t1"          "t2"          "t3"          "t4"         
## [51] "t5"          "t6"          "t7"          "t8"          "t9"         
## [56] "t10"         "n1"          "n2"          "n3"          "n4"         
## [61] "n5"          "n6"          "n7"          "n8"          "n9"         
## [66] "n10"         "n11"         "n12"         "n13"         "tuoi"       
## [71] "nhomtuoi"    "dtbloiich"   "dtbhanche"   "dtbtacdong"  "dtbnhucau"  
## [76] "loiich100"   "hanche100"   "tacdong100"  "nhucau100"   "xlloiich"   
## [81] "xlhanche"    "xltacdong"   "xlnhucau"

# hoi qui logistic da thuc MRL 
# outcom co 3 gia tri 

library(nnet)
library(table1)

## 
## Attaching package: 'table1'

## The following objects are masked from 'package:base':
## 
##     units, units<-

table1(~ xlnhucau + nganh)

	Overall (N=362)
xlnhucau
Mean (SD)	0.823 (0.455)
Median [Min, Max]	1.00 [0, 2.00]
nganh
ktha	60 (16.6%)
y da khoa	248 (68.5%)
rang ham mat	54 (14.9%)

library(forcats)

tai$nhucau = factor(tai$xlnhucau, levels = c(0,1,2), labels=c("kem", "trungb", "tot"))

tai$xlnhucau.1 = relevel(tai$nhucau, ref="kem")

m = multinom(xlnhucau.1 ~ factor(gioitinh), data=tai)

## # weights:  9 (4 variable)
## initial  value 397.697648 
## iter  10 value 224.466565
## final  value 224.466560 
## converged

summary(m)

## Call:
## multinom(formula = xlnhucau.1 ~ factor(gioitinh), data = tai)
## 
## Coefficients:
##        (Intercept) factor(gioitinh)nam
## trungb    1.011810           0.8476004
## tot      -2.926037           2.0608997
## 
## Std. Errors:
##        (Intercept) factor(gioitinh)nam
## trungb   0.1560552           0.2918457
## tot      0.5924634           0.7270734
## 
## Residual Deviance: 448.9331 
## AIC: 456.9331

exp(coef(m))

##        (Intercept) factor(gioitinh)nam
## trungb  2.75057604            2.334039
## tot     0.05360909            7.853032

zvalues = summary(m)$coefficients/summary(m)$standard.errors
pnorm(abs(zvalues), lower.tail = F)*2

##         (Intercept) factor(gioitinh)nam
## trungb 8.951756e-11         0.003681041
## tot    7.861960e-07         0.004589544

predicted = fitted(m)
head(predicted)

##         kem    trungb        tot
## 1 0.2628684 0.7230395 0.01409214
## 2 0.1275356 0.8187727 0.05369172
## 3 0.2628684 0.7230395 0.01409214
## 4 0.2628684 0.7230395 0.01409214
## 5 0.1275356 0.8187727 0.05369172
## 6 0.1275356 0.8187727 0.05369172

tai$xlnhucau.1 = relevel(tai$nhucau, ref="kem")
m.1 = multinom(xlnhucau.1 ~ nganh, data = tai) 

## # weights:  12 (6 variable)
## initial  value 397.697648 
## iter  10 value 228.215276
## iter  20 value 228.064743
## final  value 228.054455 
## converged

summary(m.1)

## Call:
## multinom(formula = xlnhucau.1 ~ nganh, data = tai)
## 
## Coefficients:
##        (Intercept) nganhy da khoa nganhrang ham mat
## trungb    1.493400     -0.3038171         0.1886166
## tot      -9.613468      7.6674797         8.6334162
## 
## Std. Errors:
##        (Intercept) nganhy da khoa nganhrang ham mat
## trungb    0.333589      0.3668425         0.5094784
## tot      36.877861     36.8797985        36.8840732
## 
## Residual Deviance: 456.1089 
## AIC: 468.1089

exp(coef(m.1))

##         (Intercept) nganhy da khoa nganhrang ham mat
## trungb 4.452206e+00      0.7379958          1.207578
## tot    6.682269e-05   2137.6869740       5616.231734

zvalues = summary(m.1)$coefficients/summary(m.1)$standard.errors
pnorm(abs(zvalues), lower.tail = F)*2

##         (Intercept) nganhy da khoa nganhrang ham mat
## trungb 7.578217e-06      0.4075601         0.7112223
## tot    7.943362e-01      0.8353035         0.8149315

predicted.1 = fitted(m.1)
head(predicted.1)

##         kem    trungb        tot
## 1 0.2258073 0.7419371 0.03225566
## 2 0.2258073 0.7419371 0.03225566
## 3 0.2258073 0.7419371 0.03225566
## 4 0.2258073 0.7419371 0.03225566
## 5 0.2258073 0.7419371 0.03225566
## 6 0.2258073 0.7419371 0.03225566