1102hw1

input data

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

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)

with(dta, table(ill, gender))
with(dta, table(ill, university))
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)

plot by gender ~university

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))

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"))

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

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