R/05_vital-rates-linearmodel.R

In phytomosaic/pnw: Tree functional strategies in the US Pacific Northwest

######################################################################
#
# GRM -- vital rate regressions
#
#  Rob Smith, phytomosaic@gmail.com, 10 Apr 2020
#
##  CC-BY-SA 4.0 License (Creative Commons Attribution-ShareAlike 4.0)

### TODO: use GAMs not linear models (random effect = plot = unique popn)
### TODO: use mixed models (random effect = plot = unique popn)
### TODO: permit predictor interactions (x1 * x2 * x3 * ...)

### preamble
rm(list=ls())
require(ecole)
require(data.table)
require(mgcv)

### load data
getwd()
load('./data/pnw_tree.rda', verbose=T)
d <- pnw_tree  ;  rm(pnw_tree)
d$nffd <- rowSums(cbind(d$nffd_wt, d$nffd_sp, d$nffd_sm, d$nffd_at))

### CAUTION: tranformations ? ? ? ?
d$dia_y1 <- log10(1 + d$dia_y1)
d$dia_y2 <- log10(1 + d$dia_y2)
set_par(4)
plot(d$dia_y1, d$dial_y1, cex=0.4, pch=16)
plot(d$dia_y1, d$dia_y2, cex=0.4, pch=16)
plot(d$dia_y2, d$dial_y2, cex=0.4, pch=16)
plot(d$dial_y1, d$dial_y2, cex=0.4, pch=16)
abline(0,1, col=2)

### vital rate regressions per species per EACH plot or plot-CLUSTER
`vital_rates` <- function(d_i, covar=c('pc1','pc2'), ...) {
    ### TODO: variable selection (GAM or NPMR)
    ### TODO: flexible number of covariates
    if (length(covar) != 2)
        stop('requires exactly 2 covariates')
    params <- data.frame(surv_int=NA)
    f_s <- as.formula(paste('surv ~ dia_y1 +',
                            paste(covar,collapse='+')))
    f_g <- as.formula(paste('dia_y2 ~ dia_y1 +',
                            paste(covar,collapse='+')))
    data.table::setDF(d_i) # coerce data.table to data.frame
    # survival regression - LOGIT
    reg_surv <- glm(f_s, data=d_i, family=binomial(link='logit'),
                    na.action='na.exclude')
    params$surv_int      <- coefficients(reg_surv)[1]
    params$surv_slope    <- coefficients(reg_surv)[2]
    params$surv_slope_x1 <- coefficients(reg_surv)[3]
    params$surv_slope_x2 <- coefficients(reg_surv)[4]
    # growth regression - LINEAR
    reg_growth             <- lm(f_g, data=d_i, na.action='na.exclude')
    params$growth_int      <- coefficients(reg_growth)[1]
    params$growth_slope    <- coefficients(reg_growth)[2]
    params$growth_slope_x1 <- coefficients(reg_growth)[3]
    params$growth_slope_x2 <- coefficients(reg_growth)[4]
    params$growth_sd       <- sd(resid(reg_growth), na.rm=T)
    # diameter distribution of new recruits
    params$recruitsizemean <- mean(d_i$dia_y2[!is.na(d_i$dia_y1)],
                                   na.rm=T)
    params$recruitsizesd   <-   sd(d_i$dia_y2[!is.na(d_i$dia_y1)],
                                   na.rm=T)
    return(list(params = params,
                models = list(reg_surv   = reg_surv,
                              reg_growth = reg_growth),
                fits = data.frame(tre_cn     = d_i$tre_cn,
                                  fit_surv   = fitted(reg_surv),
                                  fit_growth = fitted(reg_growth))))
}

### vital rates for MULTIPLE species, from most to least frequent
nm <- names(rev(sort(table(d$spp))))[1:45]
nm <- nm[nm!='quercus_douglasii']
vr <- lapply(nm, function(x) {
    cat(x,'\n')
    keep <- d$spp %in% x
    vital_rates(d_i = d[keep,], covar = c('pc1','pc2'))
})
param <- do.call(rbind, lapply(vr, function(x) { x$params }))
dimnames(param)[[1]] <- nm

### compare vital rates among species
`h` <- function(x, ...) {
    hist(x, bre=22, col='grey', xaxs='i', yaxs='i', las=1, main='', ...)
}
set_par(12)
sapply(names(param), function(x) h(param[,x], xlab=x))

### compare fitted values among plots
fits <- lapply(vr, function(x) { data.table::data.table(x$fits) })
fits <- do.call(rbind, fits)
d    <- cbind(d, fits[match(d$tre_cn, fits$tre_cn)])
# mean ACROSS species per plot:
p <- d[, list(plot_surv   = mean(fit_surv, na.rm=T),
              plot_growth = mean(fit_growth, na.rm=T),
              obs_surv    = mean(surv, na.rm=T),
              obs_growth  = mean(dia_y2 - dia_y1, na.rm=T),
              startdiam   = mean(dia_y1, na.rm=T)
), by=plt_cn]
j <- !names(d) %in% c('tre_cn','tre_cn.1','spp','dist','x','y',
                      'is_seed','is_edge','dia_y1','dia_y2',
                      'surv','fit_surv','fit_growth')
p <- cbind(p, d[match(p$plt_cn, d$plt_cn), ..j]) ; rm(j) # plot match

