# Analysis of the 2008 WFFC data # Concentrates on the log-linear modelling via ZIP, ZAP, ZINB, ZANB, etc. # T. W. Yee (2013) # *** Applied to all sectors, not just to the rivers .*** # ====================================================================== # Read the raw data in library("VGAM") # Requires Version 0.9-2 or later library("VGAMdata") # Requires Version 0.9-2 or later # Check the input stage is ok summary(wffc) dim(wffc) head(wffc.nc) # ====================================================================== # Use quasipoissonff() to fit a log-linear model. # Nb. "fnc" stands for factor version of "wffc.nc" fnc <- transform(wffc.nc, finame = factor(iname), fsector = factor(sector), day = ceiling(session / 2), fday = factor(ceiling(session / 2)), mornaft = 1 - (session %% 2), fbeatboat = factor(beatboat)) head(fnc) nrow(fnc) # Delete people who fished less than 5 sessions? fnc <- fnc[with(fnc, !is.element(comid, c(99,72,80,93,45,71,97,78))),] nrow(fnc) # Should be 455 # ----------------------------------------------- # Notes: # 1. Whanganui is the baseline river. # 2. Could use fday or day. # 3. Just to get phi-hat fit.qp <- vglm(numbers ~ fsector + mornaft + fday + finame, quasipoissonff, data=fnc, trace=TRUE) summary(fit.qp) # Get phi-hat. AIC(fit.qp) length(fit.qp@y) length(coef(fit.qp)) # Ratio of # of observations to # of parameters (> 10 is good?) length(fit.qp@y) / length(coef(fit.qp)) # ====================================================================== # Fit a negbinomial model. fit.nb <- vglm(numbers ~ fsector + mornaft + fday + finame, negbinomial, data = fnc, trace = TRUE) coef(fit.nb, matrix = TRUE) summary(fit.nb) AIC(fit.nb) length(fit.nb@y) length(coef(fit.nb)) # Ratio of # of observations to # of parameters (> 10 is good?) length(fit.nb@y) / length(coef(fit.nb)) # ===================================================================== # Try fitting a ZANB, ZINB etc. # Uncomment out the appropriate line to select between ZIP, ZAP, ZINB, ZANB, etc. fit.zaipnb <- vglm(numbers ~ fsector + mornaft + fday + finame, # Carefully think!! fam = zipoisson(zero = 1), # fam = zapoisson(zero = 1), # fam = zinegbinomial(zero = c(1, 3)), # fam = zanegbinomial(zero = c(1, 3)), data = fnc, maxit = 34, trace = TRUE) head(coef(fit.zaipnb, matrix = TRUE), 8) # summary(fit.zaipnb) # AIC(fit.zaipnb) # length(fit.zaipnb@y) length(coef(fit.zaipnb)) # Ratio of # of observations to # of parameters (> 10 is good?) length(fit.zaipnb@y) / length(coef(fit.zaipnb)) # ===================================================================== # LRT pchisq(2*(logLik(fit.nb) - logLik(fit.qp)), df = length(coef(fit.nb)) - length(coef(fit.qp)), lower.tail = FALSE) length(coef(fit.qp)) length(coef(fit.nb)) # ======================================================================