# Analysis of the 2008 WFFC data # Concentrates on density estimation analyses # T. W. Yee (2013) # Note: argument "equalsd" used to be "ESD". # ====================================================================== # 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) # Convert lengths to cm wffc <- transform(wffc, lengthcm=length/10) # Subsets; one data frame for each sector waihou <- subset(wffc, water=="Waihou") waimak <- subset(wffc, water=="Waimakariri") whang <- subset(wffc, water=="Whanganui") otam <- subset(wffc, water=="Otamangakau") roto <- subset(wffc, water=="Rotoaira") # ====================================================================== # Run this code to set up constants. myxlab <- "Length (cm)" myylab <- "Density" epsilon <-- 0.01 mybreaks <- seq(18-epsilon, 72-epsilon, by=3) # by 6 is too coarse mybreaks <- seq(18-epsilon, 70-epsilon, by=3) # by 6 is too coarse myylim <- c(0, 0.15) myxlim <- c(min(mybreaks), max(mybreaks)) xx <- seq(myxlim[1], myxlim[2], by=0.1) # cm mylwd <- 1.2 mycol <- "black" sectorNames <- c("Whanganui"="Whanganui", "Otamangakau"="Otamangakau", "Rotoaira"="Rotoaira", "Waihou"="Waihou", "Waimakariri"="Waimakariri") # ====================================================================== # Fit mix2normal() and gamma2() to the sectors. # Later overlay these onto the histograms. # (i) Waihou ------------------------------------------------------ dewaihou <- transform(waihou, y0.18cm = (length-179.99)/10) with(dewaihou, hist(y0.18cm)) fit1.waihou = vglm(y0.18cm ~ 1, fam=gamma2.ab, dewaihou, trace=TRUE) fit2.waihou = vglm(y0.18cm ~ 1, fam=gamma2, dewaihou, trace=TRUE) ifit1.waihou = vglm(y0.18cm ~ 1, gamma2.ab(lrate="identity",lshape="identity", irate=0.15, ishape=0.73), data=dewaihou, trace=TRUE) summary(ifit1.waihou) coef(fit1.waihou, mat=TRUE) Coef(fit1.waihou) coef(fit2.waihou, mat=TRUE) Coef(fit2.waihou) max(abs(fitted(fit1.waihou) - fitted(fit2.waihou))) # 0 == same! pdf.waihou = pdf1.waihou = dgamma(xx-min(xx)+epsilon, rate=Coef(fit1.waihou)["rate"], shape=Coef(fit1.waihou)["shape"]) pdf2.waihou = dgamma(xx-min(xx)+epsilon, scale=Coef(fit2.waihou)["mu"] / Coef(fit2.waihou)["shape"], shape=Coef(fit2.waihou)["shape"]) # Plot the results with(waihou, hist(lengthcm, xlab="", breaks=mybreaks, prob=TRUE, ylim=myylim, border="blue", ylab=myylab, main=sectorNames["Waihou"], lwd=2)) lines(xx, pdf1.waihou, col=mycol, lwd=mylwd) # lines(xx, pdf2.waihou, col=mycol, lwd=mylwd) max(abs(pdf1.waihou - pdf2.waihou)) # 0 == same! # (ii) Waimak ------------------------------------------------------ dewaimak = transform(waimak, y0.18cm = (length-179.99)/10) with(dewaimak, hist(y0.18cm)) fit.waimak = vglm(y0.18cm ~ 1, fam=gamma2.ab, dewaimak, trace=TRUE) ifit.waimak = vglm(y0.18cm ~ 1, fam=gamma2.ab(lrate="identity",lshape="identity", irate=0.13, ishape=0.66), dewaimak, trace=TRUE) coef(fit.waimak, mat=TRUE) Coef(fit.waimak) summary(ifit.waimak) pdf.waimak = dgamma(xx-min(xx)+epsilon, rate=Coef(fit.waimak)["rate"], shape=Coef(fit.waimak)["shape"]) # Plot the results with(waimak, hist(lengthcm, xlab="", breaks=mybreaks, prob=TRUE, ylim=myylim, border="blue", ylab=myylab, main=sectorNames["Waimak"], lwd=2)) lines(xx, pdf.waimak, col=mycol, lwd=mylwd) # ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, # LRT to see if Waihou and Waimakariri are equal wrt the gamma2 parameters # Results shows that the two parameters can be considered equal. bothwaiwai = data.frame(rbind(dewaihou, dewaimak)) bothwaiwai[1:5,] fit.gamma2.waiwai = vglm(y0.18cm ~ 1, gamma2.ab, bothwaiwai, trace=TRUE) fit.gamma2.waiwai coef(fit1.waihou) coef(fit.waimak) coef(fit.gamma2.waiwai) Coef(fit1.waihou) Coef(fit.waimak) Coef(fit.gamma2.waiwai) # p-value pchisq(2*(logLik(fit1.waihou) + logLik(fit.waimak) - logLik(fit.gamma2.waiwai)), df=length(coef(fit1.waihou)) + length(coef(fit.waimak)) - length(coef(fit.gamma2.waiwai)), lower.tail=FALSE) length(fit1.waihou@y) + length(fit.waimak@y) + length(fit.gamma2.waiwai@y) # ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, # (iii) Whanganui ------------------------------------------------------ fit.whang = vglm(lengthcm ~ 1, fam=mix2normal(imu1=20, imu2=38, iphi=0.7, isd1=4, isd2=10, equalsd=FALSE), data=whang, trace=TRUE) # Plot the results with(whang, hist(lengthcm, xlab="", breaks=mybreaks, prob=TRUE, ylim=myylim, border="blue", ylab=myylab, main=sectorNames["Whanganui"], lwd=2)) Cfit.whang = Coef(fit.whang) phi.est = Cfit.whang["phi"] mu1.est = Cfit.whang["mu1"] mu2.est = Cfit.whang["mu2"] sd1.est = Cfit.whang["sd1"] sd2.est = Cfit.whang["sd2"] pdf.whang = phi.est * dnorm(xx, mu1.est, sd1.est) + (1-phi.est) * dnorm(xx, mu2.est, sd2.est) lines(xx, pdf.whang, col=mycol, lwd=mylwd) # (iv) Otamangakau ------------------------------------------------------ fit.otam = vglm(lengthcm ~ 1, mix2normal(imu1=20, imu2=50, iphi=0.55, isd1=5, isd2=14, equalsd=FALSE), data=otam, trace=TRUE) # Plot the results with(otam, hist(lengthcm, xlab="", breaks=mybreaks, prob=TRUE, ylim=myylim, border="blue", ylab=myylab, main=sectorNames["Otamangakau"], lwd=2)) Cfit.otam = Coef(fit.otam) phi.est = Cfit.otam["phi"] mu1.est = Cfit.otam["mu1"] mu2.est = Cfit.otam["mu2"] sd1.est = Cfit.otam["sd1"] sd2.est = Cfit.otam["sd2"] pdf.otam = phi.est * dnorm(xx, mu1.est, sd1.est) + (1-phi.est) * dnorm(xx, mu2.est, sd2.est) lines(xx, pdf.otam, col=mycol, lwd=mylwd) # (v) Rotoaira ------------------------------------------------------ fit.roto = vglm(lengthcm ~ 1, mix2normal(imu2=45, imu1=25, iphi=0.1, isd1=5, isd2=7,equalsd=TRUE), data=roto, trace=TRUE) ifit.roto = vglm(lengthcm ~ 1, mix2normal(imu2=45, imu1=25, lphi="identity", lsd="identity", iphi=0.1, isd1=5, isd2=7,equalsd=TRUE), data=roto, trace=TRUE) # Plot the results with(roto, hist(lengthcm, xlab="", breaks=mybreaks, prob=TRUE, ylim=myylim, border="blue", ylab=myylab, main=sectorNames["Rotoaira"], lwd=2)) phi.est = logit(coef(fit.roto)[1], inverse=TRUE) mu1.est = Coef(fit.roto, mat=TRUE)[1,"mu1"] mu2.est = Coef(fit.roto, mat=TRUE)[1,"mu2"] sd.est = exp(coef(fit.roto)[3]) pdf.roto = phi.est * dnorm(xx, mu1.est, sd.est) + (1-phi.est) * dnorm(xx, mu2.est, sd.est) lines(xx, pdf.roto, col=mycol, lwd=mylwd) c(phi.est,mu1.est ,mu2.est, sd.est) coef(ifit.roto, mat=TRUE) constraints(ifit.roto) # Look at the SEs round(sqrt(diag(vcov(fit.whang, untransform=TRUE))), dig=3) round(sqrt(diag(vcov(fit.otam, untransform=TRUE))), dig=3) round(sqrt(diag(vcov(ifit.roto, untransform=FALSE))), dig=3) # ====================================================================== # Plot all the histograms side by side # Nb. Run the above code above first # Nb. the variable "overlaymydensities" is logical; choose T or F overlaymydensities = T # par(las=1, mfrow=c(2,3), mar=c(4,3,1.2,1)+0.1, mgp=c(0,1,0)) with(waihou, hist(length/10, xlab="", breaks=mybreaks, prob=TRUE, ylim=myylim, border="blue", ylab="", main=sectorNames["Waihou"])) if(overlaymydensities) lines(xx, pdf.waihou, col=mycol, lwd=mylwd) with(otam, hist(length/10, xlab="", breaks=mybreaks, prob=TRUE, ylim=myylim, border="blue", ylab="", main=sectorNames["Otamangakau"])) if(overlaymydensities) lines(xx, pdf.otam, col=mycol, lwd=mylwd) with(waimak, hist(length/10, xlab="", breaks=mybreaks, prob=TRUE, ylim=myylim, border="blue", ylab="", main=sectorNames["Waimakariri"])) if(overlaymydensities) lines(xx, pdf.waimak, col=mycol, lwd=mylwd) with(roto, hist(length/10, xlab="", breaks=mybreaks, prob=TRUE, ylim=myylim, border="blue", ylab="", main=sectorNames["Rotoaira"])) if(overlaymydensities) lines(xx, pdf.roto, col=mycol, lwd=mylwd) with(whang, hist(length/10, xlab="", breaks=mybreaks, prob=TRUE, ylim=myylim, border="blue", ylab="", main=sectorNames["Whanganui"])) if(overlaymydensities) lines(xx, pdf.whang, col=mycol, lwd=mylwd) mtext(text=myxlab, side=1, line= 0.0, outer=TRUE, las=0, cex=1.3) mtext(text=myylab, side=2, line= 0.0, outer=TRUE, las=0, cex=1.3) # ======================================================================