hspider[, 1:6] <- scale(hspider[, 1:6]) # Good idea when I.tolerances = TRUE
set.seed(111)
r1 <- cqo(cbind(Alopacce, Alopcune, Alopfabr, Arctlute, Arctperi,
Auloalbi, Pardmont, Pardnigr, Pardpull, Trocterr) ~
WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux,
poissonff, data = hspider, trace = FALSE, I.tolerances = TRUE)
set.seed(111) # r2 below is an ill-conditioned model
r2 <- cqo(cbind(Alopacce, Alopcune, Alopfabr, Arctlute, Arctperi,
Auloalbi, Pardmont, Pardnigr, Pardpull, Trocterr) ~
WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux,
isd.lv = c(2.4, 1.0), Muxfactor = 3.0, trace = FALSE,
poissonff, data = hspider, Rank = 2, eq.tolerances = TRUE)
sort(r1@misc$deviance.Bestof) # A history of the fits
sort(r2@misc$deviance.Bestof) # A history of the fits
if (deviance(r2) > 857) stop("suboptimal fit obtained")
persp(r1, xlim = c(-6,5), col = 1:4, label = TRUE)
# Involves all species
persp(r2, xlim = c(-6,5), ylim = c(-4, 5), theta = 10, phi = 20, zlim = c(0, 220))
# Omit the two dominant species to see what is behind them
persp(r2, xlim = c(-6,5), ylim = c(-4, 5), theta = 10, phi = 20, zlim = c(0, 220),
which = (1:10)[-c(8, 10)]) # Use zlim to retain the original z-scale
Run the code above in your browser using DataLab