R/defunct/05_vital-rates_DEFUNCT.R
In phytomosaic/pnw: Tree functional strategies in the US Pacific Northwest
Defines functions
`m`
`get_lambda`
`m`
######################################################################
#
# GRM -- vital rates and IPMs
#
# Rob Smith, phytomosaic@gmail.com, 26 Mar 2020
#
## CC-BY-SA 4.0 License (Creative Commons Attribution-ShareAlike 4.0)
### TODO: IPM use mixed models (random effect = plot = unique popn)
### TODO: permit predictor interactions (x1 * x2 * x3 * ...)
### preamble
rm(list=ls())
require(ecole)
require(viridis)
require(mgcv)
require(rgdal)
require(raster)
cc <- viridis::viridis(99)
# ### load data
getwd()
# d <- read.csv('./data/pnw_tree.csv', stringsAsFactors=F)
load('./data/pnw_tree.rda', verbose=T)
# load('./data/pnw_traits.rda', verbose=T)
d <- pnw_tree
rm(pnw_tree)
# ### CAUTION: tranformations ? ? ? ?
# d$dia_y1 <- log10(1 + d$dia_y1)
# d$dia_y2 <- log10(1 + d$dia_y2)
### explore: plot periodic ann increment for each species
tmp <- d
frq <- aggregate(tmp$spp, list(tmp$spp), length)
(frq <- frq[rev(order(frq$x)),])
nm <- as.character(frq[,1])
tmp$ann_incr <- c(tmp$dia_y2 - tmp$dia_y1) ^ 0.8 # rescale PAI
### get mean per plot
a <- aggregate(list(ann_incr = tmp$ann_incr),
list(plt_cn=tmp$plt_cn, spp=tmp$spp),
mean, na.rm=T, drop=F)
a <- cbind(a, tmp[match(a$plt_cn, tmp$plt_cn),c('lat','lon','elev')])
### black maps
`m` <- function(d, x, cex=1, main='', ...) {
require(mapproj)
maps::map('county', c('wa','or','ca'),
fill = TRUE, col = 'grey80',
ylim=c(38.5, 49.001), xlim=c(-125,-116),
bg='transparent', lforce='e', mar=c(0,0,2,0),
proj='albers', par=c(0,34) ###############
)
points(d$lon, d$lat, pch=16, cex=cex,
col=colvec(d[,x],alpha=0.9, begin=0.05,end=0.99), ...)
title(main=main, line=0.1) #col.main='white')
}
# png('C:/Users/Rob/Desktop/fig_00_topdozen_PAI.png',
# wid=9, hei=5, unit='in', bg='transparent', res=700)
set_par(12)
par(oma=c(0,0,0,0),mfrow=c(2,6))
for (i in 1:12) m(a[a$spp == nm[i],], x='ann_incr',main=nm[i],cex=0.7)
# dev.off()
# rm(tmp, a, m, i, nm, frq)
#################################################################
#### BEGIN: IPM for single species in EACH plot or plot-CLUSTER
### TODO: variable selection (GAM or NPMR)
### pick only one SINGLE SPECIES
d <- d[d$spp %in% c('pseudotsuga_menziesii'),]
### IPM for EACH FIA plot ---> MAT + NFFD
`get_lambda` <- function(d_i, ...) {
# A. vital rate regressions ##############
params <- data.frame(surv_int=NA)
# 1. survival regression
fmla_s <- surv ~ dia_y1 + mat + nffd
surv_reg <- glm(fmla_s, data=d_i, family=binomial(link='logit'))
# summary(surv_reg)
params$surv_int <- coefficients(surv_reg)[1]
params$surv_slope <- coefficients(surv_reg)[2]
params$surv_slope_mat <- coefficients(surv_reg)[3]
params$surv_slope_nffd <- coefficients(surv_reg)[4]
# 2. growth regression - LINEAR
fmla_g <- dia_y2 ~ dia_y1 + mat + nffd
growth_reg <- lm(fmla_g, data=d_i)
# summary(growth_reg)
params$growth_int <- coefficients(growth_reg)[1]
params$growth_slope <- coefficients(growth_reg)[2]
params$growth_slope_mat <- coefficients(growth_reg)[3]
params$growth_slope_nffd <- coefficients(growth_reg)[4]
params$growth_sd <- sd(resid(growth_reg))
# 3. dia_y1 distribution of recruits
params$recruitareamean <- mean(d_i$dia_y2[!is.na(d_i$dia_y1)],
na.rm=T)
params$recruitareasd <- sd(d_i$dia_y2[!is.na(d_i$dia_y1)],
na.rm=T)
# B. build vital rate functions ##############
# 1. probability of surviving
s_x <- function(x, mat, nffd, params) {
u <- exp(params$surv_int +
params$surv_slope * x +
params$surv_slope_mat * mat +
params$surv_slope_nffd * nffd)
return(u / (1 + u))
}
# 2. growth function - LINEAR model
g_yx <- function(y, x, mat, nffd, params) {
dnorm(y,
mean = params$growth_int +
params$growth_slope * x +
params$growth_slope_mat * mat +
params$growth_slope_nffd * nffd,
sd = params$growth_sd)
}
# C. make kernel (for each site) ##############
minsz <- 0.9 * min(c(d_i$dia_y1,d_i$dia_y2),na.rm=T) # lwr
maxsz <- 1.1 * max(c(d_i$dia_y1,d_i$dia_y2),na.rm=T) # upr
n <- 100 # n cells in kernel
b <- minsz+c(0:n)*(maxsz-minsz)/n # boundary points = edges
y <- 0.5*(b[1:n]+b[2:(n+1)]) # mesh points y = midpoints
h <- y[2]-y[1] # step dia_y1 h = width
mat <- mean(d_i$mat, na.rm=T)
nffd <- mean(d_i$nffd, na.rm=T)
G <- h*outer(y,y,g_yx,params=params,mat,nffd) # growth kernel
S <- s_x(y,params=params,mat,nffd) # survival kernel
P <- G # initialize growth/survival
for(i in 1:n) {P[,i] <- G[,i]*S[i]} # growth/survival kernel
K <- P # complete FULL kernel
lam <- Re(eigen(K)$values[1]) # lambda, popn growth rate
return(lam)
# #### COULD ALSO add sum of all tree diameters based on
# #### SSD = stable size distribution
# w_eigen <- Re(eigen(K)$vectors[,1]) # right eigenvector
# stable_dist<- w_eigen/sum(w_eigen) # stable size distribution
# h <- hist(stable_dist, breaks=199, plot=F)
# sum_diam <- sum(h$mids * h$counts, na.rm=T)
# return(sum_diam)
}
### IPM for EACH FIA plot --> MAT + PAS
# `get_lambda` <- function(d_i, ...) {
# # A. vital rate regressions ##############
# params <- data.frame(surv_int=NA)
# # 1. survival regression
# fmla_s <- surv ~ dia_y1 + mat + pas
# surv_reg <- glm(fmla_s, data=d_i, family=binomial(link='logit'))
# # summary(surv_reg)
# params$surv_int <- coefficients(surv_reg)[1]
# params$surv_slope <- coefficients(surv_reg)[2]
# params$surv_slope_mat <- coefficients(surv_reg)[3]
# params$surv_slope_pas <- coefficients(surv_reg)[4]
# # 2. growth regression - LINEAR
# fmla_g <- dia_y2 ~ dia_y1 + mat + pas
# growth_reg <- lm(fmla_g, data=d_i)
# # summary(growth_reg)
# params$growth_int <- coefficients(growth_reg)[1]
# params$growth_slope <- coefficients(growth_reg)[2]
# params$growth_slope_mat <- coefficients(growth_reg)[3]
# params$growth_slope_pas <- coefficients(growth_reg)[4]
# params$growth_sd <- sd(resid(growth_reg))
# # 3. dia_y1 distribution of recruits
# params$recruitareamean <- mean(d_i$dia_y2[!is.na(d_i$dia_y1)],
# na.rm=T)
# params$recruitareasd <- sd(d_i$dia_y2[!is.na(d_i$dia_y1)],
# na.rm=T)
# # B. build vital rate functions ##############
# # 1. probability of surviving
# s_x <- function(x, mat, pas, params) {
# u <- exp(params$surv_int +
# params$surv_slope * x +
# params$surv_slope_mat * mat +
# params$surv_slope_pas * pas)
# return(u / (1 + u))
# }
# # 2. growth function - LINEAR model
# g_yx <- function(y, x, mat, pas, params) {
# dnorm(y,
# mean = params$growth_int +
# params$growth_slope * x +
# params$growth_slope_mat * mat +
# params$growth_slope_pas * pas,
# sd = params$growth_sd)
# }
# # C. make kernel (for each site) ##############
# minsz <- 0.9 * min(c(d_i$dia_y1,d_i$dia_y2),na.rm=T) # lwr
# maxsz <- 1.1 * max(c(d_i$dia_y1,d_i$dia_y2),na.rm=T) # upr
# n <- 100 # n cells in kernel
# b <- minsz+c(0:n)*(maxsz-minsz)/n # boundary points = edges
# y <- 0.5*(b[1:n]+b[2:(n+1)]) # mesh points y = midpoints
# h <- y[2]-y[1] # step dia_y1 h = width
# mat <- mean(d_i$mat, na.rm=T)
# pas <- mean(d_i$pas, na.rm=T)
# G <- h*outer(y,y,g_yx,params=params,mat,pas) # growth kernel
# S <- s_x(y,params=params,mat,pas) # survival kernel
# P <- G # initialize growth/survival
# for(i in 1:n) {P[,i] <- G[,i]*S[i]} # growth/survival kernel
# K <- P # complete FULL kernel
# lam <- Re(eigen(K)$values[1]) # lambda, popn growth rate
# return(lam)
# }
#### end get_lambda ####
#### do it
get_lambda(d) # WORKS across ALL plots; now split by plot-cluster:
set.seed(2027)
### k-means clustering (to increase n trees per 'population')
# must ITERATE to avoid eigen-analysis errors
nn <- c(15,20,25) # varying n plots per cluster
lamn <- data.frame(matrix(NA,nrow=NROW(d),ncol=length(nn),
dimnames=list(NULL,paste0('k',nn))),
stringsAsFactors=F)
for (i in 1:length(nn)) {
n <- nn[i] # n plots per cluster
cat(n,'\n') # n plots per cluster
## k-means clustering (get climate clusters of n plots)
iter <- 0
nmin <- FALSE # reject if < min number of trees per cluster
while(!nmin & iter < 100) {
k <- kmeans(d[,c('mat','pas','lon')],
floor(length(unique(d$plt_cn))/n))
nmin <- all(tapply(k$cluster, k$cluster, length) >=30)
iter <- iter + 1
}
### calculate lambda per cluster
d$k <- k$cluster
sdat <- split(d, d['k']) # split data
ntree<- sapply(sdat, function(x)dim(x)[1])
lams <- sapply(sdat, function(x) get_lambda(x))
lamn[,i] <- lams[match(d$k, names(lams))]
}
# d[,nn] <- NULL
d <- data.frame(d, lamn, stringsAsFactors=F) # combine columns
### some ad-hoc corrections of *extreme* values
d$k15[d$k15>5] <- NA
d$k25[d$k25>4] <- NA
### rescale rgr and annual increment (for plotting)
d$rgr <- log10(1+d$ann_dia_growth /
ifelse(d$dia_begin==0,0.5,d$dia_begin)*100) ^ 1.1
d$ann_incr <- d$ann_dia_growth ^ 1.2
### get mean per plot
a <- aggregate(list(ann_incr = d$ann_incr,
rgr = d$rgr,
k15 = d$k15,
k20 = d$k20,
k25 = d$k25),
list(plt_cn=d$plt_cn),
mean, na.rm=T, drop=F)
a <- cbind(a,d[match(a$plt_cn, d$plt_cn),
c('lat','lon','el','pas','mat','map','cmd')])
### black maps
`m` <- function(d, x, cex=1, ...) {
maps::map('county', 'wa', fill = TRUE, col = 'grey',
ylim=c(45.5, 49.001), xlim=c(-125,-116),
bg='black', lforce='e', mar=c(0,0,0,0))
points(d$lon, d$lat, pch=16, cex=cex,
col=colvec(d[,x],alpha=0.9, begin=0,end=0.99), ...)
}
# png('./fig/fig_01_MAI_RGR_LAMDA.png',
# wid=9, hei=3.0, unit='in', bg='transparent', res=700)
set_par(3)
m(a, x='ann_incr')
m(a, x='rgr')
m(a, x='k25')
# dev.off()
####################### to GAM ####################################
# plot_heatmap(
# cor(d[,c('rgr','ann_incr','mat','mwmt','mcmt','td','map','msp',
# 'ahm','shm','nffd','ffp','pas','emt','ext',
# 'eref','cmd','rh')],
# use='complete.obs', method='kendall')
# )
### lambda vs climate
set_par(4)
plot(d$lon, d$k25)
plot(d$mat, d$k25)
plot(d$nffd, d$k25)
plot(d$el, d$k25)
### single MAI map, for presentation
png('./fig/fig_02_MAI.png',
wid=7.5, hei=4.25, unit='in', bg='transparent', res=700)
set_par(1)
m(a, x='ann_incr', cex=1.5)
dev.off()
### single LAMBDA map, for presentation
png('./fig/fig_03_LAMBDA.png',
wid=7.5, hei=4.25, unit='in', bg='transparent', res=700)
a$k25[a$k25 == max(a$k25)] <- NA
set_par(1)
m(a, x='k25', cex=1.5)
dev.off()
#### END IPMs #############################################
#################################################################
#################################################################
#################################################################
#################################################################
#################################################################
#################################################################
#################################################################
#################################################################
#################################################################
#################################################################
# #### BEGIN PLOTTING #####################################
#
# ### size-distribution of survivor/recruit/mortality classes
# hist(d$ann_dia_growth[is_survivor], breaks=44, col='#00000050',
# main='', xlab='Annual diameter growth')
# hist(d$ann_dia_growth[is_naturaldead], breaks=44, col='#FF000050', add=T)
# hist(d$ann_dia_growth[is_newrecruit], breaks=44, col='#00000050', add=T)
#
# # dia_end from TREE_GRM_COMPONENT vs d$dia from TREE
# # (should match exactly, per documentation)
# plot(d$dia_end, d$dia, cex=0.4)
#
# ### explore abs vs relative growth rates
# set_par(2)
# plot(d$rgr_dia, d$ann_dia_growth, cex=0.3)
# plot(d$rgr_ht, d$ann_ht_growth, cex=0.3)
#
# ### explore RGR vs initial size
# set_par(2)
# plot(d$dia_begin, d$rgr_dia, cex=0.3)
# plot(d$ht_begin, d$rgr_ht, cex=0.3)
#
# ### explore ht vs diameter
# set_par(2)
# plot(d$ann_ht_growth, d$ann_dia_growth, cex=0.3)
# plot(d$rgr_ht, d$rgr_dia, cex=0.3)
#
# ### explore measurement years
# set_par(4)
# hist(d$measyr_2)
# hist(d$measyr_1)
# hist(d$invyr)
# hist(d$remper)
#
# ### WHY isnt diff/yr the same as ann_diam_growth ? ? ? ? ? ?
# plot((d$dia_begin - d$dia_mid) / d$remper, d$ann_dia_growth, cex=0.1)
# abline(0,1)
# ##################################################################
#### END ####
phytomosaic/pnw documentation built on April 16, 2020, 7:29 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions `m` `get_lambda` `m`
######################################################################
#
# GRM -- vital rates and IPMs
#
# Rob Smith, phytomosaic@gmail.com, 26 Mar 2020
#
## CC-BY-SA 4.0 License (Creative Commons Attribution-ShareAlike 4.0)
### TODO: IPM use mixed models (random effect = plot = unique popn)
### TODO: permit predictor interactions (x1 * x2 * x3 * ...)
### preamble
rm(list=ls())
require(ecole)
require(viridis)
require(mgcv)
require(rgdal)
require(raster)
cc <- viridis::viridis(99)
# ### load data
getwd()
# d <- read.csv('./data/pnw_tree.csv', stringsAsFactors=F)
load('./data/pnw_tree.rda', verbose=T)
# load('./data/pnw_traits.rda', verbose=T)
d <- pnw_tree
rm(pnw_tree)
# ### CAUTION: tranformations ? ? ? ?
# d$dia_y1 <- log10(1 + d$dia_y1)
# d$dia_y2 <- log10(1 + d$dia_y2)
### explore: plot periodic ann increment for each species
tmp <- d
frq <- aggregate(tmp$spp, list(tmp$spp), length)
(frq <- frq[rev(order(frq$x)),])
nm <- as.character(frq[,1])
tmp$ann_incr <- c(tmp$dia_y2 - tmp$dia_y1) ^ 0.8 # rescale PAI
### get mean per plot
a <- aggregate(list(ann_incr = tmp$ann_incr),
list(plt_cn=tmp$plt_cn, spp=tmp$spp),
mean, na.rm=T, drop=F)
a <- cbind(a, tmp[match(a$plt_cn, tmp$plt_cn),c('lat','lon','elev')])
### black maps
`m` <- function(d, x, cex=1, main='', ...) {
require(mapproj)
maps::map('county', c('wa','or','ca'),
fill = TRUE, col = 'grey80',
ylim=c(38.5, 49.001), xlim=c(-125,-116),
bg='transparent', lforce='e', mar=c(0,0,2,0),
proj='albers', par=c(0,34) ###############
)
points(d$lon, d$lat, pch=16, cex=cex,
col=colvec(d[,x],alpha=0.9, begin=0.05,end=0.99), ...)
title(main=main, line=0.1) #col.main='white')
}
# png('C:/Users/Rob/Desktop/fig_00_topdozen_PAI.png',
# wid=9, hei=5, unit='in', bg='transparent', res=700)
set_par(12)
par(oma=c(0,0,0,0),mfrow=c(2,6))
for (i in 1:12) m(a[a$spp == nm[i],], x='ann_incr',main=nm[i],cex=0.7)
# dev.off()
# rm(tmp, a, m, i, nm, frq)
#################################################################
#### BEGIN: IPM for single species in EACH plot or plot-CLUSTER
### TODO: variable selection (GAM or NPMR)
### pick only one SINGLE SPECIES
d <- d[d$spp %in% c('pseudotsuga_menziesii'),]
### IPM for EACH FIA plot ---> MAT + NFFD
`get_lambda` <- function(d_i, ...) {
# A. vital rate regressions ##############
params <- data.frame(surv_int=NA)
# 1. survival regression
fmla_s <- surv ~ dia_y1 + mat + nffd
surv_reg <- glm(fmla_s, data=d_i, family=binomial(link='logit'))
# summary(surv_reg)
params$surv_int <- coefficients(surv_reg)[1]
params$surv_slope <- coefficients(surv_reg)[2]
params$surv_slope_mat <- coefficients(surv_reg)[3]
params$surv_slope_nffd <- coefficients(surv_reg)[4]
# 2. growth regression - LINEAR
fmla_g <- dia_y2 ~ dia_y1 + mat + nffd
growth_reg <- lm(fmla_g, data=d_i)
# summary(growth_reg)
params$growth_int <- coefficients(growth_reg)[1]
params$growth_slope <- coefficients(growth_reg)[2]
params$growth_slope_mat <- coefficients(growth_reg)[3]
params$growth_slope_nffd <- coefficients(growth_reg)[4]
params$growth_sd <- sd(resid(growth_reg))
# 3. dia_y1 distribution of recruits
params$recruitareamean <- mean(d_i$dia_y2[!is.na(d_i$dia_y1)],
na.rm=T)
params$recruitareasd <- sd(d_i$dia_y2[!is.na(d_i$dia_y1)],
na.rm=T)
# B. build vital rate functions ##############
# 1. probability of surviving
s_x <- function(x, mat, nffd, params) {
u <- exp(params$surv_int +
params$surv_slope * x +
params$surv_slope_mat * mat +
params$surv_slope_nffd * nffd)
return(u / (1 + u))
}
# 2. growth function - LINEAR model
g_yx <- function(y, x, mat, nffd, params) {
dnorm(y,
mean = params$growth_int +
params$growth_slope * x +
params$growth_slope_mat * mat +
params$growth_slope_nffd * nffd,
sd = params$growth_sd)
}
# C. make kernel (for each site) ##############
minsz <- 0.9 * min(c(d_i$dia_y1,d_i$dia_y2),na.rm=T) # lwr
maxsz <- 1.1 * max(c(d_i$dia_y1,d_i$dia_y2),na.rm=T) # upr
n <- 100 # n cells in kernel
b <- minsz+c(0:n)*(maxsz-minsz)/n # boundary points = edges
y <- 0.5*(b[1:n]+b[2:(n+1)]) # mesh points y = midpoints
h <- y[2]-y[1] # step dia_y1 h = width
mat <- mean(d_i$mat, na.rm=T)
nffd <- mean(d_i$nffd, na.rm=T)
G <- h*outer(y,y,g_yx,params=params,mat,nffd) # growth kernel
S <- s_x(y,params=params,mat,nffd) # survival kernel
P <- G # initialize growth/survival
for(i in 1:n) {P[,i] <- G[,i]*S[i]} # growth/survival kernel
K <- P # complete FULL kernel
lam <- Re(eigen(K)$values[1]) # lambda, popn growth rate
return(lam)
# #### COULD ALSO add sum of all tree diameters based on
# #### SSD = stable size distribution
# w_eigen <- Re(eigen(K)$vectors[,1]) # right eigenvector
# stable_dist<- w_eigen/sum(w_eigen) # stable size distribution
# h <- hist(stable_dist, breaks=199, plot=F)
# sum_diam <- sum(h$mids * h$counts, na.rm=T)
# return(sum_diam)
}
### IPM for EACH FIA plot --> MAT + PAS
# `get_lambda` <- function(d_i, ...) {
# # A. vital rate regressions ##############
# params <- data.frame(surv_int=NA)
# # 1. survival regression
# fmla_s <- surv ~ dia_y1 + mat + pas
# surv_reg <- glm(fmla_s, data=d_i, family=binomial(link='logit'))
# # summary(surv_reg)
# params$surv_int <- coefficients(surv_reg)[1]
# params$surv_slope <- coefficients(surv_reg)[2]
# params$surv_slope_mat <- coefficients(surv_reg)[3]
# params$surv_slope_pas <- coefficients(surv_reg)[4]
# # 2. growth regression - LINEAR
# fmla_g <- dia_y2 ~ dia_y1 + mat + pas
# growth_reg <- lm(fmla_g, data=d_i)
# # summary(growth_reg)
# params$growth_int <- coefficients(growth_reg)[1]
# params$growth_slope <- coefficients(growth_reg)[2]
# params$growth_slope_mat <- coefficients(growth_reg)[3]
# params$growth_slope_pas <- coefficients(growth_reg)[4]
# params$growth_sd <- sd(resid(growth_reg))
# # 3. dia_y1 distribution of recruits
# params$recruitareamean <- mean(d_i$dia_y2[!is.na(d_i$dia_y1)],
# na.rm=T)
# params$recruitareasd <- sd(d_i$dia_y2[!is.na(d_i$dia_y1)],
# na.rm=T)
# # B. build vital rate functions ##############
# # 1. probability of surviving
# s_x <- function(x, mat, pas, params) {
# u <- exp(params$surv_int +
# params$surv_slope * x +
# params$surv_slope_mat * mat +
# params$surv_slope_pas * pas)
# return(u / (1 + u))
# }
# # 2. growth function - LINEAR model
# g_yx <- function(y, x, mat, pas, params) {
# dnorm(y,
# mean = params$growth_int +
# params$growth_slope * x +
# params$growth_slope_mat * mat +
# params$growth_slope_pas * pas,
# sd = params$growth_sd)
# }
# # C. make kernel (for each site) ##############
# minsz <- 0.9 * min(c(d_i$dia_y1,d_i$dia_y2),na.rm=T) # lwr
# maxsz <- 1.1 * max(c(d_i$dia_y1,d_i$dia_y2),na.rm=T) # upr
# n <- 100 # n cells in kernel
# b <- minsz+c(0:n)*(maxsz-minsz)/n # boundary points = edges
# y <- 0.5*(b[1:n]+b[2:(n+1)]) # mesh points y = midpoints
# h <- y[2]-y[1] # step dia_y1 h = width
# mat <- mean(d_i$mat, na.rm=T)
# pas <- mean(d_i$pas, na.rm=T)
# G <- h*outer(y,y,g_yx,params=params,mat,pas) # growth kernel
# S <- s_x(y,params=params,mat,pas) # survival kernel
# P <- G # initialize growth/survival
# for(i in 1:n) {P[,i] <- G[,i]*S[i]} # growth/survival kernel
# K <- P # complete FULL kernel
# lam <- Re(eigen(K)$values[1]) # lambda, popn growth rate
# return(lam)
# }
#### end get_lambda ####
#### do it
get_lambda(d) # WORKS across ALL plots; now split by plot-cluster:
set.seed(2027)
### k-means clustering (to increase n trees per 'population')
# must ITERATE to avoid eigen-analysis errors
nn <- c(15,20,25) # varying n plots per cluster
lamn <- data.frame(matrix(NA,nrow=NROW(d),ncol=length(nn),
dimnames=list(NULL,paste0('k',nn))),
stringsAsFactors=F)
for (i in 1:length(nn)) {
n <- nn[i] # n plots per cluster
cat(n,'\n') # n plots per cluster
## k-means clustering (get climate clusters of n plots)
iter <- 0
nmin <- FALSE # reject if < min number of trees per cluster
while(!nmin & iter < 100) {
k <- kmeans(d[,c('mat','pas','lon')],
floor(length(unique(d$plt_cn))/n))
nmin <- all(tapply(k$cluster, k$cluster, length) >=30)
iter <- iter + 1
}
### calculate lambda per cluster
d$k <- k$cluster
sdat <- split(d, d['k']) # split data
ntree<- sapply(sdat, function(x)dim(x)[1])
lams <- sapply(sdat, function(x) get_lambda(x))
lamn[,i] <- lams[match(d$k, names(lams))]
}
# d[,nn] <- NULL
d <- data.frame(d, lamn, stringsAsFactors=F) # combine columns
### some ad-hoc corrections of *extreme* values
d$k15[d$k15>5] <- NA
d$k25[d$k25>4] <- NA
### rescale rgr and annual increment (for plotting)
d$rgr <- log10(1+d$ann_dia_growth /
ifelse(d$dia_begin==0,0.5,d$dia_begin)*100) ^ 1.1
d$ann_incr <- d$ann_dia_growth ^ 1.2
### get mean per plot
a <- aggregate(list(ann_incr = d$ann_incr,
rgr = d$rgr,
k15 = d$k15,
k20 = d$k20,
k25 = d$k25),
list(plt_cn=d$plt_cn),
mean, na.rm=T, drop=F)
a <- cbind(a,d[match(a$plt_cn, d$plt_cn),
c('lat','lon','el','pas','mat','map','cmd')])
### black maps
`m` <- function(d, x, cex=1, ...) {
maps::map('county', 'wa', fill = TRUE, col = 'grey',
ylim=c(45.5, 49.001), xlim=c(-125,-116),
bg='black', lforce='e', mar=c(0,0,0,0))
points(d$lon, d$lat, pch=16, cex=cex,
col=colvec(d[,x],alpha=0.9, begin=0,end=0.99), ...)
}
# png('./fig/fig_01_MAI_RGR_LAMDA.png',
# wid=9, hei=3.0, unit='in', bg='transparent', res=700)
set_par(3)
m(a, x='ann_incr')
m(a, x='rgr')
m(a, x='k25')
# dev.off()
####################### to GAM ####################################
# plot_heatmap(
# cor(d[,c('rgr','ann_incr','mat','mwmt','mcmt','td','map','msp',
# 'ahm','shm','nffd','ffp','pas','emt','ext',
# 'eref','cmd','rh')],
# use='complete.obs', method='kendall')
# )
### lambda vs climate
set_par(4)
plot(d$lon, d$k25)
plot(d$mat, d$k25)
plot(d$nffd, d$k25)
plot(d$el, d$k25)
### single MAI map, for presentation
png('./fig/fig_02_MAI.png',
wid=7.5, hei=4.25, unit='in', bg='transparent', res=700)
set_par(1)
m(a, x='ann_incr', cex=1.5)
dev.off()
### single LAMBDA map, for presentation
png('./fig/fig_03_LAMBDA.png',
wid=7.5, hei=4.25, unit='in', bg='transparent', res=700)
a$k25[a$k25 == max(a$k25)] <- NA
set_par(1)
m(a, x='k25', cex=1.5)
dev.off()
#### END IPMs #############################################
#################################################################
#################################################################
#################################################################
#################################################################
#################################################################
#################################################################
#################################################################
#################################################################
#################################################################
#################################################################
# #### BEGIN PLOTTING #####################################
#
# ### size-distribution of survivor/recruit/mortality classes
# hist(d$ann_dia_growth[is_survivor], breaks=44, col='#00000050',
# main='', xlab='Annual diameter growth')
# hist(d$ann_dia_growth[is_naturaldead], breaks=44, col='#FF000050', add=T)
# hist(d$ann_dia_growth[is_newrecruit], breaks=44, col='#00000050', add=T)
#
# # dia_end from TREE_GRM_COMPONENT vs d$dia from TREE
# # (should match exactly, per documentation)
# plot(d$dia_end, d$dia, cex=0.4)
#
# ### explore abs vs relative growth rates
# set_par(2)
# plot(d$rgr_dia, d$ann_dia_growth, cex=0.3)
# plot(d$rgr_ht, d$ann_ht_growth, cex=0.3)
#
# ### explore RGR vs initial size
# set_par(2)
# plot(d$dia_begin, d$rgr_dia, cex=0.3)
# plot(d$ht_begin, d$rgr_ht, cex=0.3)
#
# ### explore ht vs diameter
# set_par(2)
# plot(d$ann_ht_growth, d$ann_dia_growth, cex=0.3)
# plot(d$rgr_ht, d$rgr_dia, cex=0.3)
#
# ### explore measurement years
# set_par(4)
# hist(d$measyr_2)
# hist(d$measyr_1)
# hist(d$invyr)
# hist(d$remper)
#
# ### WHY isnt diff/yr the same as ann_diam_growth ? ? ? ? ? ?
# plot((d$dia_begin - d$dia_mid) / d$remper, d$ann_dia_growth, cex=0.1)
# abline(0,1)
# ##################################################################
#### END ####
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.