1 input data

data(addiction, package = "catdata")
dta<-addiction
head(dta)

##   ill gender age university
## 1   1      1  61          0
## 2   0      1  43          0
## 3   2      0  44          0
## 4   0      1  21          1
## 5   0      0  33          0
## 6   1      0  83          0

str(dta)

## 'data.frame':    712 obs. of  4 variables:
##  $ ill       : int  1 0 2 0 0 1 0 1 1 1 ...
##  $ gender    : int  1 1 0 1 0 0 0 0 1 0 ...
##  $ age       : int  61 43 44 21 33 83 29 61 37 19 ...
##  $ university: int  0 0 0 1 0 0 0 0 0 0 ...

dta$gender <- factor(dta$gender, 0:1, labels = c("M", "F"))
dta$ill <- factor(dta$ill, 0:2, labels = c("weak-willed","disease","both"))
dta$university <- factor(dta$university, 0:1, labels = c("No","Yes"))
str(dta)

## 'data.frame':    712 obs. of  4 variables:
##  $ ill       : Factor w/ 3 levels "weak-willed",..: 2 1 3 1 1 2 1 2 2 2 ...
##  $ gender    : Factor w/ 2 levels "M","F": 2 2 1 2 1 1 1 1 2 1 ...
##  $ age       : int  61 43 44 21 33 83 29 61 37 19 ...
##  $ university: Factor w/ 2 levels "No","Yes": 1 1 1 2 1 1 1 1 1 1 ...

with(dta, table(ill, gender))

##              gender
## ill             M   F
##   weak-willed  99 103
##   disease     137 179
##   both         79  90

with(dta, table(ill, university))

##              university
## ill            No Yes
##   weak-willed 178  23
##   disease     194 121
##   both        123  44

HH::likert(t(with(dta, table(ill, gender))), 
           main="", 
           ylab="gender", 
           as.percent=TRUE)

HH::likert(t(with(dta, table(ill, university))), 
           main="", 
           ylab="university", 
           as.percent=TRUE)

2 plot by gender ~university

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)
dta3 <- dta %>% mutate(resp_n=ifelse(ill=='weak-willed', 1, 0),
                       resp_f=ifelse(ill=='disease', 1, 0),
                       resp_p=ifelse(ill=='both', 1, 0))
ggplot(dta3, aes(age, resp_n)) + 
 stat_smooth(method="glm", formula=y ~ x, 
             method.args=list(family=binomial),
             alpha=.2, se=FALSE, col=3) +
 stat_smooth(aes(y=resp_f),
             method="glm", formula=y ~ x, 
             method.args=list(family=binomial),
             alpha=.2, se=FALSE, col=2) +
 stat_smooth(aes(y=resp_p), method="glm", 
             formula=y ~ x, 
             method.args=list(family=binomial),
             alpha=.2, se=FALSE, col=4) +
 facet_grid( gender ~university )+
 scale_y_continuous(limits=c(0,1), breaks=seq(0, 1, by=.1))+
 guides(color=guide_legend(reverse=T))+
 labs(x="age", 
      y="Probability of ill",
      title="weak-willed=Green, disease=Red, both=Blue") +
 theme_minimal()+
 theme(legend.position = c(.9, .2))

## Warning: Removed 26 rows containing non-finite values (stat_smooth).

## Warning: Removed 26 rows containing non-finite values (stat_smooth).

## Warning: Removed 26 rows containing non-finite values (stat_smooth).

3 model

m0 <- nnet::multinom(ill ~ gender*age*university, data=dta, Hess=TRUE, trace=F)
m1 <- nnet::multinom(ill ~ gender+age+university, data=dta, Hess=TRUE, trace=F)
knitr::kable(anova(m0, m1, test="Chisq"))

Model	Resid. df	Resid. Dev	Test	Df	LR stat.	Pr(Chi)
gender + age + university	1356	1350.417		NA	NA	NA
gender * age * university	1348	1342.698	1 vs 2	8	7.718734	0.4614161

沒有顯著，用簡單的模型m1即可

sjPlot::tab_model(m1, show.r2=FALSE, show.obs=FALSE)

	ill
Predictors	Odds Ratios	CI	p	Response
(Intercept)	0.31	0.19 – 0.53	<0.001	disease
gender [F]	1.55	1.06 – 2.26	0.024	disease
age	1.03	1.02 – 1.04	<0.001	disease
university [Yes]	5.06	3.08 – 8.32	<0.001	disease
(Intercept)	0.13	0.07 – 0.24	<0.001	both
gender [F]	1.33	0.87 – 2.06	0.192	both
age	1.04	1.03 – 1.06	<0.001	both
university [Yes]	2.91	1.65 – 5.12	<0.001	both

1102hw1

1 input data

2 plot by gender ~university

3 model