R语言与统计分析遥感数据科学与R语言

R - 广义加性模型（GAM）构建详解：以归一化植被指数NDVI

2019-12-22  本文已影响0人  TroyShen

目录

0. 问题导入

今天我们通过建立归一化植被指数（NDVI）与气象及生态要素（包括长波辐射，短波辐射，温度，降水，0-200cm土壤水，根系处土壤水，蒸散发以及总初级生产力（GPP））的关系来详细阐述广义加性模型在面向地理及时间序列数据建模方面的应用。

1. 示例数据

本数据为基于GLDAS再分析数据集与NDVI遥感影像在随机点的时间序列，时间长度为2002-4 ~ 2015-12，时间分辨率为月。
点我下载示例数据
示例数据预览：

head(df)
          ndvi       soil1      soil2     soil3      soil4        rain
1  0.237984169 -0.01210715 0.03579731 0.1269299 0.07318894 -0.01543584
2  0.370455335  0.38147139 0.31089661 0.2241396 0.10204067  0.20701857
3  0.331657733  0.41044975 0.48385978 0.4471074 0.25112199  0.62105802
4  0.216662956  0.32583872 0.41198999 0.4231082 0.42613716  0.37216417
5  0.054132382  0.24177292 0.20540345 0.2979310 0.43549429  0.06553887
6 -0.005636952  0.41268755 0.29207486 0.2508858 0.37087816  0.25502620
      longwave    shortwave    root_sm        evap      temper         gpp
1  0.059414987  0.215758745 0.06890271 -0.07747205 0.009909431 -0.04072053
2  0.009142641  0.244385277 0.31129426  0.23793998 0.172678808  0.18329118
3 -0.097150000  0.353491078 0.48706357  0.59985033 0.314583437  0.24478460
4 -0.031527285  0.355970841 0.42948289  0.51469995 0.348457057  0.47172457
5 -0.291598633  0.297255464 0.27643746  0.38420614 0.326270291  0.37802032
6  0.010729154 -0.009709685 0.32028271  0.31684306 0.119881673  0.15986512

1. 数据导入及训练组/验证组拆分（70%/30%）

setwd('L:\\JianShu\\20191222')
df = read.csv('data.csv',header = T)
df = df[,-1]
df = as.data.frame(df)

train = df[1:115,]
test = df[116:165,]

2. 构建并训练广义加性模型（GAM）

注意：本例基于mgcv包的gam函数展开
通过模型训练结果，我们可以判断在众多因子中，在显著性水平小于0.01 的情况下，soil2，gpp，rain与evaporation的时间演变过程对该点归一化植被指数（NDVI）的变化规律会产生显著影响。

library(mgcv)
fit = mgcv::gam(ndvi ~s(soil1)+s(soil2)+s(soil3)+s(soil4)+s(gpp)+
                  s(rain)+s(longwave)+s(shortwave)+s(root_sm)+s(evap)+s(temper),data = train,
          trace = TRUE)
summary(fit)

Family: gaussian 
Link function: identity 

Formula:
ndvi ~ s(soil1) + s(soil2) + s(soil3) + s(soil4) + s(gpp) + s(rain) + 
    s(longwave) + s(shortwave) + s(root_sm) + s(evap) + s(temper)

Parametric coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.011837   0.006971  -1.698   0.0936 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
               edf Ref.df      F  p-value    
s(soil1)     3.916  4.885  2.146  0.07025 .  
s(soil2)     8.391  8.813  3.212  0.00224 ** 
s(soil3)     1.000  1.000  6.147  0.01532 *  
s(soil4)     4.255  5.079  2.601  0.03054 *  
s(gpp)       6.819  7.762  3.107  0.00524 ** 
s(rain)      1.463  1.776 21.553 2.45e-06 ***
s(longwave)  1.000  1.000  0.008  0.92901    
s(shortwave) 1.879  2.356  4.550  0.01028 *  
s(root_sm)   2.679  3.510  3.750  0.01658 *  
s(evap)      1.000  1.000  9.893  0.00234 ** 
s(temper)    5.103  6.322  1.751  0.12357    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.804   Deviance explained = 86.9%
GCV = 0.008401  Scale est. = 0.0055881  n = 115

3. 训练集NDVI~训练结果的相关性与时间演进分析

3.1 训练集NDVI~训练结果的相关性分析（图1）

pl_df_cor = data.frame(SIMU =fit$fitted.values, REAL = train$ndvi)

label_df = data.frame(x = -0.2, y = 0.4, label = paste0('R: ',round(cor(pl_df_cor$SIMU,pl_df_cor$REAL),2)))

p1_cor = ggplot()+
  geom_point(data = pl_df_cor,aes(x= SIMU,y = REAL),color = 'blue',size = 5,
             shape = 1)+
  geom_abline(intercept = 0,slope = 1,size = 1)+
  geom_text(data = label_df,aes(x = x,y = y,label = label),size = 6,color = 'black')+
  theme_bw()+
  theme(
    axis.text =  element_text(face = 'bold',color = 'black',size = 12, hjust = .5),
    axis.title =  element_text(face = 'bold',color = 'black',size = 14, hjust = .5)
  )+
  xlab('SIMULATION RESULT')+
  ylab('REAL NDVI')

png('plot2.png',
    height = 10,
    width = 20,
    units = 'cm',
    res = 800)
print(p1_cor)
dev.off()
图1 训练集NDVI~训练结果的相关性分析图

3.2 训练集NDVI~训练结果的时间演进分析（图2）

date = seq(as.Date('2002-04-01'),
           as.Date('2015-12-01'),
           '1 month')

date_train = date[1:115]



