# NOT RUN { # Example 1; Fit an unequal tolerances model to the hunting spiders data hspider[,1:6] <- scale(hspider[,1:6]) # Standardized environmental variables set.seed(1234) # For reproducibility of the results p1ut <- cqo(cbind(Alopacce, Alopcune, Alopfabr, Arctlute, Arctperi, Auloalbi, Pardlugu, Pardmont, Pardnigr, Pardpull, Trocterr, Zoraspin) ~ WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux, fam = poissonff, data = hspider, Crow1positive = FALSE, eq.tol = FALSE) sort(deviance(p1ut, history = TRUE)) # A history of all the iterations if (deviance(p1ut) > 1177) warning("suboptimal fit obtained") S <- ncol(depvar(p1ut)) # Number of species clr <- (1:(S+1))[-7] # Omits yellow lvplot(p1ut, y = TRUE, lcol = clr, pch = 1:S, pcol = clr, las = 1) # Ordination diagram legend("topright", leg = colnames(depvar(p1ut)), col = clr, pch = 1:S, merge = TRUE, bty = "n", lty = 1:S, lwd = 2) (cp <- Coef(p1ut)) (a <- latvar(cp)[cp@latvar.order]) # Ordered site scores along the gradient # Names of the ordered sites along the gradient: rownames(latvar(cp))[cp@latvar.order] (aa <- Opt(cp)[, cp@Optimum.order]) # Ordered optimums along the gradient aa <- aa[!is.na(aa)] # Delete the species that is not unimodal names(aa) # Names of the ordered optimums along the gradient trplot(p1ut, which.species = 1:3, log = "xy", type = "b", lty = 1, lwd = 2, col = c("blue","red","green"), label = TRUE) -> ii # Trajectory plot legend(0.00005, 0.3, paste(ii$species[, 1], ii$species[, 2], sep = " and "), lwd = 2, lty = 1, col = c("blue", "red", "green")) abline(a = 0, b = 1, lty = "dashed") S <- ncol(depvar(p1ut)) # Number of species clr <- (1:(S+1))[-7] # Omits yellow persp(p1ut, col = clr, label = TRUE, las = 1) # Perspective plot # Example 2; Fit an equal tolerances model. Less numerically fraught. set.seed(1234) p1et <- 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) sort(deviance(p1et, history = TRUE)) # A history of all the iterations if (deviance(p1et) > 1586) warning("suboptimal fit obtained") S <- ncol(depvar(p1et)) # Number of species clr <- (1:(S+1))[-7] # Omits yellow persp(p1et, col = clr, label = TRUE, las = 1) # Example 3: A rank-2 equal tolerances CQO model with Poisson data # This example is numerically fraught... need I.toler = TRUE too. set.seed(555) p2 <- 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, I.toler = TRUE, Rank = 2, Bestof = 3, isd.latvar = c(2.1, 0.9)) sort(deviance(p2, history = TRUE)) # A history of all the iterations if (deviance(p2) > 1127) warning("suboptimal fit obtained") lvplot(p2, ellips = FALSE, label = TRUE, xlim = c(-3,4), C = TRUE, Ccol = "brown", sites = TRUE, scol = "grey", pcol = "blue", pch = "+", chull = TRUE, ccol = "grey") # Example 4: species packing model with presence/absence data set.seed(2345) n <- 200; p <- 5; S <- 5 mydata <- rcqo(n, p, S, fam = "binomial", hi.abundance = 4, eq.tol = TRUE, es.opt = TRUE, eq.max = TRUE) myform <- attr(mydata, "formula") set.seed(1234) b1et <- cqo(myform, binomialff(multiple.responses = TRUE, link = "cloglog"), data = mydata) sort(deviance(b1et, history = TRUE)) # A history of all the iterations lvplot(b1et, y = TRUE, lcol = 1:S, pch = 1:S, pcol = 1:S, las = 1) Coef(b1et) # Compare the fitted model with the 'truth' cbind(truth = attr(mydata, "concoefficients"), fitted = concoef(b1et)) # Example 5: Plot the deviance residuals for diagnostic purposes set.seed(1234) p1et <- cqo(cbind(Alopacce, Alopcune, Alopfabr, Arctlute, Arctperi, Auloalbi, Pardlugu, Pardmont, Pardnigr, Pardpull, Trocterr, Zoraspin) ~ WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux, poissonff, data = hspider, eq.tol = TRUE, trace = FALSE) sort(deviance(p1et, history = TRUE)) # A history of all the iterations if (deviance(p1et) > 1586) warning("suboptimal fit obtained") S <- ncol(depvar(p1et)) par(mfrow = c(3, 4)) for (ii in 1:S) { tempdata <- data.frame(latvar1 = c(latvar(p1et)), sppCounts = depvar(p1et)[, ii]) tempdata <- transform(tempdata, myOffset = -0.5 * latvar1^2) # For species ii, refit the model to get the deviance residuals fit1 <- vglm(sppCounts ~ offset(myOffset) + latvar1, poissonff, data = tempdata, trace = FALSE) # For checking: this should be 0 # print("max(abs(c(Coef(p1et)@B1[1,ii],Coef(p1et)@A[ii,1])-coef(fit1)))") # print( max(abs(c(Coef(p1et)@B1[1,ii],Coef(p1et)@A[ii,1])-coef(fit1))) ) # Plot the deviance residuals devresid <- resid(fit1, type = "deviance") predvalues <- predict(fit1) + fit1@offset ooo <- with(tempdata, order(latvar1)) plot(predvalues + devresid ~ latvar1, data = tempdata, col = "red", xlab = "latvar1", ylab = "", main = colnames(depvar(p1et))[ii]) with(tempdata, lines(latvar1[ooo], predvalues[ooo], col = "blue")) } # }
Run the code above in your browser using DataCamp Workspace