mistnet-ordination

mistnet-ordination.R

# Install the most recent version of mistnet (if needed) using the devtools package
# devtools::install_github("davharris/mistnet")
library(mistnet)

# Grab a small data set from another package
library(ade4)
data(aviurba)

# Reduce the fauna table to ones and zeros
y = as.matrix(aviurba$fau > 0)
storage.mode(y) = "integer"

# Build a mistnet model
net = mistnet(
  x = matrix(NA, nrow = nrow(y), ncol = 0),
  y = y,
  layer.definitions = list(
    defineLayer(
      nonlinearity = rectify.nonlinearity(),
      size = 10,
      prior = gaussian.prior(mean = 0, var = 100)
    ),
    defineLayer(
      nonlinearity = sigmoid.nonlinearity(),
      size = ncol(y),
      prior = gaussian.prior(mean = 0, var = 100)
    )
  ),
  loss = bernoulliLoss(),
  updater = adagrad.updater(learning.rate = 0.1),
  n.importance.samples = 50
)

# Currently, mistnet does not initialize the coefficients automatically.
# This gets it started with nonzero values.
for(layer in net$layers){
  layer$coefficients[ , ] = rnorm(length(layer$coefficients), sd = .1)
  
}

# initialize the intercepts
net$layers[[1]]$biases[] = .1
net$layers[[2]]$biases[] = qlogis(colMeans(y))

starting.coefs = net$layers[[2]]$coefficients

# Fit the model
net$fit(1000)

# Confirm that the coefficients have actually moved
plot(net$layers[[2]]$coefficients ~ starting.coefs, asp = 1)
abline(0,1)

# Plot the resuls ---------------------------------------------------------

# Responses to the transformed latent variables
plot(prcomp(net$layers[[2]]$coefficients))
biplot(prcomp(net$layers[[2]]$coefficients))
lattice::levelplot(net$layers[[2]]$coefficients)

# Plot two random species' occurrence probabilities against each other:
spp = sample.int(ncol(y), 2)

plot(
  c(net$layers[[2]]$outputs[ , spp[1], ]), 
  c(net$layers[[2]]$outputs[ , spp[2] , ]),
  xlim = c(0, 1),
  ylim = c(0, 1),
  asp = 1
)
abline(v = c(0, 1), h = c(0, 1))