# png('./fig/fig_02_linearmodeldiagnostics.png',9,3,'in',bg='transparent',res=1000)
set_par(6) ; par(mfcol=c(2,3))
plo(p$lon,p$lat,cex=0.5,col=colvec(p$obs_surv),asp=1.6,ylab='',xlab='')
add_text(0.03,0.95,'Obsvd\nsurvival', cex=1.0)
plo(p$lon,p$lat,cex=0.5,col=colvec(p$obs_growth),asp=1.6,ylab='',xlab='')
add_text(0.03,0.95,'Obsvd\ngrowth', cex=1.0)
plo(p$lon,p$lat,cex=0.5,col=colvec(p$plot_surv),asp=1.6,ylab='',xlab='')
add_text(0.03,0.95,'Fitted\nsurvival', cex=1.0)
plo(p$lon,p$lat,cex=0.5,col=colvec(p$plot_growth),asp=1.6,ylab='',xlab='')
add_text(0.03,0.95,'Fitted\ngrowth', cex=1.0)
plo(p$obs_surv, p$plot_surv, cex=0.5, ylab='Fitted survival', xlab='Obsvd survival') ; abline(0,1)
plo(p$obs_growth, p$plot_growth, cex=0.5, ylab='Fitted growth', xlab='Obsvd growth') ; abline(0,1)
# dev.off()

# png('./fig/fig_03_linearmodels-vs-climate.png',9,3,'in',bg='transparent',res=1000)
set_par(12)
plo(p$pc1, p$obs_surv, cex=0.5, xlab='PC1', ylab='Obsvd survival')
plo(p$pc1, p$plot_surv, cex=0.5, xlab='PC1', ylab='Fitted survival')
plo(p$pc1, p$obs_growth, cex=0.5, xlab='PC1', ylab='Obsvd growth')
plo(p$pc1, p$plot_growth, cex=0.5, xlab='PC1', ylab='Fitted growth')
# set_par(4)
plo(p$pc2, p$obs_surv, cex=0.5, xlab='PC2', ylab='Obsvd survival')
plo(p$pc2, p$plot_surv, cex=0.5, xlab='PC2', ylab='Fitted survival')
plo(p$pc2, p$obs_growth, cex=0.5, xlab='PC2', ylab='Obsvd growth')
plo(p$pc2, p$plot_growth, cex=0.5, xlab='PC2', ylab='Fitted growth')
# set_par(4)
plo(p$startdiam, p$obs_surv, cex=0.5, xlab='Diam Y1', ylab='Obsvd survival')
plo(p$startdiam, p$plot_surv, cex=0.5, xlab='Diam Y1', ylab='Fitted survival')
plo(p$startdiam, p$obs_growth, cex=0.5, xlab='Diam Y1', ylab='Obsvd growth')
plo(p$startdiam, p$plot_growth, cex=0.5, xlab='Diam Y1', ylab='Fitted growth')
# dev.off()

### plot regression effect sizes, by species
`plot_se` <- function(xval, yval, se, x_order=TRUE, ylim, ...) {
    o <- if (isTRUE(x_order)) order(yval) else 1:length(yval)
    if (is.character(xval)) {
        nm     <- xval[o]
        yval   <- yval[o]
        se     <- se[o]
        xval   <- 1:length(xval)
    }
    if(missing(ylim)) {
        rng  <- range(c(yval-se, yval+se), na.rm=T)
        buf  <- diff(rng) * 0.1
        ylim <- c(rng[1]-buf, rng[2]+buf)
    }
    par(mar=c(8,2.5,0.5,0.5), las=2, bty='L', yaxs='i', tcl=-0.2,
        mgp=c(1.5,0.4,0), pty='s')
    plot(xval, yval, pch=19, cex.axis=0.7, cex.lab=0.9, xaxt='n',
         ylim=ylim, xlab='', ...)
    axis(1, at=xval, labels=nm, cex.axis=0.7)
    arrows(as.numeric(xval), yval-se, as.numeric(xval), yval+se,
           len=0.02, ang=90, code=3)
}
### recruit sizes
set_par(1)
plot_se(xval = rownames(param), yval = param$recruitsizemean,
        se   = param$recruitsizesd, ylab = 'Recruit size (DBH)')
### growth effects
set_par(4); par(mfrow=c(1,4))
plot_se(xval = rownames(param), yval = param$growth_int,
        se   = param$growth_sd, ylab = 'Growth intercept')
abline(h=0)
plot_se(xval = rownames(param), yval = param$growth_slope,
        se   = param$growth_sd, ylab = 'Effect on growth (Y1 size)')
abline(h=0)
plot_se(xval = rownames(param), yval = param$growth_slope_x1,
        se   = param$growth_sd, ylab = 'Effect on growth (thermal PC)')
abline(h=0)
plot_se(xval = rownames(param), yval = param$growth_slope_x2,
        se   = param$growth_sd, ylab = 'Effect on growth (aridity PC)')
abline(h=0)
### survival effects
set_par(4); par(mfrow=c(1,4))
plot_se(xval = rownames(param), yval = param$surv_int,
        se   = param$growth_sd, ylab = 'Survival intercept')
abline(h=0)
plot_se(xval = rownames(param), yval = param$surv_slope,
        se   = param$growth_sd, ylab = 'Effect on survival (Y1 size)')
abline(h=0)
plot_se(xval = rownames(param), yval = param$surv_slope_x1,
        se   = param$growth_sd, ylab = 'Effect on survival (thermal PC)')
abline(h=0)
plot_se(xval = rownames(param), yval = param$surv_slope_x2,
        se   = param$growth_sd, ylab = 'Effect on survival (aridity PC)')
abline(h=0)

####    END    ####