STAT 590 project presentation

1. Marginal association

all$Year=as.factor(all$Year)
par(mfrow=c(1,2))
for(i in c(1,3,5:7,12:13)){
  y1=car::logit(all$carbon/100)
  y2=all$carbon
  x=all[,i]
  # adding r square and p-value
  a = boxplot(y1~x,outline=F,xlab=colnames(all)[i],col=i,ylab = "logit(carbon/100)",ylim=c(min(y1),max(y1)+1))
  text(c(1:nlevels(x)) , a$stats[nrow(a$stats) , ]+1, paste0("p < ", round(summary(lm(y1~x))$coef[,4], 3)))
  b = boxplot(y2~x,outline=F,xlab=colnames(all)[i],ylim=c(min(y2),max(y2)+2),col=i+1,ylab = "carbon")
  text(c(1:nlevels(x)) , b$stats[nrow(b$stats) , ]+1, paste0("p < ", round(summary(lm(y2~x))$coef[,4], 3)))
}

library(ggplot2)
library(ggpubr) 
p1<-ggscatter(all, x = "pH", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
                add = "reg.line", 
                add.params = list(color = "blue", fill = "lightgray"), 
                conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "pH", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,
                    ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "Elevation", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "Elevation", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "ARI_S2", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "ARI_S2", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "DPOL_S1", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "DPOL_S1", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "IRECI_S2", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "IRECI_S2", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "NDVI_S2", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "NDVI_S2", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "NDVIdiff_S2", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "NDVIdiff_S2", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "NDWI1_S2", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "NDWI1_S2", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "NDWIdiff_S2", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "NDWIdiff_S2", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "REIP_S2", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "REIP_S2", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "SWI_LiDAR", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "SWI_LiDAR", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "TPI250_LIDAR", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "TPI250_LIDAR", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "TPI500_LiDAR", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "TPI500_LiDAR", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "TPI750_LiDAR", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "TPI750_LiDAR", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "TRI_LIDAR", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "TRI_LIDAR", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "VBF_LIDAR", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "VBF_LIDAR", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "VH_S1", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "VH_S1", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "VHdiff_S1", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "VHdiff_S1", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "VV_S1", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "VV_S1", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

p1<-ggscatter(all, x = "VVdiff_S1", y = "carbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
p2<-ggscatter(all, x = "VVdiff_S1", y = "logitcarbon",
              palette = "jco",size = 3, alpha = 0.6,color = 'lightblue',
              add = "reg.line", 
              add.params = list(color = "blue", fill = "lightgray"), 
              conf.int = TRUE # Add confidence interval
)+ stat_cor(method = "pearson",label.x.npc = "center")
figure <- ggarrange(p1,p2,ncol = 2, nrow = 1)
figure

Modeling

1. linear regression model: logit(carbon)~.

lm1 = lm(logit(carbon/100)~ Year+nutrient+species+stage+position+direction+Landcover+
           pH+Elevation+ARI_S2+DPOL_S1+IRECI_S2+
           NDVI_S2+NDVIdiff_S2+NDWI1_S2+NDWIdiff_S2+REIP_S2+
      SWI_LiDAR+`TPI250_LIDAR` +`TPI500_LiDAR`+`TPI750_LiDAR`+`TRI_LIDAR`       +`VBF_LIDAR`  +`VH_S1`+`VHdiff_S1` +`VV_S1` +`VVdiff_S1`
      ,data = all1)
summary(lm1)$adj.r.squared

[1] 0.2732792

2. Generalized Additive Model: logit(carbon)~.

gam = gam(logit(carbon)/100~ s(pH)+s(Elevation)+s(ARI_S2)+s(DPOL_S1)+s(IRECI_S2)+s(NDVI_S2)+s(NDVIdiff_S2)+s(NDWI1_S2)+s(NDWIdiff_S2)+s(REIP_S2)+
            s(SWI_LiDAR)+s(TPI250_LIDAR) +s(TPI500_LiDAR)+s(TPI750_LiDAR)+s(TRI_LIDAR)+s(VBF_LIDAR)  +s(VH_S1)+s(VHdiff_S1) +s(VV_S1) +s(VVdiff_S1) +s(VVdiff_S1)+Year+nutrient+species+stage+position+direction+Landcover,data=all1)
summary(gam)$r.sq#

[1] 0.3213097

3. Principal Component Analysis

pca.out = prcomp(all1[,c(11,14:32)],center = T, scale=T )
pve =100* pca.out$sdev ^2/ sum(pca.out$sdev ^2)
par(mfrow =c(1,2))
plot(pve , type ="o", ylab="PVE ", xlab=" Principal Component ",col ="blue")
plot(cumsum(pve), type="o", ylab =" Cumulative PVE", xlab="
Principal Component ", col ="brown3 ")

#pca score
gam1 = gam(logit(carbon)/100 ~ s(PC1)+s(PC2)+s(PC3)+s(PC4)+s(PC5)+s(PC6)+s(PC7)+s(PC8)+s(PC9)+s(PC10)+s(PC11)+s(PC12)+s(PC13)+s(PC14)+s(PC15)+Year+nutrient+species+stage+position+direction+Landcover,data=cbind(pca.out$x,all1))
summary(gam1)$r.sq

[1] 0.3114894

#pcr: numerical variable only
set.seed(2)
pcr.fit=pcr(logitcarbon~.,data=all1[,c(11,14:33)] ,scale=TRUE,center=T)
validationplot(pcr.fit ,val.type="R2")

summary(pcr.fit)

Data:   X dimension: 3350 20 
    Y dimension: 3350 1
Fit method: svdpc
Number of components considered: 20
TRAINING: % variance explained
             1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
X             26.979   41.738   54.921    64.93   72.801   78.357
logitcarbon    0.566    1.355    2.526     2.53    4.124    4.695
             7 comps  8 comps  9 comps  10 comps  11 comps  12 comps
X             83.168   86.936   90.584    93.095    95.027     96.34
logitcarbon    5.778    6.335    8.547     8.758     8.872      8.88
             13 comps  14 comps  15 comps  16 comps  17 comps  18 comps
X              97.420    98.446     99.08     99.46     99.77     99.90
logitcarbon     8.921     9.193      9.30     11.22     11.27     11.75
             19 comps  20 comps
X               99.98    100.00
logitcarbon     11.77     11.94

#pls: numerical variable only
pls.fit=plsr(logitcarbon~., data=all1[,c(11,14:33)] ,scale=T ,center=T)
summary(pls.fit)

Data:   X dimension: 3350 20 
    Y dimension: 3350 1
Fit method: kernelpls
Number of components considered: 20
TRAINING: % variance explained
             1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
X             18.035   35.281   45.566   53.131    62.34    67.91
logitcarbon    5.874    8.224    9.421    9.941    10.42    11.09
             7 comps  8 comps  9 comps  10 comps  11 comps  12 comps
X              71.93    76.03    83.31     87.18     90.33     91.71
logitcarbon    11.42    11.56    11.62     11.67     11.71     11.79
             13 comps  14 comps  15 comps  16 comps  17 comps  18 comps
X               93.47     96.36     97.25     98.27     98.57     99.38
logitcarbon     11.84     11.85     11.86     11.88     11.93     11.93
             19 comps  20 comps
X               99.91    100.00
logitcarbon     11.94     11.94

validationplot(pls.fit ,val.type="R2")

4. Random Forest

library(randomForest)
bag.carbon =randomForest(logitcarbon~.,data=data1,mtry=6,importance=T, proximity=TRUE)
varImpPlot(bag.carbon)

plot(bag.carbon)

print(bag.carbon)

Call:
 randomForest(formula = logitcarbon ~ ., data = data1, mtry = 6,      importance = T, proximity = TRUE) 
               Type of random forest: regression
                     Number of trees: 500
No. of variables tried at each split: 6

          Mean of squared residuals: 0.6548261
                    % Var explained: 38.68

5. SVM

library(e1071) 
svmmodel<-svm(logitcarbon~., data1,cost=100,gamma=0.5)
predsvm<- predict(svmmodel, data1)
R2 <- 1 - (sum((all1$logitcarbon-predsvm)^2)/sum((all1$logitcarbon-mean(all1$logitcarbon))^2))#0.4497074
svmmodel

Call:
svm(formula = logitcarbon ~ ., data = data1, cost = 100, 
    gamma = 0.5)


Parameters:
   SVM-Type:  eps-regression 
 SVM-Kernel:  radial 
       cost:  100 
      gamma:  0.5 
    epsilon:  0.1 


Number of Support Vectors:  2960

R2

[1] 0.9907579