Using GAM effectively with characters vectors

Question

Perhaps I require a simplifcation of GAM models regarding character vectors, although, I cannot seem to use s() to smooth these vectors, it usually returns an error like:

'names' attribute [2] must be the same length as the vector 1

Essentially, I'm using a GAM because my data has a skewed representation and non-normality. Essentially, my data is of encounter rates that are the probability of encountering any specific bird in a particular region in north America.

The covaraites are landcover; pland_change; ly_name- landcover represent the type of habitat type, pland_change is the difference in pland relative to each year from 2010, and ly_name are specific region names.

For further detail of my methodology, I link here to a question I've asked before that's gone into this with more detail.

My goal is to use the results of the GAM model to produce predictions on a dataset that has the original habitat values (pland), to determine the probability of encounter when pland decreases over the year.

Extra detail for Sir. Gavin Simpson:

I've gone through some articles recently to better understand what model I'm aiming for.

I'm new to GAM, so please go easy on my understanding of the terminology. I wish to develop a model where the predictor (Y) is the encounter probability, which has been developed from randomforests as a likely detection of presence across the landcover variables I'm working with. There are 16 covariate variables, as factors from 1:16; these are across 17 regions, over 2011-2019. pland_change is represented as the difference of the landscape metric from each year since 2010.

My model wishes propose a 'prediction of predicted probabilities', given that the encounter rates are already 'predicted' from presence and absence. My theory is to understand, given these calculations, and the changes in covariates, what's the likelihood of predicting these probabilities over a decreasing landscape in North America, for my bird of interest.

I plan on creating a prediction-model over my original dataframe, with the original values of pland and habitat, to overlay these progressive changes.

I'm planning on using a stepforward selection (Leathwick et al, 2006) to keep only significant covariates, and Becker et al. 2020, mentions using a restricted maximum likelihood to optimize the parameter estimates, which I'm thinking of using alongside stepforward selection to modify the smoothing penalty. Can this be done in gam?

Furthermore, to account for observed geographic differences between regions, and habitat type, pland_change can be modied using a tensor product smooth of latitude and longitude. I plan on calculating the nearest point to each point, non-equivalence so a=b; but b =/ a (=/ means not equal to); to log-transform distance, as an offset for the encounter rate models.

Following Gavins comment, I have performed a gam model as:

m0 <- gam(encounter ~ te(pland_change, landcover, year, ly_name, bs=c('tp', 're','re','re')), data=test.all.2)

After converting landcover and ly_name into factors.

The output is given as this, though I'm hoping for an interpretation of te(pland_change, landcover, year, ly_name, bs=c('tp', 're','re','re') and it's significance to the models output :

Family: gaussian 
Link function: identity 

Formula:
encounter ~ te(pland_change, landcover, year, ly_name, bs = c("tp", 
    "re", "re", "re"))

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.20274    0.04668   4.343 1.41e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                                          edf Ref.df    F p-value    
te(pland_change,landcover,year,ly_name) 636.1    775 21.4  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.311   Deviance explained = 32.4%
GCV = 9.649e-05  Scale est. = 9.4757e-05  n = 35464
> 

References:

1. Leathwick, J.R., Elith, J. and Hastie, T. 2006. Comparative performance of generalized additive models and multivariate adaptive regression splines for statistical modelling of species distributions.
2. Becker et al. 2020. Performance evaluation of cetacean species distribution models developed using generalized additive models and boosted regression trees.

dataset glimpse:

# A tibble: 6 x 7
          X        Y  year landcover pland_change encounter ly_name         
      <dbl>    <dbl> <int> <fct>            <dbl>     <dbl> <fct>           
1 -7615907. 5862248.  2011 barren         -0.143     0.0253 Prairie_Potholes
2 -7215137. 5862248.  2011 barren         -0.143     0.0392 Prairie_Potholes
3 -7215137. 5860168.  2011 barren         -0.143     0.0392 Prairie_Potholes
4 -7390337. 5858088.  2011 barren         -0.0556    0.0246 Prairie_Potholes
5 -7399097. 5851848.  2011 barren         -0.0833    0.0608 Prairie_Potholes
6 -7696937. 5849768.  2011 barren         -0.143     0.0326 Prairie_Potholes

Answer 1

All character variables need to be coded as factors, so coerce them to be factors before fitting the model.

Also, you can't use s() for this as this will make a 3-d thin plate spline smooth and those aren't defined for factor covariates.

What you want instead is a tensor product smooth where the margins can be any smooth basis:

te(x1, x2, x3, bs = c('tp', 're', 're'))

where x2 and x3 are factor covariates for which we use the random effect basis.

You may instead want s(pland_change, by = land over) etc where you have a separate univariate smooth for each land over type, etc.

You'll need to describe more fully what model you want to fit for me to advise more specifically on how to write the formula.