options(prompt = "> ", continue=" ", useFancyQuotes = FALSE) options(width = 80) options(digits = 5) ## ------------------------------------------------------------------- library("VGAM") library("VGAMdata") ps.options(pointsize = 12, family = "Times") pdf.options(pointsize = 12, family = "Times") # rm(list = ls()) ## ------------------------------------------------------------------- hspider[, 1:6] <- scale(hspider[, 1:6]) # Standardized environmental variables set.seed(1234) # For reproducibility of the results p1ut.hs <- cqo(cbind(Alopacce, Alopcune, Alopfabr, Arctlute, Arctperi, Auloalbi, Pardlugu, Pardmont, Pardnigr, Pardpull, Trocterr, Zoraspin) ~ WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux, poissonff, data = hspider, eq.toler = FALSE, trace = FALSE) ## ------------------------------------------------------------------- Tol(p1ut.hs)[1, 1, ] ## ------------------------------------------------------------------- par(mfrow = c(2, 1), mar = c(5.2, 4.4, 1.1, 1.1) + 0.1, las = 1) set.seed(1234) # For reproducibility of the results p1ut.hs.2 <- cqo(cbind(Alopacce, Alopcune, Alopfabr, Arctlute, Arctperi, Auloalbi, Pardmont, Pardnigr, Pardpull, Trocterr, Zoraspin) ~ WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux, poissonff, data = hspider, Crow1positive = FALSE, eq.toler = FALSE, trace = FALSE) S <- ncol(depvar(p1ut.hs.2)) # Number of species clr <- (1:(S+1))[-7] # Omits yellow lvplot(p1ut.hs.2, main = "(a)", y = TRUE, lcol = clr, pch = 1:S, pcol = clr) legend("topright", leg = colnames(depvar(p1ut.hs.2)), col = clr, pch = 1:S, merge = TRUE, bty = "n", lty = 1:S, lwd = 2) persp(p1ut.hs.2, main = "(b)", col = clr, label = TRUE) # Perspective plot par(las = 0) # R default ## ------------------------------------------------------------------- round(concoef(p1ut.hs.2), digits = 2) ## ------------------------------------------------------------------- set.seed(1234) # For reproducibility of the results p1et.hs <- cqo(cbind(Alopacce, Alopcune, Alopfabr, Arctlute, Arctperi, Auloalbi, Pardlugu, Pardmont, Pardnigr, Pardpull, Trocterr, Zoraspin) ~ WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux, poissonff, data = hspider, Crow1positive = FALSE, eq.toler = TRUE, trace = FALSE) Coef(p1et.hs) ## ------------------------------------------------------------------- sort(deviance(p1et.hs, history = TRUE)) # A history of all the iterations ## ------------------------------------------------------------------- par(mfrow = c(1, 1), mar = c(5.2, 4.4, 1.1, 1.1) + 0.1, las = 1) set.seed(555) # For reproducibility of the results p2et.hs <- cqo(cbind(Alopacce, Alopcune, Alopfabr, Arctlute, Arctperi, Auloalbi, Pardlugu, Pardmont, Pardnigr, Pardpull, Trocterr, Zoraspin) ~ WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux, poissonff, data = hspider, Crow1positive = FALSE, Rank = 2, I.toler = TRUE, Bestof = 3, isd.latvar = c(2.1, 0.9)) if (deviance(p2et.hs) > 1127) warning("suboptimal fit obtained") persp(p2et.hs, xlim = c(-6, 5), ylim = c(-6, 3), theta = 120, phi = 20) ## ------------------------------------------------------------------- par(mfrow = c(2, 1), mar = c(5.2, 4.4, 1.1, 1.1) + 0.1, las = 1) lvplot(p1et.hs, main = "(a)", y = TRUE, lcol = clr, pch = 1:S, pcol = clr, las = 1) legend("topright", leg = colnames(depvar(p1et.hs)), col = clr, pch = 1:S, merge = TRUE, bty = "n", lty = 1:S, lwd = 2) persp(p1et.hs, main = "(b)", col = clr, label = TRUE, las = 1) # Perspective plot ## ------------------------------------------------------------------- par(mfrow = c(2, 1), mar = c(4.8, 3.9, 1.6, 0.1) + 0.1, las = 1) lvplot(p2et.hs, ellipse = 0.95, label = TRUE, xlim = c(-3, 5.7), C = FALSE, Ccol = "brown", sites = TRUE, scol = "gray50", pcol = "blue", pch = "+", chull = TRUE, ccol = "gray50", main = "(a)") lvplot(p2et.hs, ellipse = FALSE, label = TRUE, xlim = c(-3, 5.7), C = TRUE, Ccol = "brown", sites = TRUE, scol = "gray50", pcol = "blue", pch = "+", chull = TRUE, ccol = "gray50", main = "(b)") # 20140804; aspect ratio manually computed above to make the ellipses circular ## ------------------------------------------------------------------- (allvariance <- diag(latvar(p2et.hs))) sqrt(allvariance) ## ------------------------------------------------------------------- par(mfrow = c(1, 1), mar = c(4.1, 4, 0.2, 0.2) + 0.1) mycols <- c("blue", "orange", "green") myxlim <- c(0.5e-4, 20) tr.hs <- trplot(p1ut.hs.2, which.species = 1:3, log = "xy", type = "b", lty = 1, col = mycols, lwd = 2, label = TRUE, xlim = myxlim, ylim = myxlim) legend("left", lwd = 2, lty = 1, col = mycols, with(tr.hs, paste(species.names[, 1], species.names[, 2], sep = " and "))) abline(a = 0, b = 1, lty = "dashed", col = "gray50") # A useful reference line ## ------------------------------------------------------------------- par(mfrow = c(1, 1), mar = c(4, 4, 0.1, 0.9) + 0.1) set.seed(1234) # For reproducibility hspider[, 1:6] <- scale(hspider[, 1:6]) # Standardize the environmental variables N <- 2 # Number of sites to calibrate for test.index <- sample(nrow(hspider), size = N) # N randomly chosen sites p1et.hs.m2 <- cqo(cbind(Alopacce, Alopcune, Alopfabr, Arctlute, Arctperi, Auloalbi, Pardlugu, Pardmont, Pardnigr, Pardpull, Trocterr, Zoraspin) ~ WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux, quasipoissonff, data = hspider[-test.index, ], trace = FALSE) cal.p1et.hs <- calibrate(p1et.hs.m2, hspider[test.index, ]) # Add the calibrated sites scores to a perspective plot S <- ncol(depvar(p1et.hs)) # Number of species clr <- rep(c(1:6, 8), len = S) # Omits yellow persp(p1et.hs, las = 1, col = clr, label = TRUE) abline(v = cal.p1et.hs, col = 1:N, lty = 1:N, lwd = 1) # Calibrated values C.matrix <- concoef(p1et.hs) # Constrained coefficients C abline(v = as.matrix(hspider[test.index, rownames(C.matrix)]) %*% C.matrix, col = 1:N, lty = 1:N, lwd = 2) # From CQO fit ## ------------------------------------------------------------------- hspider[test.index, colnames(depvar(p1et.hs))] ## ------------------------------------------------------------------- trapO <- transform(trapO, sc.doy = scale(doy), sc.LevelTW = scale(LevelTW), sc.MinAT = scale(MinAT), sc.MaxAT = scale(MaxAT), f.year = factor(Year)) set.seed(123) p1ut.to <- cqo(cbind(BFTW, BMTW, RFTW, RMTW) ~ sc.doy + f.Year + sc.MinAT + sc.MaxAT + sc.LevelTW, eq.tolerances = FALSE, family = poissonff, trace = FALSE, Bestof = 10, Crow1positive = TRUE, data = trapO) ## ------------------------------------------------------------------- round(concoef(p1ut.to), digits = 2) # sc.doy is by far the largest coefficient ## ------------------------------------------------------------------- par(las = 1, mfrow = c(1, 1), mar = c(4.4, 4.0, 0.2, 1.5) + 0.1, cex = 1) persp(p1ut.to, col = 1:4, label = TRUE) ## ------------------------------------------------------------------- tol.ratio.trapO <- sqrt(range(Tol(p1ut.to)[1, 1, ])) tol.ratio.trapO <- max(tol.ratio.trapO) / min(tol.ratio.trapO) ## ------------------------------------------------------------------- Tol(p1ut.to)[1, 1, ] sqrt(range(Tol(p1ut.to)[1, 1, ])) ## ------------------------------------------------------------------- sort(round(Max(p1ut.to), digits = 1)) ## ------------------------------------------------------------------- set.seed(111) n <- 100; p <- 5; S <- 5 pdata <- rcqo(n, p, S, es.opt = FALSE, eq.max = FALSE, eq.toler = TRUE, sd.latvar = 3/4) true.nu <- attr(pdata, "latvar") # The 'truth' attr(pdata, "tolerances")[, 1] # The tolerances attr(pdata, "optimums")[, 1] # The optimums Y <- Select(pdata, "y") # Y matrix (n x S) uqo.rcim1 <- rcim(Y, Rank = 1) # Traditional parameterization uqo.grc1 <- grc(Y) # An equivalent simpler call ## ------------------------------------------------------------------- par(mfrow = c(2, 2), mar = c(5.2, 4.4, 1.1, 1.1) + 0.1, las = 1) max(abs( fitted(uqo.grc1) - fitted(uqo.rcim1))) # Should be 0 max(abs(predict(uqo.grc1) - predict(uqo.rcim1))) # Should be 0 # Plot 1 plot(attr(pdata, "optimums"), Coef(uqo.grc1)@A, col = 1:S, type = "p", main = "(a)") mylm <- lm(Coef(uqo.grc1)@A ~ attr(pdata, "optimums")) abline(coef = coef(mylm), col = "orange", lty = "dashed") # Plot 2 fill.val <- 0 # Choose this for the traditional parameterization plot(attr(pdata, "latvar"), c(fill.val, Coef(uqo.grc1)@C), las = 1, col = "blue", type = "p", main = "(b)") mylm <- lm(c(fill.val, Coef(uqo.grc1)@C) ~ attr(pdata, "latvar")) abline(coef = coef(mylm), col = "orange", lty = "dashed") ## ------------------------------------------------------------------- myform <- attr(pdata, "formula") p1ut <- cqo(myform, family = poissonff, eq.toler = FALSE, trace = FALSE, data = pdata) c1ut <- cqo(cbind(y1, y2, y3, y4, y5) ~ scale(latvar(uqo.rcim1)), family = poissonff, eq.toler = FALSE, trace = FALSE, data = pdata) lvplot(p1ut, lcol = 1:S, y = TRUE, pcol = 1:S, pch = 1:S, pcex = 0.5, main = "(c)") lvplot(c1ut, lcol = 1:S, y = TRUE, pcol = 1:S, pch = 1:S, pcex = 0.5, main = "(d)") ## ------------------------------------------------------------------- data.frame(BFTW = 10, BMTW = 5, RFTW = 20, RMTW = 10)