fhui28 / CBFM

Spatio-temporal joint species distribution modeling using community-level basis functions
13 stars 2 forks source link

Simple example of ordination for CBFM #7

Open fhui28 opened 3 years ago

fhui28 commented 3 years ago
rm(list = ls())
library(autoFRK)
library(FRK)
library(MASS)
library(tidyverse)
library(GGally)
library(Hmsc)
library(devtools)

# use HMSC data generation mechanism
source_url("https://raw.githubusercontent.com/hmsc-r/HMSC/master/vignettes/makedata.R")

# Generate data
tmp = makedata(ns=50, ny=500, spatial=TRUE)
all.data=tmp[[1]]
all.parameters=tmp[[2]]

L1 = all.parameters$L
Y2 = 1*(L1 + matrix(rnorm(n = nrow(all.data$Y)*ncol(all.data$Y)), ncol = ncol(all.data$Y)) > 0)
all.data$Y = Y2

# Set up spatial basis functions for CBFM -- Most practitioners will start here! 
num_basisfunctions <- 30 # Number of spatial basis functions to use
basisfunctions <- mrts(all.data$xy, num_basisfunctions) %>% 
as.matrix %>%
{.[,-(1)]} # Remove the first intercept column

# Fit CBFM 
# Note also that Hmsc generates and fits models assuming a probit link, but CBFM uses a logit link
fitcbfm <- CBFM(y = all.data$Y, formula_X = all.data$X.formula, data = all.data$X.data, 
                B_space = basisfunctions, family = binomial(), control = list(trace = 1))

# Compare ordinations from CBFM to actual LVs and loadings from data generation mechanism 
getords <- ordinate(fitcbfm, num_comp = 2)
data.frame(true = all.parameters$eta, cbfmestimates = getords$scores) %>% 
     ggpairs

data.frame(true = t(all.parameters$lambda), cbfmestimates = getords$loadings) %>% 
     ggpairs

But please also see the help file for ordinate.CBFM