Métodos e Técnicas de MDS
Predição de atributos do solo usando modelos lineares e não lineares
Elias Mendes Costa
2018-10-15
## Loading required package: lattice## Loading required package: ggplot21 Modelos lineares
1.1 Regressão linear
A equação geral de regressão linear múltipla é dada por: Yi=β0, β1X1, β2X2, β3X3… βiXi + ε onde Yi é o valor da variável dependente na i-ésima observação, β0, β1, β2, β3… βi são parâmetros, X1, X2, X3 … Xi, são constantes conhecidas, ou seja, o valor das variáveis preditoras na i-ésima observação, ε é o erro aleatório com média E{ε} = 0 e variância σ²{εi} = σ², e i = 1, … , n (Hastie et al., 2009).
1.1.1 Dentre as principais vantagens (Potencialidades) tem-se:
São modelos relativamente simples, e sua interpretação é facilitada pois é possível construir equações de predição. Ainda é possível avaliar a contribuição de cada variável para o ajuste do modelo de regressão, bem como selecionar (por exemplo método Stepwise) as variáveis mais relevantes ao modelo. Não exige grande capacidade computacional comparando-se a exigência computacional de modelos não paramétricos, que normalmente são mais complexos. As variáveis preditoras podem ser qualitativas ou variáveis dummy (Chagas et al., 2017; Dotto et al., 2014; Hastie et al., 2009; Moura-Bueno, 2014; Samuel-Rosa, 2012).
1.1.2 Dentre as principais desvantagens (Limitações) tem-se:
As principais limitações da RLM estão relacionadas aos pressupostos estatísticos que são exigidos na modelagem utilizando o modelo de regressão linear. - 1. E ε i = 0 para todo i = 1, 2, … , n. (Os erros possuem média 0). - 2. Os erros são homocedásticos, isto é, Var ε i = σ 2 para todo i = 1, 2, … , n. (Homocedasticidade) - 3. Os erros são independentes, isto é, cov ε i , ε j = 0 para todo i ≠ j. (Independência dos erros) - 4. Os erros têm distribuição normal. (Normalidade) - 5. Exigência de linearidade entre as variáveis regressoras (preditoras) e a variável dependente (propriedade do solo) (Hastie et al., 2009) Ainda estão sujeitos aos efeitos da multicolinearidade. No caso de ocorrência de multicolinearidade das variáveis preditoras há aumento da variância dos parâmetros ajustados (Hair et al., 2009)
2 Modelos não Lineares
2.1 Modelo Aditivo Generalizado (GAM)
A equação geral do GAM é dada por: E(Y|X1,X2, . . . ,Xi)=α + f1(X1) + f2(X2) + · · · + fi(Xi) Como o usual X1,X2, . . . ,Xi representa os preditores (covariáveis) e Y a saída (variável resposta), α é um parâmetro a ser estimado, f são as funções de suavização não especificadas (não paramétricas) (Hastie et al., 2009).
2.1.1 Dentre as principais vantagens (Potencialidades) tem-se:
Possibilita o ajuste dos modelos considerando as relações não lienares entre as variváveis através das funçpões de suavização. Sua interpretação é relativamente simples, comparado a modelos de aprendizado de máquina. Não tem as mesmas exigências dos presupostos estatisticos assim como nos modelos de Regressão Linear. As variáveis preditoras podem ser qualitativas ou variáveis dummy (Hastie et al., 2009; Wood, 2006).
2.1.2 Dentre as principais desvantagens (Limitações) tem-se:
São modelos computacionalmente mais exigentes. Há uma limitação no uso de grande número covariáveis e seleção por métodos como Stepwise quando o número de amostras é pequeno. Pode ser complexa a selação da melhor curva de suavização (Hastie et al., 2009; Wood, 2006).
2.2 Ajustando modelos de Regressão Linear
### carregando os dados
dados <- foreign::read.dbf("../data/Treino1.dbf")
names(dados) # ver os nomes dos campos na tabela## [1] "Perfil" "X" "Y" "Areia" "Silte" "Argila" "CC"
## [8] "PMP" "B1" "B2" "B3" "B4" "B5" "B7"
## [15] "NDVI" "B3_B2" "B3_B7" "B5_B7" "GSI" "ID"str(dados) # ver a estrutura da tabela de dados## 'data.frame': 319 obs. of 20 variables:
## $ Perfil: num 947 1377 404 569 570 ...
## $ X : num 334043 338316 338521 339144 338736 ...
## $ Y : num 8911157 8911314 8911944 8912044 8912440 ...
## $ Areia : num 355 408 295 306 280 ...
## $ Silte : num 235 174 160 174 208 ...
## $ Argila: num 410 418 545 521 512 520 415 216 487 565 ...
## $ CC : num 20.6 20.8 27.1 22.4 25.4 ...
## $ PMP : num 14.8 17.3 20.4 16.2 17.4 ...
## $ B1 : int 99 123 108 111 101 106 118 109 118 104 ...
## $ B2 : int 48 61 51 54 48 51 56 51 59 50 ...
## $ B3 : int 57 73 60 65 55 60 66 59 72 60 ...
## $ B4 : int 65 73 60 64 61 61 65 64 66 60 ...
## $ B5 : int 144 164 140 145 124 139 139 143 141 133 ...
## $ B7 : int 76 88 71 75 61 72 74 71 75 69 ...
## $ NDVI : num 0.06557 0 0 -0.00775 0.05172 ...
## $ B3_B2 : num 1.19 1.2 1.18 1.2 1.15 ...
## $ B3_B7 : num 0.75 0.83 0.845 0.867 0.902 ...
## $ B5_B7 : num 1.89 1.86 1.97 1.93 2.03 ...
## $ GSI : num -0.206 -0.195 -0.219 -0.2 -0.225 ...
## $ ID : int 1 2 4 5 6 8 10 11 12 13 ...
## - attr(*, "data_types")= chr "N" "N" "N" "N" ...summary(dados)## Perfil X Y Areia
## Min. : 9.0 Min. :325003 Min. :8911157 Min. :116.0
## 1st Qu.: 316.5 1st Qu.:332218 1st Qu.:8918697 1st Qu.:211.5
## Median : 599.0 Median :336478 Median :8925102 Median :280.0
## Mean : 667.7 Mean :336446 Mean :8924694 Mean :319.5
## 3rd Qu.: 988.5 3rd Qu.:340396 3rd Qu.:8930513 3rd Qu.:398.2
## Max. :1438.0 Max. :347401 Max. :8938130 Max. :801.5
## Silte Argila CC PMP
## Min. : 41.0 Min. :133.0 Min. : 4.04 Min. : 0.60
## 1st Qu.:189.2 1st Qu.:396.5 1st Qu.:19.30 1st Qu.:13.63
## Median :228.0 Median :486.0 Median :23.10 Median :16.20
## Mean :220.1 Mean :460.4 Mean :22.25 Mean :15.69
## 3rd Qu.:256.0 3rd Qu.:545.0 3rd Qu.:25.94 3rd Qu.:18.60
## Max. :344.0 Max. :647.0 Max. :32.38 Max. :25.19
## B1 B2 B3 B4
## Min. : 80.0 Min. :34.00 Min. : 33.00 Min. :51.00
## 1st Qu.: 97.0 1st Qu.:49.00 1st Qu.: 58.00 1st Qu.:61.00
## Median :107.0 Median :54.00 Median : 65.00 Median :65.00
## Mean :107.4 Mean :53.56 Mean : 64.24 Mean :65.88
## 3rd Qu.:116.0 3rd Qu.:59.00 3rd Qu.: 72.00 3rd Qu.:70.00
## Max. :159.0 Max. :81.00 Max. :100.00 Max. :91.00
## B5 B7 NDVI B3_B2
## Min. : 86.0 Min. : 37.00 Min. :-0.08943 Min. :0.9706
## 1st Qu.:123.0 1st Qu.: 62.00 1st Qu.:-0.04133 1st Qu.:1.1698
## Median :134.0 Median : 70.00 Median : 0.00000 Median :1.1964
## Mean :135.7 Mean : 70.29 Mean : 0.01774 Mean :1.1951
## 3rd Qu.:146.0 3rd Qu.: 77.00 3rd Qu.: 0.05263 3rd Qu.:1.2293
## Max. :191.0 Max. :112.00 Max. : 0.28261 Max. :1.4200
## B3_B7 B5_B7 GSI ID
## Min. :0.7143 Min. :1.705 Min. :-0.32432 Min. : 1.0
## 1st Qu.:0.8530 1st Qu.:1.873 1st Qu.:-0.21576 1st Qu.: 98.0
## Median :0.9041 Median :1.934 Median :-0.19535 Median :202.0
## Mean :0.9182 Mean :1.946 Mean :-0.19413 Mean :199.4
## 3rd Qu.:0.9831 3rd Qu.:2.017 3rd Qu.:-0.17121 3rd Qu.:297.5
## Max. :1.1818 Max. :2.350 Max. :-0.09434 Max. :399.0### Ajustando alguns modelos teste
form0 = stats::formula (Argila~1)
form1 = stats::formula (Argila~B1+B2+B3+B4+B5+B7)
form2 = stats::formula (Argila~X+Y+B1+B2+B3+B4+B5+B7+NDVI+B3_B2+B3_B7+B5_B7)
mod0 <- stats::lm (formula = form0, data=dados)
mod1 <- stats::lm (formula = form1, data=dados)
mod2 <- stats::lm (formula = form2, data=dados)
# Ver os modelos
mod0##
## Call:
## stats::lm(formula = form0, data = dados)
##
## Coefficients:
## (Intercept)
## 460.4mod1##
## Call:
## stats::lm(formula = form1, data = dados)
##
## Coefficients:
## (Intercept) B1 B2 B3 B4
## 725.676 -6.621 19.313 1.979 -6.666
## B5 B7
## 1.790 -7.392mod2##
## Call:
## stats::lm(formula = form2, data = dados)
##
## Coefficients:
## (Intercept) X Y B1 B2
## -2.960e+04 -3.136e-03 3.170e-03 -3.949e+00 6.245e+01
## B3 B4 B5 B7 NDVI
## -5.286e+01 2.858e+01 -1.826e+01 1.105e+01 -4.507e+03
## B3_B2 B3_B7 B5_B7
## 1.858e+03 -1.405e+03 1.169e+03# Test F pela ANOVA
anova(mod0,mod1)## Analysis of Variance Table
##
## Model 1: Argila ~ 1
## Model 2: Argila ~ B1 + B2 + B3 + B4 + B5 + B7
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 318 3727856
## 2 312 2190757 6 1537099 36.485 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mod0,mod2)## Analysis of Variance Table
##
## Model 1: Argila ~ 1
## Model 2: Argila ~ X + Y + B1 + B2 + B3 + B4 + B5 + B7 + NDVI + B3_B2 +
## B3_B7 + B5_B7
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 318 3727856
## 2 306 1718287 12 2009569 29.823 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ANOVA Particionada:
anova(mod0,mod1,mod2)## Analysis of Variance Table
##
## Model 1: Argila ~ 1
## Model 2: Argila ~ B1 + B2 + B3 + B4 + B5 + B7
## Model 3: Argila ~ X + Y + B1 + B2 + B3 + B4 + B5 + B7 + NDVI + B3_B2 +
## B3_B7 + B5_B7
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 318 3727856
## 2 312 2190757 6 1537099 45.622 < 2.2e-16 ***
## 3 306 1718287 6 472470 14.023 4.204e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# test T
summary(mod1)##
## Call:
## stats::lm(formula = form1, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -278.843 -49.413 6.575 60.377 202.345
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 725.676 73.394 9.887 < 2e-16 ***
## B1 -6.621 1.312 -5.048 7.61e-07 ***
## B2 19.313 5.203 3.712 0.000243 ***
## B3 1.979 2.921 0.677 0.498688
## B4 -6.666 1.037 -6.425 4.93e-10 ***
## B5 1.790 1.465 1.222 0.222602
## B7 -7.392 2.276 -3.247 0.001291 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 83.8 on 312 degrees of freedom
## Multiple R-squared: 0.4123, Adjusted R-squared: 0.401
## F-statistic: 36.48 on 6 and 312 DF, p-value: < 2.2e-16summary(mod2)##
## Call:
## stats::lm(formula = form2, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -242.868 -46.856 4.242 50.405 198.722
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.960e+04 7.774e+03 -3.808 0.000169 ***
## X -3.136e-03 1.363e-03 -2.301 0.022062 *
## Y 3.170e-03 8.699e-04 3.643 0.000316 ***
## B1 -3.949e+00 1.343e+00 -2.940 0.003536 **
## B2 6.245e+01 2.771e+01 2.254 0.024915 *
## B3 -5.286e+01 1.684e+01 -3.139 0.001860 **
## B4 2.858e+01 9.445e+00 3.026 0.002686 **
## B5 -1.826e+01 5.766e+00 -3.167 0.001694 **
## B7 1.105e+01 8.244e+00 1.340 0.181253
## NDVI -4.507e+03 1.197e+03 -3.765 0.000199 ***
## B3_B2 1.858e+03 1.166e+03 1.593 0.112121
## B3_B7 -1.405e+03 5.046e+02 -2.785 0.005692 **
## B5_B7 1.169e+03 3.957e+02 2.955 0.003373 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 74.94 on 306 degrees of freedom
## Multiple R-squared: 0.5391, Adjusted R-squared: 0.521
## F-statistic: 29.82 on 12 and 306 DF, p-value: < 2.2e-16No primeiro momento, pelo test F, foi testada se é significativa a inclusão das variáveis, ideal que fosse testado cada variável separadamente, mas não se aplica por exemplo em casos que há centenas e até mesmo milhares de variáveis regressoras (covariáveis). Após ter testado a significância dos modelos em relação ao molode apenas com o intercepto, foi feito o teste para verificar se a contribuição de uma determinada covariável é significativa pelo test T.
## Warning in lm(formula = form2, method = "forward", data = dados): method =
## 'forward' is not supported. Using 'qr'## Start: AIC=2766.74
## Argila ~ X + Y + B1 + B2 + B3 + B4 + B5 + B7 + NDVI + B3_B2 +
## B3_B7 + B5_B7
##
## Df Sum of Sq RSS AIC
## - B7 1 10082 1728369 2766.6
## <none> 1718287 2766.7
## - B3_B2 1 14255 1732542 2767.4
## - B2 1 28524 1746811 2770.0
## - X 1 29732 1748019 2770.2
## - B3_B7 1 43543 1761830 2772.7
## - B1 1 48524 1766811 2773.6
## - B5_B7 1 49023 1767310 2773.7
## - B4 1 51428 1769715 2774.1
## - B3 1 55334 1773621 2774.8
## - B5 1 56329 1774616 2775.0
## - Y 1 74543 1792830 2778.3
## - NDVI 1 79611 1797898 2779.2## Warning in lm(formula = Argila ~ X + Y + B1 + B2 + B3 + B4 + B5 + NDVI + :
## method = 'forward' is not supported. Using 'qr'##
## Step: AIC=2766.6
## Argila ~ X + Y + B1 + B2 + B3 + B4 + B5 + NDVI + B3_B2 + B3_B7 +
## B5_B7
##
## Df Sum of Sq RSS AIC
## - B3_B2 1 9262 1737631 2766.3
## <none> 1728369 2766.6
## - B2 1 21407 1749776 2768.5
## - X 1 28532 1756901 2769.8
## - B3_B7 1 39987 1768356 2771.9
## - B1 1 42800 1771169 2772.4
## - B4 1 42896 1771265 2772.4
## - B3 1 45702 1774071 2772.9
## - B5_B7 1 49629 1777998 2773.6
## - NDVI 1 69790 1798160 2777.2
## - B5 1 70314 1798683 2777.3
## - Y 1 77829 1806198 2778.7## Warning in lm(formula = Argila ~ X + Y + B1 + B2 + B3 + B4 + B5 + NDVI + :
## method = 'forward' is not supported. Using 'qr'##
## Step: AIC=2766.31
## Argila ~ X + Y + B1 + B2 + B3 + B4 + B5 + NDVI + B3_B7 + B5_B7
##
## Df Sum of Sq RSS AIC
## <none> 1737631 2766.3
## - B3_B7 1 32716 1770348 2770.3
## - X 1 36533 1774165 2770.9
## - B1 1 37134 1774765 2771.1
## - B5_B7 1 41351 1778982 2771.8
## - B5 1 62059 1799691 2775.5
## - Y 1 77611 1815242 2778.2
## - B2 1 82292 1819923 2779.1
## - B3 1 97868 1835500 2781.8
## - B4 1 141489 1879121 2789.3
## - NDVI 1 211605 1949237 2801.0## [1] "NDVI" "B3_B2" "B3_B7" "B5_B7" "B2" "B3" "B4" "B5"
## [9] "B7" "B1" "Y" "X"##
## Call:
## lm(formula = y ~ ., data = tmp)
##
## Coefficients:
## (Intercept) NDVI B3_B2 B3_B7 B5_B7
## 4.346e+04 5.467e+03 -2.437e+03 2.055e+03 -1.353e+03
## B2 B3 B4 B5 B7
## -8.102e+01 6.271e+01 -3.631e+01 2.243e+01 -8.192e+00
## B1 Y X
## 5.917e+00 -4.535e-03 1.473e-03### Algumas informações diagnosticas sobre os modelos ajustados (MLR)
stats::shapiro.test(resid(mod1)) # Testar os resíduos de normalidade do modelo de regressão##
## Shapiro-Wilk normality test
##
## data: resid(mod1)
## W = 0.98555, p-value = 0.002755lmtest::bptest(mod1) # Homoscedasticidade de variância##
## studentized Breusch-Pagan test
##
## data: mod1
## BP = 25.824, df = 6, p-value = 0.0002401par(mfrow=c(2,2))
plot(mod1,which=1:4)
stats::shapiro.test(resid(mod3)) # Testar os resíduos de normalidade do modelo de regressão##
## Shapiro-Wilk normality test
##
## data: resid(mod3)
## W = 0.98947, p-value = 0.02119lmtest::bptest(mod3) # Homoscedasticidade de variância##
## studentized Breusch-Pagan test
##
## data: mod3
## BP = 27.272, df = 10, p-value = 0.002358par(mfrow=c(2,2))
plot(mod3,which=1:4)
# Escrevendo os resultados de ajuste dos modelos numa tabela
resultados = as.data.frame(matrix(NA,4,4)) # 4 linha e 4 colunas
names(resultados) = c("Modelo","R2", "n_covars","covars")
#mod1
resultados[1,1] = paste(as.character(formula(mod1)[2]),"mod1", sep=" ")
resultados[1,2] = summary(mod1)$r.squared
resultados[1,3] = nrow(summary(mod1)$coefficients)-1
resultados[1,4] = as.character(formula(mod1)[3])
#mod2
resultados[2,1] = paste(as.character(formula(mod2)[2]),"mod2", sep=" ")
resultados[2,2] = summary(mod2)$r.squared
resultados[2,3] = nrow(summary(mod2)$coefficients)-1
resultados[2,4] = as.character(formula(mod2)[3])
#mod3
resultados[3,1] = paste(as.character(formula(mod3)[2]),"mod3", sep=" ")
resultados[3,2] = summary(mod3)$r.squared
resultados[3,3] = nrow(summary(mod3)$coefficients)-1
resultados[3,4] = as.character(formula(mod3)[3])
#mod4
resultados[4,1] = paste(as.character(formula(mod3)[2]),"mod4", sep=" ")
resultados[4,2] = mod4$results$Rsquared[9]
resultados[4,3] = mod4$bestSubset
resultados[4,4] = as.character(formula(mod2)[3])
resultados## Modelo R2 n_covars
## 1 Argila mod1 0.4123280 6
## 2 Argila mod2 0.5390683 12
## 3 Argila mod3 0.5338791 10
## 4 Argila mod4 0.4976679 12
## covars
## 1 B1 + B2 + B3 + B4 + B5 + B7
## 2 X + Y + B1 + B2 + B3 + B4 + B5 + B7 + NDVI + B3_B2 + B3_B7 + B5_B7
## 3 X + Y + B1 + B2 + B3 + B4 + B5 + NDVI + B3_B7 + B5_B7
## 4 X + Y + B1 + B2 + B3 + B4 + B5 + B7 + NDVI + B3_B2 + B3_B7 + B5_B7write.csv(resultados, file="../res/tab/ajuste_modelos.MLR.csv")### Vazendo a validação dos resultados usando um conjunto de amostra independente
# Lwndo os dados de validação
dados.valida <- read.dbf("../data/Teste1.dbf")
## Fazendo as predições nos dados de validação usando os 4 modelos
valida1 <- predict(mod1, newdata = dados.valida)
valida2 <- predict(mod2, newdata = dados.valida)
valida3 <- predict(mod3, newdata = dados.valida)
valida4 <- predict(mod4, newdata = dados.valida)
dados.valida=cbind(dados.valida, valida1, valida2, valida3, valida4) # juntar os dataframes das prediçpões com os valores obervados
## Calculando o R2 da validação
R.mod1=summary(lm(Argila~valida1,data=dados.valida))
R.mod2=summary(lm(Argila~valida2,data=dados.valida))
R.mod3=summary(lm(Argila~valida3,data=dados.valida))
R.mod4=summary(lm(Argila~valida4,data=dados.valida))
## Calculando RMSE e MSE
RMSE1 <- rmse(dados.valida$Argila, dados.valida$valida1); MSE1 <- mse(dados.valida$Argila, dados.valida$valida1)
RMSE2 <- rmse(dados.valida$Argila, dados.valida$valida2); MSE2 <- mse(dados.valida$Argila, dados.valida$valida2)
RMSE3 <- rmse(dados.valida$Argila, dados.valida$valida3); MSE3 <- mse(dados.valida$Argila, dados.valida$valida3)
RMSE4 <- rmse(dados.valida$Argila, dados.valida$valida4); MSE4 <- mse(dados.valida$Argila, dados.valida$valida4) ### Criando uma tabela de resultados com o resultado da validação externa
## Pegando os resultados de R2, RMSE e MSE de cada modelo
rlm.df.rsq=c(R.mod1$r.squared,R.mod2$r.squared,R.mod3$r.squared,R.mod4$r.squared)
rlm.df.RMSE=c(RMSE1, RMSE2, RMSE3, RMSE4)
rlm.df.MSE=c(MSE1, MSE2, MSE3, MSE4)
mods=c("RLM1","RLM2","RLM3", "RLM4")
## escolhendo o número de dígitos
rlm.df.rsq=round(rlm.df.rsq, digits = 2)
rlm.df.RMSE=round(rlm.df.RMSE, digits = 3)
rlm.df.MSE=round(rlm.df.MSE, digits = 3)
## criando um dataframe com os resultados
rlm.df.external=cbind(Modelo=mods,R2=rlm.df.rsq, RMSE=rlm.df.RMSE, MSE=rlm.df.MSE)
rlm.df.external=as.data.frame(rlm.df.external)
rlm.df.external## Modelo R2 RMSE MSE
## 1 RLM1 0.43 80.923 6548.593
## 2 RLM2 0.53 73.659 5425.719
## 3 RLM3 0.54 72.385 5239.58
## 4 RLM4 0.51 234.704 55085.835write.csv(rlm.df.external, file = "../res/tab/RLM.df.external.csv")## png
## 2## png
## 2
Figure 2.1: validação
Figure 2.2: validação
### Gerando o mapa do teor de argila usano o nosso melhor modelo
# importar o grid da area de estudo
gride <- read.dbf("../data/Pontos1.dbf")
gride$X=gride$POINT_X; gride$Y=gride$POINT_Y
names (gride) # ver os nomes das colunas na tabela## [1] "POINTID" "POINT_X" "POINT_Y" "B1" "B2" "B3" "B4"
## [8] "B5" "B7" "B3_B2" "B3_B7" "B5_B7" "NDVI" "X"
## [15] "Y"pred<-predict(mod3, gride, interval="conf", level=0.95)
pred[1:10,] # mostrar as 10 primeiras linhas da predicao## fit lwr upr
## 1 445.8679 416.2506 475.4853
## 2 399.9032 364.0972 435.7091
## 3 460.8458 417.7639 503.9277
## 4 474.6659 439.9820 509.3498
## 5 414.5493 383.2960 445.8027
## 6 322.2474 277.9847 366.5100
## 7 385.7663 344.1542 427.3783
## 8 394.2849 355.5932 432.9766
## 9 402.4877 366.0066 438.9688
## 10 410.0618 374.8560 445.2676predicao=cbind(pred[,1],gride[,2:3]) # este data.frame pode ser salvo como CSV e aberto no ArcGIS
predicao[1:10,] # ver os nomes das colunas na tabela## pred[, 1] POINT_X POINT_Y
## 1 445.8679 337807.4 8938236
## 2 399.9032 337537.4 8938206
## 3 460.8458 337567.4 8938206
## 4 474.6659 337597.4 8938206
## 5 414.5493 337627.4 8938206
## 6 322.2474 337657.4 8938206
## 7 385.7663 337687.4 8938206
## 8 394.2849 337717.4 8938206
## 9 402.4877 337747.4 8938206
## 10 410.0618 337777.4 8938206summary(predicao$fit)## Length Class Mode
## 0 NULL NULLgridded(predicao) <- ~POINT_X+POINT_Y # transformar os dados de dataframe para spatialpixeldataframe
plot(predicao)
#### PACOTE ESPECIAL PARA EXPORTAR ####
library(rgdal)## rgdal: version: 1.3-4, (SVN revision 766)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: /usr/share/gdal/2.2
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.3-1writeGDAL(predicao[,1], "../res/fig/Argila.RLM.sdat", drivername="SAGA",mvFlag=-99999)2.3 Ajustando modelos GAM
GAM1 = mgcv::gam (Argila~s(X,Y, bs="ts")+s(B1, bs="ts")+s(B2, bs="ts")+
s(B3, bs="ts")+s(B4, bs="ts")+ s(B5, bs="ts")+
s(B7, bs="ts"), data=dados)
summary(GAM1)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Argila ~ s(X, Y, bs = "ts") + s(B1, bs = "ts") + s(B2, bs = "ts") +
## s(B3, bs = "ts") + s(B4, bs = "ts") + s(B5, bs = "ts") +
## s(B7, bs = "ts")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 460.408 3.448 133.5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X,Y) 1.955e+01 29 3.594 3.51e-14 ***
## s(B1) 5.362e-01 9 0.166 0.07881 .
## s(B2) 8.772e+00 9 8.784 4.97e-14 ***
## s(B3) 2.425e-09 9 0.000 0.60469
## s(B4) 1.105e+00 9 2.089 5.09e-06 ***
## s(B5) 7.086e+00 9 2.051 0.00151 **
## s(B7) 3.115e+00 9 0.757 0.04427 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.677 Deviance explained = 71.7%
## GCV = 4354 Scale est. = 3792.1 n = 319GAM2 = mgcv::gam (Argila~s(X,Y, bs="ts")+s(B1, bs="ts")+s(B2, bs="ts")+
s(B3, bs="ts")+s(B4, bs="ts")+ s(B5, bs="ts")+
s(B7, bs="ts")+s(NDVI, bs="ts")+s(B3_B2, bs="ts")+
s(B3_B7, bs="ts")+s(B5_B7, bs="ts"),
method="REML", select=TRUE, data=dados)
summary(GAM2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Argila ~ s(X, Y, bs = "ts") + s(B1, bs = "ts") + s(B2, bs = "ts") +
## s(B3, bs = "ts") + s(B4, bs = "ts") + s(B5, bs = "ts") +
## s(B7, bs = "ts") + s(NDVI, bs = "ts") + s(B3_B2, bs = "ts") +
## s(B3_B7, bs = "ts") + s(B5_B7, bs = "ts")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 460.408 3.594 128.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X,Y) 1.767e+01 29 3.440 1.09e-14 ***
## s(B1) 1.084e+00 9 2.359 8.74e-07 ***
## s(B2) 3.835e-04 9 0.000 0.458
## s(B3) 1.195e-04 9 0.000 0.628
## s(B4) 1.559e-04 9 0.000 0.546
## s(B5) 3.150e-04 9 0.000 0.466
## s(B7) 6.358e-05 9 0.000 0.853
## s(NDVI) 3.284e+00 9 7.033 < 2e-16 ***
## s(B3_B2) 1.030e+00 9 1.664 3.15e-05 ***
## s(B3_B7) 2.434e+00 9 1.935 2.72e-05 ***
## s(B5_B7) 7.456e-02 9 0.009 0.281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.649 Deviance explained = 67.7%
## -REML = 1821.9 Scale est. = 4119.9 n = 319GAM3 <- mgcv::gam (Argila~s(X,Y, bs="ts")+s(B1, bs="ts")+s(B2, bs="ts")+
s(B3, bs="ts")+s(B4, bs="ts")+ s(B5, bs="ts")+
s(B7, bs="ts")+s(NDVI, bs="ts")+s(B3_B2, bs="ts")+
s(B3_B7, bs="ts")+s(B5_B7, bs="ts"),
select=TRUE, data=dados)
summary(GAM3)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Argila ~ s(X, Y, bs = "ts") + s(B1, bs = "ts") + s(B2, bs = "ts") +
## s(B3, bs = "ts") + s(B4, bs = "ts") + s(B5, bs = "ts") +
## s(B7, bs = "ts") + s(NDVI, bs = "ts") + s(B3_B2, bs = "ts") +
## s(B3_B7, bs = "ts") + s(B5_B7, bs = "ts")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 460.408 3.227 142.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X,Y) 1.970e+01 29 3.203 3.63e-12 ***
## s(B1) 8.884e+00 9 3.557 0.000130 ***
## s(B2) 8.931e+00 9 3.788 1.74e-05 ***
## s(B3) 9.590e-10 9 0.000 0.728406
## s(B4) 1.063e-08 9 0.000 0.611516
## s(B5) 8.540e+00 9 3.852 2.18e-06 ***
## s(B7) 1.381e+00 9 0.534 0.004967 **
## s(NDVI) 4.362e+00 9 4.589 4.27e-10 ***
## s(B3_B2) 3.328e+00 9 0.373 0.206321
## s(B3_B7) 2.026e+00 9 0.783 0.005827 **
## s(B5_B7) 9.246e-01 9 0.333 0.000612 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.717 Deviance explained = 76.8%
## GCV = 4076.6 Scale est. = 3321.7 n = 319## Esse demora bastante
gam <- caret::rfeControl (functions=gamFuncs, method= "LOOCV") # gam function
GAM4 <- caret::rfe(dados[,c(2,3,9:18)], dados[,4], sizes=c(1:8), rfeControl=gam)### Algumas informações diagnosticas sobre os modelos ajustados (GAM)
#GAM1
par(mfrow=c(2,2))
mgcv::gam.check(GAM1,pch=19,cex=.3, old.style=T)
##
## Method: GCV Optimizer: magic
## Smoothing parameter selection converged after 86 iterations.
## The RMS GCV score gradient at convergence was 0.0001789549 .
## The Hessian was positive definite.
## Model rank = 84 / 84
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(X,Y) 2.90e+01 1.96e+01 0.99 0.36
## s(B1) 9.00e+00 5.36e-01 0.98 0.34
## s(B2) 9.00e+00 8.77e+00 1.00 0.54
## s(B3) 9.00e+00 2.42e-09 1.02 0.58
## s(B4) 9.00e+00 1.11e+00 1.07 0.89
## s(B5) 9.00e+00 7.09e+00 0.97 0.27
## s(B7) 9.00e+00 3.11e+00 1.02 0.65#GAM3
par(mfrow=c(2,2))
mgcv::gam.check(GAM3,pch=19,cex=.3, old.style=T)
##
## Method: GCV Optimizer: magic
## Smoothing parameter selection converged after 76 iterations.
## The RMS GCV score gradient at convergence was 8.869209e-05 .
## The Hessian was positive definite.
## Model rank = 120 / 120
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(X,Y) 2.90e+01 1.97e+01 1.00 0.47
## s(B1) 9.00e+00 8.88e+00 1.01 0.55
## s(B2) 9.00e+00 8.93e+00 1.02 0.66
## s(B3) 9.00e+00 9.59e-10 1.04 0.77
## s(B4) 9.00e+00 1.06e-08 1.05 0.80
## s(B5) 9.00e+00 8.54e+00 0.96 0.24
## s(B7) 9.00e+00 1.38e+00 1.04 0.68
## s(NDVI) 9.00e+00 4.36e+00 1.01 0.49
## s(B3_B2) 9.00e+00 3.33e+00 1.03 0.70
## s(B3_B7) 9.00e+00 2.03e+00 1.06 0.83
## s(B5_B7) 9.00e+00 9.25e-01 1.08 0.94## Escrevendo os resultados de ajuste dos modelos numa tabela
resultados = as.data.frame(matrix(NA,4,4)) # 4 linha e 4 colunas
names(resultados) = c("Modelo","R2", "n_covars","covars")
#GAM1
resultados[1,1] = paste(as.character(formula(mod1)[2]),"GAM1", sep=" ")
resultados[1,2] = summary(GAM1)$r.sq
resultados[1,3] = summary(GAM1)$m+1
resultados[1,4] = as.character(formula(GAM1)[3])
#GAM2
resultados[2,1] = paste(as.character(formula(mod2)[2]),"GAM2", sep=" ")
resultados[2,2] = summary(GAM2)$r.sq
resultados[2,3] = summary(GAM2)$m+1
resultados[2,4] = as.character(formula(GAM2)[3])
#GAM3
resultados[3,1] = paste(as.character(formula(mod3)[2]),"GAM3", sep=" ")
resultados[3,2] = summary(GAM3)$r.sq
resultados[3,3] = summary(GAM3)$m+1
resultados[3,4] = as.character(formula(GAM3)[3])
#GAM4
resultados[4,1] = paste(as.character(formula(mod3)[2]),"GAM4", sep=" ")
resultados[4,2] = GAM4$results$Rsquared[4]
resultados[4,3] = GAM4$bestSubset
resultados[4,4] = as.character(GAM4$fit$pred.formula[2])
resultados## Modelo R2 n_covars
## 1 Argila GAM1 0.6765161 8
## 2 Argila GAM2 0.6485601 12
## 3 Argila GAM3 0.7166480 12
## 4 Argila GAM4 0.5181448 3
## covars
## 1 s(X, Y, bs = "ts") + s(B1, bs = "ts") + s(B2, bs = "ts") + s(B3, bs = "ts") + s(B4, bs = "ts") + s(B5, bs = "ts") + s(B7, bs = "ts")
## 2 s(X, Y, bs = "ts") + s(B1, bs = "ts") + s(B2, bs = "ts") + s(B3, bs = "ts") + s(B4, bs = "ts") + s(B5, bs = "ts") + s(B7, bs = "ts") + s(NDVI, bs = "ts") + s(B3_B2, bs = "ts") + s(B3_B7, bs = "ts") + s(B5_B7, bs = "ts")
## 3 s(X, Y, bs = "ts") + s(B1, bs = "ts") + s(B2, bs = "ts") + s(B3, bs = "ts") + s(B4, bs = "ts") + s(B5, bs = "ts") + s(B7, bs = "ts") + s(NDVI, bs = "ts") + s(B3_B2, bs = "ts") + s(B3_B7, bs = "ts") + s(B5_B7, bs = "ts")
## 4 Y + B2 + B3_B7write.csv(resultados, file="../res/tab/ajuste_modelos.GAM.csv")### Vazendo a validação dos resultados usando um conjunto de amostra independente
# Lwndo os dados de validação
dados.valida <- read.dbf("../data/Teste1.dbf")
## Fazendo as predições nos dados de validação usando os 4 modelos
valida1 <- predict(GAM1, newdata = dados.valida)
valida2 <- predict(GAM2, newdata = dados.valida)
valida3 <- predict(GAM3, newdata = dados.valida)
valida4 <- predict(GAM4, newdata = dados.valida)
dados.valida=cbind(dados.valida, valida1, valida2, valida3, valida4);names(dados.valida)[24]<-c("valida4") # juntar os dataframes das prediçpões com os valores obervados
## Calculando o R2 da validação
R.mod1=summary(lm(Argila~valida1,data=dados.valida))
R.mod2=summary(lm(Argila~valida2,data=dados.valida))
R.mod3=summary(lm(Argila~valida3,data=dados.valida))
R.mod4=summary(lm(Argila~valida4,data=dados.valida))
## Calculando RMSE e MSE
RMSE1 <- rmse(dados.valida$Argila, dados.valida$valida1); MSE1 <- mse(dados.valida$Argila, dados.valida$valida1)
RMSE2 <- rmse(dados.valida$Argila, dados.valida$valida2); MSE2 <- mse(dados.valida$Argila, dados.valida$valida2)
RMSE3 <- rmse(dados.valida$Argila, dados.valida$valida3); MSE3 <- mse(dados.valida$Argila, dados.valida$valida3)
RMSE4 <- rmse(dados.valida$Argila, dados.valida$valida4); MSE4 <- mse(dados.valida$Argila, dados.valida$valida4) ### Criando uma tabela de resultados com o resultado da validação externa
## Pegando os resultados de R2, RMSE e MSE de cada modelo
gam.df.rsq=c(R.mod1$r.squared,R.mod2$r.squared,R.mod3$r.squared,R.mod4$r.squared)
gam.df.RMSE=c(RMSE1, RMSE2, RMSE3, RMSE4)
gam.df.MSE=c(MSE1, MSE2, MSE3, MSE4)
mods=c("GAM1","GAM2","GAM3", "GAM4")
## escolhendo o número de dígitos
gam.df.rsq=round(gam.df.rsq, digits = 2)
gam.df.RMSE=round(gam.df.RMSE, digits = 3)
gam.df.MSE=round(gam.df.MSE, digits = 3)
## criando um dataframe com os resultados
gam.df.external=cbind(Modelo=mods,R2=gam.df.rsq, RMSE=gam.df.RMSE, MSE=gam.df.MSE)
gam.df.external=as.data.frame(gam.df.external)
gam.df.external## Modelo R2 RMSE MSE
## 1 GAM1 0.68 60.654 3678.933
## 2 GAM2 0.67 62.119 3858.731
## 3 GAM3 0.65 63.041 3974.206
## 4 GAM4 0.53 242.104 58614.289write.csv(gam.df.external, file = "../res/tab/GAM.df.external.csv")# Fazendo alguns plotes com curvas de suavização para algumas covariáveis
p1 <- qplot(x = NDVI, y = Argila, data = dados , geom = c("point", "smooth"))
p2 <- qplot(x = B1, y = Argila, data = dados , geom = c("point", "smooth"))
p3 <- qplot(x = B2, y = Argila, data = dados , geom = c("point", "smooth"))
p4 <- qplot(x = B3, y = Argila, data = dados , geom = c("point", "smooth"))
p5 <- qplot(x = B4, y = Argila, data = dados , geom = c("point", "smooth"))
p6 <- qplot(x = B5, y = Argila, data = dados , geom = c("point", "smooth"))
p7 <- qplot(x = B7, y = Argila, data = dados , geom = c("point", "smooth"))
grid.arrange(p1, p2, p3, p4, p5, p6, p7, nrow=3, ncol = 3)## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'
### Gerando figuras da validação
# FAZER FIGURA DE ESTIMADOS x OBSERVADOS COM R2 #
#GAM1
jpeg(filename = "../res/fig/valida.GAM1.jpg", width = 600, height = 400, units = "px", pointsize = 12,
quality = 100, bg = "white", res = NA, family = "")
par(mfrow=c(1,1))
{plot(valida1~Argila,xlab="Observedo (g kg-1)", data=dados.valida, ylab="Estimado (g kg-1)", main="Argila (GAM1) ", cex.lab=1.5, cex.main=1.5, cex.axis=1.5,pch = 19)
abline(lm(valida1~Argila, data=dados.valida))
text(200,500,"R2 =",cex=1)
text(230,500,round(summary(lm(Argila~valida1, data=dados.valida))$r.squared,2),cex=1)}
dev.off()## png
## 2#GAM3
jpeg(filename = "../res/fig/valida.GAM3.jpg", width = 600, height = 400, units = "px", pointsize = 12,
quality = 100, bg = "white", res = NA, family = "")
par(mfrow=c(1,1))
{plot(valida3~Argila,xlab="Observedo (g kg-1)", data=dados.valida, ylab="Estimado (g kg-1)", main="Argila (GAM3) ", cex.lab=1.5, cex.main=1.5, cex.axis=1.5,pch = 19)
abline(lm(valida3~Argila, data=dados.valida))
text(200,500,"R2 =",cex=1)
text(230,500,round(summary(lm(Argila~valida3, data=dados.valida))$r.squared,2),cex=1)}
dev.off()## png
## 2
Figure 2.3: validação
Figure 2.4: validação
### Gerando o mapa do teor de argila usano o nosso melhor modelo
# importar o grid da area de estudo
gride <- read.dbf("../data/Pontos1.dbf")
gride$X=gride$POINT_X; gride$Y=gride$POINT_Y
names (gride) # ver os nomes das colunas na tabela## [1] "POINTID" "POINT_X" "POINT_Y" "B1" "B2" "B3" "B4"
## [8] "B5" "B7" "B3_B2" "B3_B7" "B5_B7" "NDVI" "X"
## [15] "Y"pred<-predict(GAM1, gride, interval="conf", level=0.95)
predicao=cbind(pred,gride[,2:3]) # este data.frame pode ser salvo como CSV e aberto no ArcGIS
predicao[1:10,] # ver os nomes das colunas na tabela## pred POINT_X POINT_Y
## 1 379.3606 337807.4 8938236
## 2 353.8333 337537.4 8938206
## 3 380.2189 337567.4 8938206
## 4 410.8855 337597.4 8938206
## 5 306.7033 337627.4 8938206
## 6 192.0280 337657.4 8938206
## 7 226.7614 337687.4 8938206
## 8 283.5475 337717.4 8938206
## 9 335.6467 337747.4 8938206
## 10 380.6050 337777.4 8938206summary(predicao$fit)## Length Class Mode
## 0 NULL NULLgridded(predicao) <- ~POINT_X+POINT_Y # transformar os dados de dataframe para spatialpixeldataframe
plot(predicao)
#### PACOTE ESPECIAL PARA EXPORTAR ####
library(rgdal)
writeGDAL(predicao[,1], "../res/fig/Argila.GAM.sdat", drivername="SAGA",mvFlag=-99999)3 A era das máquinas começou
Figure 3.1: Máquina aprendendo
Referências
Chagas, C. da S., Pinheiro, H.S.K., Carvalho Junior, W. de, Anjos, L.H.C. dos, Pereira, N.R., Bhering, S.B., 2017. Data mining methods applied to map soil units on tropical hillslopes in Rio de Janeiro, Brazil. Geoderma Regional 9, 47–55. https://doi.org/10.1016/j.geodrs.2017.03.004
Dotto, A.C., Dalmolin, R.S.D., Pedron, F. de A., Caten, A. ten, Ruiz, L.F.C., 2014. Mapeamento digital de atributos: granulometria e matéria orgânica do solo utilizando espectroscopia de reflectância difusa. Revista Brasileira de Ciencia do Solo 38, 1663–1671.
Hair, J.F., Black, W.C., Babin, B.J., Anderson, R.E., 2009. Multivariate data analysis.
Hastie, T., Tibshirani, R., Friedman, J., 2009. The elements of statistical learning: Data mining, inference and prediction, 2nd ed. Springer Series in Statistics, Stanford, California.
Moura-Bueno, J.M., 2014. Modelos digitais de elevação e predição do carbono orgânico do solo no Planalto do estado do Rio Grande do Sul (PhD thesis).
Samuel-Rosa, A., 2012. Funções de predição espacial de propriedades do solo (PhD thesis).
Wood, S., 2006. Generalized Additive Models: An Introduction with R. CRC Texts in Statistical Science.