Load package

pacman::p_load(tidyverse,janitor,raster,magrittr,DT)

Load data

setwd("D:/GIS/doc_Elo/R_figure")

#read data and remove sp=Indet and plot 17 18
data <- read.csv("tree_raw_data.csv",stringsAsFactors = F) %>% filter(sp!="Indet." & Plot<=16)

#filter pour les individus sont encore vivent
data %<>% filter(Mort==0)
datatable(data,class="compact")

Calcul AGR - Annual Growth Rate

\(agr =10* \frac{DBH_{max} - DBH_{min}}{Year_{max} - Year_{min} + 1}\)

# calcul des croissances
DBHmin <- tapply(data$DBH, data$idArbre, min)
DBHmax <- tapply(data$DBH, data$idArbre, max)
Y_min <- tapply(data$year,data$idArbre, min)
Y_max <- tapply(data$year,data$idArbre, max)
data2 <- data.frame(idArbre=names(DBHmin), DBHmin, DBHmax, Y_min, Y_max)
id_name <- data %>% dplyr::select(idArbre,Xutm,Yutm,Plot,sp) %>% distinct()
data2 <- left_join(data2,id_name,by="idArbre")

data2$agr<-10*(data2$DBHmax-data2$DBHmin)/(data2$Y_max-data2$Y_min+1)
datatable(data2,class="compact")

Get relative elevation from raster for each individual

#read raster
DEM_relel <- raster("D:/GIS/doc_Elo/matlab/DEM0/DEM2015_Vertical_Distance_to_Channel_Network_2.tif")
#plot(DEM)
data2$relel <- raster::extract(DEM_relel, data2[,c("Xutm","Yutm")])
datatable(data2,class="compact")

Calcul NCI - Neighborhood crowding index (indice de voisinage) (Lasky et al. 2014)

\[NCI_m = \sum_{j=1,m\neq j}^J\frac{DBH_j^2}{d_{mj}^2}\]

#convert DBH unit from cm to m
data2$DBHmin_enm <- data2$DBHmin / 100
data2$DBHmax_enm <- data2$DBHmax / 100


data2$NCI <- vector(nrow(data2),mode="numeric")


# distance entre les individus pour chaque parcelle
for(i in 1 : 16){
  dat <- data2[data2[,"Plot"] == i,]
  dat[,"Xutm"] <- as.numeric(as.character(dat[,"Xutm"]))
  dat[,"Yutm"] <- as.numeric(as.character(dat[,"Yutm"]))
  dist <- as.matrix(dist(dat[,c("Xutm","Yutm")],method="euclidean"))
  row.names(dist) <- dat[,"idArbre"]
  colnames(dist) <- dat[,"idArbre"]
  assign(paste0("dist",i), dist)
  rm(dat)
}



# calcul NCI pour chaque individu

for(i in 1:nrow(data2)){
  #print(i)
  
  d.text=paste0("dist",data2[i,"Plot"])
  d <- eval(parse(text=d.text))
  # if(data[i,"Plot"] == 1) {d = dist1}
  # if(data[i,"Plot"] == 6) {d = dist6}

  
  ind <- data2[i,"idArbre"]
  # selection des ind de la parcelle entre 1m et 20m autour de l'individu concerne cad ind de la ligne i
  t <- d[row.names(d) == ind,]
  t2 <- t[names(t) != ind]
  t3 <- t2[t2 > 1 & t2 < 20.000001]
  if(length(t3)==0) next
  
  dPr <- data2[data2[,"idArbre"] %in% names(t3),] # selection des colonnes pour les ind selectionnes au dessus
  dPr$crow <- vector(nrow(dPr),mode="numeric")
  
  for(n in 1:nrow(dPr)){    
    dPr[n,"crow"] <- (dPr[n,"DBHmin_enm"]^2) / (t3[names(t3) == dPr[n,"idArbre"]]^2) # calcul de dbh? / dist? pour les dPr
  }
  data2[i,"NCI"] <- sum(dPr[,"crow"])
}

We have

datatable(data2,class="compact")

Linear model

Remove outliers in variable “agr”

par(mfrow=c(1,2))
boxplot(data2$agr, main="Original")
data2 <- data2[!data2$agr %in% boxplot.stats(data2$agr)$out,]
boxplot(data2$agr, main="remove outliers")

par(mfrow=c(1,1))

Linear model analysis

Model 1:
\(AGRstand = DBHstand + Residu\)

Model 2:
\(Residu = RelEle + NCI + RelElv * NCI + residu\)

## decoupage du tableau par sp
datsp<-by(data2,data2$sp,subset)

# Standardisation des valeurs d'agr et DBH13 par percentile 95 de chaque espece (on aura valeurs entre 0 et 1)
for(i in 1:length(datsp)){
  datsp[[i]]$agrstand <- datsp[[i]]$agr/(quantile(datsp[[i]]$agr,0.95))
  datsp[[i]]$dbhstand <- datsp[[i]]$DBHmax/(quantile(datsp[[i]]$DBHmax,0.95))
}

# make function
b.get.relel.and.median <- function(df){
  #model 1
  model1 <- lm(agrstand ~ dbhstand, data = df)
  df$residu <- resid(model1)
  #model 2
  model2 <- lm(residu ~ relel+NCI+relel*NCI, data=df)
  #resultat
  res <- c(median(df$relel),model2$coefficients[2])
  names(res) <- c("median_relel","relel")
  return(res)
}

b.get.relel.and.median(datsp[[1]])
## median_relel        relel 
##  1.179626942 -0.003063603
(res<-sapply(datsp,b.get.relel.and.median) %>% t() %>% as.data.frame())
##              median_relel        relel
## obtusa           1.179627 -0.003063603
## sciadophylla     3.127301 -0.001469852

PLOT

#density
g1 <- ggplot(data2, aes(x=relel, color=sp))+
  geom_density(size=1)+theme_bw()+
  scale_color_manual(values = c("blue","red"))+
  guides(color=guide_legend(title="Species"))+
  labs(y="Densité des individus", 
       x="Altitude relative")
#Coefficient
res$sp=row.names(res)
g2 <- ggplot(res, aes(x=median_relel,y=relel, col=sp, shape=sp))+
  geom_point(size=2)+scale_color_manual(values = c("blue","red"))+
  geom_text(data=res,aes(x=median_relel,y=relel,label=sp),vjust = -1)+
  xlim(c(0,15)) +ylim(c(-0.005,0.002))+
  labs(title="Performances", 
       y="Coefficient du paramètre altitude", 
       x="Altitude médiane de la distribution d’abondance (m)")+
  theme_bw()+theme(legend.position="none")

library(ggpubr)
ggarrange(g1,g2,ncol=2)