set.seed(123); nn <- 200
cdata <- data.frame(x2 = rnorm(nn),  # Has mean 0 (needed when I.tol=TRUE)
                    x3 = rnorm(nn),  # Has mean 0 (needed when I.tol=TRUE)
                    x4 = rnorm(nn))  # Has mean 0 (needed when I.tol=TRUE)
cdata <- transform(cdata, latvar1 =  x2 + x3 - 2*x4,
                          latvar2 = -x2 + x3 + 0*x4)
# Nb. latvar2 is weakly correlated with latvar1
cdata <- transform(cdata,
            lambda1 = exp(6 - 0.5 * (latvar1-0)^2 - 0.5 * (latvar2-0)^2),
            lambda2 = exp(5 - 0.5 * (latvar1-1)^2 - 0.5 * (latvar2-1)^2),
            lambda3 = exp(5 - 0.5 * (latvar1+2)^2 - 0.5 * (latvar2-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, I.tolerances = 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(latvar(p2))  # A diagonal matrix, i.e., uncorrelated latent vars
# vvv var(latvar(p2, varI.latvar = TRUE))  # Identity matrix
# vvv Tol(p2)[, , 1:2]  # Identity matrix
# vvv Tol(p2, varI.latvar = TRUE)[, , 1:2]  # A diagonal matrixRun the code above in your browser using DataLab