Choosing a nonlinear model: GAM?

Question

I have the response variable that is an abundance matrix Y. for example, I have four columns with

• species 1

• species 2

• species 3

• and species 4

  abundances at each location (each location is separate row).

I have a matrix X which include a set of environmental parameters (e.g. Temperature and salinity).

I need to predict what will be my abundance matrix $Y_i$ (i.e. abundance of

• Species1=...,
• species 2=...,
• species 3=...,
• and species 4=...)

at a given set of environmental conditions ( for example T=20, S=35)

I saw a method that was using multivariate multiple regression approach:

model=lm(cbind(A,B) ~ c+d+e+f+g)

However, in my model, the relationship is not linear. so I was wondering whether you can suggest some other method.

I was thinking of doing the following steps:

1. Train a model to give a probability of occurrence of any combination of abundance matrix Y in a given set of environmental conditions
2. Use this model to predict what abundance matrix will look like for a new set of env conditions

I am trying it on on another dataset :

dataset(airquality)
fmla= cbind(Day, Month)~Ozone + Wind + Temp

I am not sure whether GAM model will be appropriate for this.

Answer 1

You'll need to make some assumptions about the nature of the response, but one option is to fit a multivariate normal model, which is possible for example in mgcv.

Using the airquality dataset mentioned in the question, we might have

data(airquality)
aq <- transform(airquality,
                date = as.Date(paste('1973', Month, Day, sep = '-')))
aq <- transform(aq, DoY = as.numeric(format(date, '%j')))

and we might assume (probably naively) that OzoneandTemp` are multivariate normal. In which case we can fit the model using

m <- gam(list(Ozone ~ s(DoY),
              Temp  ~ s(DoY)),
         data = aq, family = mvn(d=2))

where we specify the linear predictors for the two, in this instance, response variables. I chose to model them as smooth functions of the day of year and implicit in this model is the estimation of a covariance matrix.

> summary(m)

Family: Multivariate normal 
Link function: 

Formula:
Ozone ~ s(DoY)
Temp ~ s(DoY)

Parametric coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)    42.1293     2.6470   15.92   <2e-16 ***
(Intercept).1  77.8707     0.5704  136.53   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df Chi.sq  p-value    
s(DoY)   3.546  4.385  34.42 1.37e-06 ***
s.1(DoY) 5.711  6.832 149.04  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Deviance explained = 41.9%
-REML =  708.4  Scale est. = 1         n = 116

That covariance matrix can be extracted via

solve(crossprod(m$family$data$R))

which results in

> solve(crossprod(m$family$data$R))
         [,1]      [,2]
[1,] 812.7657 104.67443
[2,] 104.6744  37.73835

or as a correlation matrix

> cov2cor(solve(crossprod(m$family$data$R)))
          [,1]      [,2]
[1,] 1.0000000 0.5976769
[2,] 0.5976769 1.0000000

If your response is a count or similar abundance measure, you'll need to transform them if you are going to try the multivariate normal as abundances tend to be skewed and exhibit heterogeneous variance.

Alternatively, if the response matrix $Y$ contains counts from a total, you might use the multinomial response model. This can be done in mgcv also via the multinom family.

There's nothing stopping you from using ns() or bs() terms in the linear model example you had found — you just won't have smoothness selection so you'll need to tune the degrees of freedom/wiggliness of the fitted smooths yourself. bs() and ns() are in the splines package, which comes with R.

Related models could be estimated using a multivariate adaptive regression spline approach, with an R implementation in the earth package, which can now fit these models as GLMs — which might be appropriate given the nature of your response data (counts?)

Alternatively, joint species distribution models (JSDMs) are currently quite popular in ecology and in active development. The idea often involves stacking the columns of your matrix $Y$ and then adding factor variables for site and species which are then introduced into the model via random effects for site:species or latent factor variables.

A recent JSDM overview can be found in Warton et al (2015).

Warton DI, Blanchet FG, O’Hara RB et al. So Many Variables: Joint Modeling in Community Ecology. Trends in ecology & evolution [Internet] 2015;Available from: http://dx.doi.org/10.1016/j.tree.2015.09.007