Cara's RVC GAM snippet

gam_excerpt.R

# for each column in data table
richness.partial.deviance.list = for(j in names(gam_data[,3])) {
  # ??? why are we looping here? 
  # ??? Is this code not the same as setting j="abundance_sum" and not looping? 
  # ??? Maybe there are multiple columns with name "abundance_sum"? Seems odd.
  if(j=="abundance_sum") {  # only this one abundance column 
    # fit a GAM with the following model: 
    # $y_i = \alpha + f_1(X_1i) + f_2(X_2i) + g_1(   X_3i  )       + g_2(X_4i ) + \epsilon_i$
    # $y_i = \alpha + f_1(YEAR) + f_2(PROT) + g_1(LONG, LAT)* YEAR + g_2(STRAT) + \epsilon_i$
    # ??? why are some of these f_* and others g_*?
    # diversity = alpha_weight 
    #           + smoother_function(YEAR) 
    #           + f_2(PROT)  # ??? what is this?
    #           + covariates_strata_smoother(LONG, LAT)*YEAR
    #           + g_2(STRAT)  # ??? what is this?
    #           + residual_error
    #
    # this model is given within the paste( ... ) statement below as:
    #   abundance_sum ~ YEAR + STRATA + PROT
    #                 + mgcv::s(LON_DEGREES, LAT_DEGREES, by=YEAR, bs='ts')
    #                 + mgcv::s(DEPTH, bs='ts')
    fullmod = mgcv::gam(
      formula(
        paste(
          j,
          "~YEAR+STRATA+PROT+mgcv::s(LON_DEGREES,LAT_DEGREES,by=YEAR, bs='ts')+mgcv::s(DEPTH, bs='ts')"
        )
      ), 
      family=poisson,  # ??? I don't see this as a valid input from https://www.rdocumentation.org/packages/mgcv/versions/1.8-36/topics/family.mgcv
      data=gam_data
    )
    # build a null-model gam
    nullmod = mgcv::gam(
      formula(
        paste(j,"~1")
      ), 
      family=poisson, 
      data=gam_data
    )
  }
}

??? why are we looping here?

We are analyzing by year and need to call the column each time we rerun the model.

??? Is this code not the same as setting j="abundance_sum" and not looping?

??? Maybe there are multiple columns with name "abundance_sum"? Seems odd.

To look at different diversity indices within the GAMs, I change the input column and j name. So to look at abundance I would change to that column gam_data[ ,4] and j == "abundance_sum".

??? why are some of these f_* and others g_*?

f_* represent the estimated mean effects for each year, strata, and no-take marine zones, and g_*s are the nonparametric smoothing functions for the continuous covariates depth, and latitude and longitude.

??? I don't see this as a valid input from https://www.rdocumentation.org/packages/mgcv/versions/1.8-36/topics/family.mgcv

This is a fairly common way to use GAMs, I am not sure why it's not in the documentation.

I am going to send you the gam_data input that should work in this model. Run this quick line of code before putting it in the model.

gam_data = read_csv("Data/Diversity_Sample_Data.csv")
gam_data = gam_data %>%
mutate(DEPTH = DEPTH_AVERAGE)%>%
select(-DEPTH_AVERAGE, -evenness, -shannon) %>%
mutate(
STRATA = #group strata
ifelse(grepl("FMLR|FSLR|FDLR", STRAT), 'LINEAR_REEF',
ifelse(grepl("OFPR|MCPR|INPR", STRAT), 'PATCH_REEF',
ifelse(grepl('HRRF', STRAT), 'HIGH_RELIEF', "IDK"))))

ALSO: I have this on github. It says you have an invitation to collaborate, but it might expired. It is RVC_publication.