library("cotram")
data("spiders", package = "cotram")
### for illustration only
OR <- 1 ### order of transformation function
### OR = 1 means log-linear, use OR ~ 6
M <- 100 ### number of Halton sequences, seem sufficient here
fastopt <- mltoptim(abstol = 1e-3, reltol = 1e-3)
### faster convergence
## fit conditional marginal count transformation models
## one for each species
m_PF <- cotram(Pardosa_ferruginea ~ Elevation + Canopy_openess,
data = spiders, method = "probit", order = OR)
m_HL <- cotram(Harpactea_lepida ~ Elevation + Canopy_openess,
data = spiders, method = "probit", order = OR)
m_CC <- cotram(Callobius_claustrarius ~ Elevation + Canopy_openess,
data = spiders, method = "probit", order = OR)
m_CT <- cotram(Coelotes_terrestris ~ Elevation + Canopy_openess,
data = spiders, method = "probit", order = OR)
m_PL <- cotram(Pardosa_lugubris ~ Elevation + Canopy_openess,
data = spiders, method = "probit", order = OR)
m_PR <- cotram(Pardosa_riparia ~ Elevation + Canopy_openess,
data = spiders, method = "probit", order = OR)
### fit dependence parameters
mm <- mcotram(m_PF, m_HL, m_CC, m_CT, m_PL, m_PR, data = spiders,
M = M, scale = TRUE, optim = fastopt)
logLik(mm)
### Kendall's tau: Dependence of species after accounting
### for elevation and canopy openess in marginal models
coef(mm, type = "Kendall")
# \donttest{
### regress dependencies on elevation and canopy openess
mmc <- mcotram(m_PF, m_HL, m_CC, m_CT, m_PL, m_PR, data = spiders,
formula = ~ Elevation + Canopy_openess, M = M,
scale = TRUE, optim = fastopt)
logLik(mmc)
### weak evidence for such effects
pchisq(2 * (logLik(mmc) - logLik(mm)), df = 30, lower.tail = FALSE)
### plot Kendall's tau for different elevations / openess levels
nd <- expand.grid(Elevation = 80:120 * 10, Canopy_openess = 1:10 * 10)
KD <- Lower_tri(coef(mmc, newdata = nd, type = "Kendall"))
f <- factor(rownames(KD))
nd <- cbind(f = rep(f, nrow(nd)), nd[rep(1:nrow(nd), each = nlevels(f)),])
nd$KD <- c(KD)
if (require("lattice"))
contourplot(KD ~ Elevation + Canopy_openess | f, data = nd,
cuts = 18, xlab = "Elevation", ylab = "Canopy openess")
### for example:
### => constant negative dependence of Pardosa_lugubris and Coelotes_terrestris
### => weak dependence of Harpactea_lepida and Pardosa_ferruginea
### for low elevations, negative dependence increasing with elevation
# }
Run the code above in your browser using DataLab