# HW 4
# 6.18
library(MASS)
library(VGAM)
Alligators2 <- read.csv("Alligators2.dat", sep="", header=TRUE)
N <- nrow(Alligators2)
head(Alligators2)
# lake gender size y1 y2 y3 y4 y5
# 1 1    1 1 7 1 0 0 5
# 2 1    1 0 4 0 0 1 2
# 3 1    0 1 16 3 2 2 3
# 4 1    0 0 3 0 1 2 3
# 5 2    1 1 2 2 0 0 1
# 6 2    1 0 13 7 6 0 0
attach(Alligators2)

formula <- cbind(y2,y3,y4,y5,y1) ~
size + factor(lake) + factor(gender) +
size:factor(lake) + size:factor(gender) + factor(lake):factor(gender)
fit <- vglm(formula,
         family=multinomial,
         data=Alligators2) # fish=1 is baseline category > summary(fit)

# fit2 <- vglm(formula = cbind(y2,y3,y4,y5,y1) ~ size,
#          family=multinomial,
#          data=Alligators2) # fish=1 is baseline category > summary(fit)
# anova(fit1, fit2, test="LRT", type="I")

step_fit <- step4vglm(fit, k=log(N))
# Step: AIC=211.54
# cbind(y2, y3, y4, y5, y1) ~ size + factor(lake)
#
#             Df Deviance AIC
# <none>          59.709 211.54
# - size       4 73.566 214.31
# - factor(lake) 12 109.327 227.89
summary(step_fit)
# Call:
# vglm(formula = cbind(y2, y3, y4, y5, y1) ~ size + factor(lake),
#    family = multinomial, data = Alligators2)
#
# Coefficients:
#             Estimate Std. Error z value Pr(>|z|)
# (Intercept):1 -3.0655 0.6485 -4.727 2.28e-06 ***
# (Intercept):2 -2.0646 0.7296 -2.830 0.00466 **
# (Intercept):3 -1.6166 0.6401 -2.526 0.01155 *
# (Intercept):4 -0.9041 0.4666 -1.938 0.05266 .
# size:1          1.2966 0.4228 3.067 0.00216 **
# size:2       -0.3604 0.6458 -0.558 0.57680
# size:3       -1.6166 0.6444 -0.405 0.68563
# size:4          0.0956 0.4591 0.208 0.83507
# factor(lake)2:1 2.9195 0.7032 4.152 3.30e-05 ***
# factor(lake)2:2 1.1457 0.8360 1.370 0.17055
# factor(lake)2:3 -1.2545 1.1961 -1.049 0.29428

6.18, 6.20, 6.21, 6.22, 6.23 Date set found in this website: Book...