set.seed(123)
nn <- 200
cdata <- data.frame(x2 = rnorm(nn), # Has mean 0 (needed when ITol=TRUE)
x3 = rnorm(nn), # Has mean 0 (needed when ITol=TRUE)
x4 = rnorm(nn)) # Has mean 0 (needed when ITol=TRUE)
cdata <- transform(cdata, lv1 = x2 + x3 - 2*x4,
lv2 = -x2 + x3 + 0*x4)
# Nb. lv2 is weakly correlated with lv1
cdata <- transform(cdata, lambda1 = exp(6 - 0.5 * (lv1-0)^2 - 0.5 * (lv2-0)^2),
lambda2 = exp(5 - 0.5 * (lv1-1)^2 - 0.5 * (lv2-1)^2),
lambda3 = exp(5 - 0.5 * (lv1+2)^2 - 0.5 * (lv2-0)^2))
cdata <- transform(cdata, spp1 = rpois(nn, lambda1),
spp2 = rpois(nn, lambda2),
spp3 = rpois(nn, lambda3))
set.seed(111)
# vvv p2 <- cqo(cbind(spp1,spp2,spp3) ~ x2 + x3 + x4, poissonff,
# vvv data = cdata,
# vvv Rank=2, ITolerances=TRUE,
# vvv Crow1positive=c(TRUE,FALSE)) # deviance = 505.81
# vvv if (deviance(p2) > 506) stop("suboptimal fit obtained")
# vvv sort(p2@misc$deviance.Bestof) # A history of the fits
# vvv Coef(p2)
lvplot(p2, sites = TRUE, spch = "*", scol = "darkgreen", scex = 1.5,
chull = TRUE, label = TRUE, Absolute = TRUE, ellipse = 140,
adj = -0.5, pcol = "blue", pcex = 1.3, las = 1,
C = TRUE, Cadj = c(-.3,-.3,1), Clwd = 2, Ccex = 1.4, Ccol = "red",
main = paste("Contours at Abundance = 140 with",
"convex hull of the site scores"))
# vvv var(lv(p2)) # A diagonal matrix, i.e., uncorrelated latent variables
# vvv var(lv(p2, varlvI = TRUE)) # Identity matrix
# vvv Tol(p2)[,,1:2] # Identity matrix
# vvv Tol(p2, varlvI = TRUE)[,,1:2] # A diagonal matrix
Run the code above in your browser using DataLab