# Structured/correlated row effects/LVs, cover data models, phylogenetic model # Section 4.2. in the paper. # Note, that minor differences are to be expected, due to randomness. set.seed(2025) library(gllvm) data("kelpforest") Xenv <- kelpforest$X SPinfo <- kelpforest$SPinfo Yabund <- kelpforest$Y # include species that have been observed on at least 10 units, sum up the rest Yinv <- Yabund[SPinfo$GROUP == "INVERT"] Yinv10 <- Yinv[,colSums(Yinv>0, na.rm = TRUE)>=10] Yinv10$RARE = rowSums(Yinv[,(colSums(Yinv>0, na.rm=TRUE)<10)], na.rm=TRUE) SPinfo10 = SPinfo[SPinfo$SP_CODE %in% colnames(Yinv10),] SPinfo10 = rbind(SPinfo10,c("RARE","INVERT", rep(NA, ncol(SPinfo10)-2))) SPinfo10 # remove empty rows: rem <- which(rowSums(Yinv10)==0) Yinv10 <- Yinv10[-rem,] Xenv <- Xenv[-rem,] SPinfo <- SPinfo[-rem,] # covariates: par(mfrow=c(1,2)) hist(Xenv$KELP_FRONDS, main = "KELP FRONDS") Xenv$logKELP_FRONDS = log(Xenv$KELP_FRONDS+1) Xenv$logKELP_FRONDSsc = scale(Xenv$logKELP_FRONDS) Xenv$PERCENT_ROCKYsc = scale(Xenv$PERCENT_ROCKY) hist(Xenv$logKELP_FRONDSsc, main = "log of KELP FRONDS") dev.off() # 4.2.1.-4.2.2.: only use sessile invertebrates Ysess <- Yinv10 X <- Xenv Xenv$TRANSECT = as.factor(Xenv$TRANSECT) # Structured row effect: setMap = list(zeta=c(1:ncol(Ysess),rep(NA,ncol(Ysess)))) shapeForm = rep(1, ncol(Ysess)) ftStrucRow <- gllvm(y = Ysess, X = Xenv, family = "orderedBeta", formula = ~ logKELP_FRONDSsc + PERCENT_ROCKYsc, studyDesign = Xenv[,c("SITE", "TRANSECT", "YEAR")], row.eff = ~ YEAR + (1|SITE/TRANSECT), num.lv = 0, method = "EVA", link = "logit", disp.formula = shapeForm, setMap = setMap, control.start = list(zetacutoff = c(0, 20))) # AR1 latent variables: ftCorrLV <- gllvm(Ysess, Xenv, family = "orderedBeta", formula = ~ logKELP_FRONDSsc + PERCENT_ROCKYsc, studyDesign = Xenv[,c("SITE", "TRANSECT", "YEAR")], row.eff = ~ SITE/TRANSECT, starting.val="zero", num.lv = 2, disp.formula = shapeForm, method="EVA", lvCor = ~corAR1(1|YEAR), link="logit", zetacutoff = c(0, 20), setMap = setMap) ftCorrLV$params$rho.lv # Variance partitioning, reproduces Fig 3. varPart <- varPartitioning(ftCorrLV, groupnames=c("Giant kelp frond density (log)", "Seabed rockiness (%)", "LV1", "LV2", "Site-transect fixed effect")) plotVP(varPart, args.legend = list(cex=0.9), col=hcl.colors(5, "TealRose"), xlab="Species", las=2, cex.names=0.7) # Fourth corner model with trait for group, i.e., INVERTEBRAE or ALGAE Xenv <- kelpforest$X SPinfo <- kelpforest$SPinfo Xenv$logKELP_FRONDS = log(Xenv$KELP_FRONDS+1) Xenv$logKELP_FRONDSsc = scale(Xenv$logKELP_FRONDS) Xenv$PERCENT_ROCKYsc = scale(Xenv$PERCENT_ROCKY) Y10 <- Yabund Y10 <- Yabund[SPinfo$GROUP != "PLANT"] Y10 <- Y10[,colSums(Y10>0, na.rm = TRUE)>=10] NAs <- which(is.na(colSums(Y10))) NAs Y10 <- Y10[,-NAs] emptyR <- which(rowSums(Y10)==0) Y10 <- Y10[-emptyR,] Xenv <- Xenv[-emptyR,] SPinfo10 = SPinfo[SPinfo$SP_CODE %in% colnames(Y10),] SPinfo10$GROUP <- droplevels(SPinfo10$GROUP) table(SPinfo10$GROUP) Ykelp = Y10 Traits = SPinfo10[,1:2, drop=FALSE] #rownames(Traits) <- SPinfo10$SP_CODE setMap = list(zeta=c(1:ncol(Ykelp),rep(NA,ncol(Ykelp)))) shapeForm = ifelse(Traits$GROUP == "INVERT", 1, 2) # Fit 4th corner model: # ftFourthC <- gllvm(y = Ykelp, X = Xenv, TR = Traits, # formula = ~ logKELP_FRONDSsc + PERCENT_ROCKYsc + # (logKELP_FRONDSsc + PERCENT_ROCKYsc) : (GROUP) # #+ (0+logKELP_FRONDSsc + PERCENT_ROCKYsc|1) # , num.lv = 0, # disp.formula = shapeForm, starting.val = "zero", # family = "orderedBeta", link = "logit", method = "EVA", # zetacutoff = c(0,20), setMap = setMap) ftFourthC <- gllvm(y = Ykelp, X = Xenv, TR = Traits, formula = ~ logKELP_FRONDSsc + PERCENT_ROCKYsc + (logKELP_FRONDSsc + PERCENT_ROCKYsc) : (GROUP), num.lv = 0, randomX = ~logKELP_FRONDSsc + PERCENT_ROCKYsc, disp.formula = shapeForm, starting.val = "zero", beta0com=TRUE, optimizer="nlminb", #randomX.start="res", family = "orderedBeta", link = "logit", method = "EVA", zetacutoff = c(0,20), setMap = setMap) # phylogenetic model: #install.packages("ape") library(ape) taxa <- data.frame(Kingdom=factor(SPinfo10$TAXON_KINGDOM), Phylum=factor(SPinfo10$TAXON_PHYLUM), Class=factor(SPinfo10$TAXON_CLASS), Order=factor(SPinfo10$TAXON_ORDER), Family=factor(SPinfo10$TAXON_FAMILY), Genus=factor(SPinfo10$TAXON_GENUS), Species=factor(SPinfo10$SP_CODE)) # remove species with NA information taxa <- na.omit(taxa) Yphyl <- Y10[,colnames(Y10) %in% taxa$Species] SPphyl = SPinfo10[SPinfo10$SP_CODE %in% taxa$Species,,-3] dim(Yphyl) phyl <- ape::as.phylo(~Kingdom/Phylum/Class/Order/Family/Genus/Species, data=taxa) tree <- ape::compute.brlen(phyl) colMat <- ape::vcv(tree) dist <- ape::cophenetic.phylo(tree) allDists <- ape::dist.nodes(tree) order0 <- as.integer(names(sort(allDists[1:length(tree$tip.label), nrow(allDists)],decreasing = TRUE))) order <- tree$tip.label[order0] Trphyl <- SPphyl[,"GROUP", drop=FALSE] row.names(Trphyl) <- SPphyl[,"SP_CODE"] shapeForm = ifelse(Trphyl[order,,drop=FALSE]$GROUP == "INVERT", 1, 2) Xenv <- Xenv[-which(rowSums(Yphyl)==0),] Yphyl <- Yphyl[-which(rowSums(Yphyl)==0),] setMap = list(zeta=c(1:ncol(Yphyl),rep(NA,ncol(Yphyl)))) ftPhylo = gllvm(y = Yphyl[,order], X = Xenv, TR = Trphyl[order,, drop=FALSE], formula = ~(logKELP_FRONDSsc + PERCENT_ROCKYsc) + (logKELP_FRONDSsc + PERCENT_ROCKYsc) : (GROUP), randomX =~logKELP_FRONDSsc + PERCENT_ROCKYsc, num.lv = 0, colMat = list(colMat[order, order], dist = dist[order, order]), beta0com = TRUE, colMat.rho.struct = "term", nn.colMat = 10, disp.formula = shapeForm, starting.val = "zero", family = "orderedBeta", link = "logit", method = "EVA", zetacutoff = c(0, 20), setMap = setMap, optim.method = "L-BFGS-B", n.init = 3) summary(ftPhylo) # Phylogeny plot, Fig. 4 phyloplot(ftPhylo, tree, cex = 0.8, tick.length = 0.5, widths = c(.5,.1), heights = c(0.4,0.4))