pl_df = data.frame(date = date_train, SIMU = fit$fitted.values, REAL = train$ndvi)
library(reshape2)
library(ggplot2)
pl_df = melt(pl_df,'date')


p1 = ggplot()+
  geom_line(data = pl_df,aes(x = date, y = value,color = variable),size = 1)+
  scale_color_manual(values= c('green','blue'))+
  theme_bw()+
  theme(
    axis.text =  element_text(face = 'bold',color = 'black',size = 12, hjust = .5),
    axis.title =  element_text(face = 'bold',color = 'black',size = 14, hjust = .5),
    legend.position = 'bottom',
    legend.direction = 'horizontal'
  )+
  xlab('Time (month)')+
  ylab('NDVI (1)')

png('plot1.png',
    height = 10,
    width = 20,
    units = 'cm',
    res = 800)
print(p1)
dev.off()
图2 训练集NDVI~训练结果的时间演进分析

4. 验证集NDVI~模型预测结果的相关性与时间演进分析

4.1 模型预测结果计算

predict_test = predict(fit, test)

4.2 验证集NDVI~模型预测结果的相关性分析（图3）
通过对比图1与图3， 我们可以发现模型在验证期与训练期相关性均大于0.7，并未在验证期出现模型预测能力显著下降的问题，证明该模型未产生过拟合。

pl_df_cor = data.frame(SIMU =predict_test, REAL = test$ndvi)

label_df = data.frame(x = -0.2, y = 0.4, label = paste0('R: ',round(cor(pl_df_cor$SIMU,pl_df_cor$REAL),2)))

p2_cor = ggplot()+
  geom_point(data = pl_df_cor,aes(x= SIMU,y = REAL),color = 'blue',size = 5,
             shape = 1)+
  geom_abline(intercept = 0,slope = 1,size = 1)+
  geom_text(data = label_df,aes(x = x,y = y,label = label),size = 6,color = 'black')+
  theme_bw()+
  theme(
    axis.text =  element_text(face = 'bold',color = 'black',size = 12, hjust = .5),
    axis.title =  element_text(face = 'bold',color = 'black',size = 14, hjust = .5)
  )+
  xlab('SIMULATION RESULT')+
  ylab('REAL NDVI')

png('plot3.png',
    height = 20,
    width = 20,
    units = 'cm',
    res = 800)
print(p2_cor)
dev.off()
图3 验证集NDVI~模型预测结果的相关性分析

4.3 验证集NDVI~模型预测结果的时间演进分析（图4）

date_test = date[116:165]
pl_df_test = data.frame(date = date_test, SIMU = predict_test, REAL = test$ndvi)
pl_df_test = melt(pl_df_test,'date')

p2 = ggplot()+
  geom_line(data = pl_df_test,aes(x = date, y = value,color = variable),size = 1)+
  scale_color_manual(values= c('green','blue'))+
  theme_bw()+
  theme(
    axis.text =  element_text(face = 'bold',color = 'black',size = 12, hjust = .5),
    axis.title =  element_text(face = 'bold',color = 'black',size = 14, hjust = .5),
    legend.position = 'bottom',
    legend.direction = 'horizontal'
  )+
  xlab('Time (month)')+
  ylab('NDVI (1)')

png('plot4.png',
    height = 10,
    width = 20,
    units = 'cm',
    res = 800)
print(p2)
dev.off()
图4 验证集NDVI~模型预测结果的时间演进分析

5. 模型训练及验证期间残差分析（图5）

根据图5，我们可以发现模型在训练及验证期间残差均大致服从与均值为1 的正太分布，间接证明模型的实用性。

train_residuals = fit$residuals
test_residuals = test$ndvi - predict_test
residuals1 = data.frame(RESIDUALS = train_residuals,type = 'TRAIN')
residuals2 = data.frame(RESIDUALS = test_residuals,type = 'TEST')
residuals = rbind(residuals1,residuals2)

print(round(mean(residuals1$RESIDUALS),2))
[1] 0
print(round(mean(residuals2$RESIDUALS),2))
[1] 0.01

p3 = ggplot(data = residuals)+
  geom_histogram(aes(x = RESIDUALS, stat(count),fill = type),
           binwidth = 0.05)+
  scale_fill_manual(values = c('green','blue'))+
  theme_bw()+
  theme(
    axis.text =  element_text(face = 'bold',color = 'black',size = 12, hjust = .5),
    axis.title =  element_text(face = 'bold',color = 'black',size = 14, hjust = .5),
    legend.text = element_text(face = 'bold',color = 'black',size = 12, hjust = .5),
    legend.title = element_text(face = 'bold',color = 'black',size = 14, hjust = .5),
    legend.position = 'bottom',
    legend.direction = 'horizontal'
  )+
  xlab('Residuals')+
  ylab('Count')

png('plot4.png',
    height = 20,
    width = 20,
    units = 'cm',
    res = 800)
print(p3)
dev.off()
图5 模型训练及验证期间残差分析

6. 总结

本文主要解决了以下问题：

  1. 如何利用广义加性模型（GAM）面向时间序列建模？
  2. 如何评估多元模型是否发生过拟合及其可用性？

7. 本文所用软件包（没有需要通过install.packages进行安装）

library(reshape2)
library(ggplot2)
library(mgcv)

8. 致谢

首先，感谢大家的持续关注，小编会继续努力，持续更新下去的！

大家如果觉得有用，还麻烦大家转发点赞加关注哈，也可以扩散到朋友圈，多谢大家啦～

大家如果在使用本文代码的过程有遇到问题的，可以留言评论，也可以私信我哈~~


小编联系方式
上一篇下一篇

猜你喜欢

热点阅读