# Analysis of the 2008 WFFC data # Concentrates on fitting extreme value models and quantile regression. # T. W. Yee (2013) # Fits an extreme value model and quantile regression to the data. # ====================================================================== # Read the raw data in library("VGAM") # Requires Version 0.9-2 or later library("VGAMdata") # Requires Version 0.9-2 or later data(wffc) summary(wffc) dim(wffc) wffc.i = wffc.indiv waihou = subset(wffc, water=="Waihou") # ====================================================================== # Fit a GEV model to the Waihou only # Overwrite 'wffc.i' by the values only pertaining to Waihou wffc.i[["placing"]] = wffc.i[["points"]] = wffc.i[["individual"]] = wffc.i[["approxmeanlength"]] = NULL # Replace wffc.i[["longest"]] by the longest only from the small streams only unames = unique(wffc.i[["iname"]]) for(ii in unames) { iptr = (1:nrow(wffc.i))[wffc.i[["iname"]] == ii] index1 = with(waihou, iname == ii) if(any(index1)) { smalldf = waihou[index1,] wffc.i[iptr, "longest"] = max(smalldf$length) wffc.i[iptr, "nooffish"] = nrow(smalldf) } else { wffc.i[iptr, "longest"] = NA wffc.i[iptr, "nooffish"] = 0 } } wffc.i = na.omit(wffc.i) # Delete those who did not fish the Waihou wffc.i = transform(wffc.i, longestcm=longest/10) # Sort by nooffish (for plotting) wffc.i = wffc.i[sort.list(wffc.i$nooffish),] # ====================================================================== # Fit a GEV model to the Waihou only (recommended) standalone = FALSE # zz set FALSE for the manuscript, else TRUE myylab = "Longest length (cm)" par(mfrow=c(1+!standalone,1), mar=c(4.9, 4.1, 0.3, 0.3), oma=c(0.0, if(standalone) 0 else 1.8, 0, 0), las=1, xpd=TRUE) fit.gev = vglm(longestcm ~ ns(nooffish, df=4), gev(perc=c(50,75,90,95), zero=2:3, lshape="identity"), data=wffc.i, trace=TRUE) coef(fit.gev, matrix=TRUE) fitted(fit.gev)[1:4,] (sfit.gev <- summary(fit.gev)) # Is it a Gumbel distribution? 2*pnorm(-abs(sfit.gev@coef3["(Intercept):3","t value"])) # 2-sided pvalue LLL = ncol(fitted(fit.gev)) with(wffc.i, plot(nooffish, longestcm, col="blue", ylim=c(18,50), xlab=ifelse(standalone, "Number of fish caught", ""), las=1, ylab=ifelse(standalone, myylab, ""))) with(wffc.i, matlines(nooffish, fitted(fit.gev), lty=c(1:(LLL-2),(LLL-0):(LLL+1)), col=c(1:(LLL-1),LLL+0), lwd=2)) if(standalone) { dev.off() } else { text(x=2.0, y=49.0, "(a)", font=1) # "helvitica" } # ====================================================================== # Quantile regression based on LMS method; run after GEV model fitted. fit.lmsbcn = vglm(longestcm ~ ns(nooffish, df=4), lms.bcn(percentile=c(50,75,90,95), zero=c(1,3)), data=wffc.i, trace=TRUE) coef(fit.lmsbcn, matrix=TRUE) fitted(fit.lmsbcn)[1:4,] fit.lmsbcn@extra LLL = ncol(fitted(fit.lmsbcn)) with(wffc.i, plot(nooffish, longestcm, col="blue", log="", ylim=c(18,50), xlab="Number of fish caught", las=1, ylab=ifelse(standalone, myylab, ""))) with(wffc.i, matlines(nooffish, fitted(fit.lmsbcn), lty=c(1:(LLL-2),(LLL-0):(LLL+1)), col=c(1:(LLL-1),LLL+0), lwd=2)) if(standalone) { } else { text(x=2.0, y=49.0, "(b)", font=1) mtext(text=myylab, side=2, line= 0.0, outer=TRUE, las=0, cex=1.3) } # ======================================================================