Definições da sessão
#-----------------------------------------------------------------------
# Prof. Dr. Walmes M. Zeviani
# leg.ufpr.br/~walmes · github.com/walmes
# walmes@ufpr.br · @walmeszeviani
# Laboratory of Statistics and Geoinformation (LEG)
# Department of Statistics · Federal University of Paraná
# 2020-jun-14 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
# Pacotes.
library(caret)
library(tidyverse)
library(readxl)
Importação e inspeção
#-----------------------------------------------------------------------
# Importação.
tb <- read_excel("./Fungo teca_W.xlsx")
str(tb)
## tibble [206 × 12] (S3: tbl_df/tbl/data.frame)
## $ Local : num [1:206] 1 1 2 2 3 3 4 4 5 5 ...
## $ Fungo : num [1:206] 1 1 1 1 1 1 1 1 1 1 ...
## $ Camada : chr [1:206] "0-0.1" "0,4-0.6" "0-0.1" "0,4-0.6" ...
## $ COS : num [1:206] NA NA 23.8 NA 23.6 ...
## $ MOS : num [1:206] NA NA 41 NA 40.6 ...
## $ Areia : num [1:206] NA 43.9 58.1 NA 53.7 ...
## $ Silte : num [1:206] NA 13.4 16.6 NA 25.1 ...
## $ Argila : num [1:206] NA 42.7 25.2 NA 21.1 ...
## $ Gradiente textura: num [1:206] NA NA NA NA NA ...
## $ pH : num [1:206] 7.27 7.38 6.7 NA 6.18 ...
## $ Umidade : num [1:206] NA 0.303 0.222 NA 0.321 ...
## $ Gradiente umidade: num [1:206] NA NA NA NA NA ...
names(tb) %>%
str_to_lower() %>%
dput()
## c("local", "fungo", "camada", "cos", "mos", "areia", "silte",
## "argila", "gradiente textura", "ph", "umidade", "gradiente umidade"
## )
names(tb) <- c("local", "fungo", "camada", "cos", "mos", "areia",
"silte", "argila", "gtextura", "ph", "umidade",
"gumidade")
#-----------------------------------------------------------------------
# Inpeção.
# Quantidade de camadas.
tb$camada %>%
unique()
## [1] "0-0.1" "0,4-0.6"
# Percentual de preenchimento das variáveis.
tb %>%
group_by(camada) %>%
summarise_all(~sum(!is.na(.))/length(.)) %>%
ungroup() %>%
gather("var", "val", cos:gumidade) %>%
ggplot(mapping = aes(x = var, y = val, fill = camada)) +
geom_col(position = "dodge") +
labs(x = "Variável", y = "Preenchimento", fill = "Camada") +
scale_fill_manual(values = c("navy", "purple"))
# Criar a tabela analítica de dados.
tb_a <- tb %>%
filter(camada == "0-0.1") %>%
select(-gtextura, -gumidade, -camada)
tb_a
## # A tibble: 103 x 9
## local fungo cos mos areia silte argila ph umidade
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 NA NA NA NA NA 7.27 NA
## 2 2 1 23.8 41.0 58.2 16.6 25.2 6.70 0.222
## 3 3 1 23.6 40.6 53.7 25.1 21.1 6.18 0.321
## 4 4 1 NA NA NA NA NA NA NA
## 5 5 1 22.6 38.9 19.9 47.3 32.8 6.52 0.434
## 6 6 1 23.6 40.6 52.1 24.6 23.3 4.94 0.303
## 7 7 1 NA NA 39.8 20.7 39.5 5.91 0.309
## 8 8 1 21.5 37.0 12.2 52.0 35.8 6.26 0.446
## 9 9 1 23.7 40.9 57.2 21.6 21.2 6.06 0.263
## 10 10 1 18.0 31.1 28.1 51.2 20.6 6.16 0.392
## # … with 93 more rows
tb_b <- tb %>%
filter(camada != "0-0.1") %>%
select(local, gtextura, gumidade)
tb_b
## # A tibble: 103 x 3
## local gtextura gumidade
## <dbl> <dbl> <dbl>
## 1 1 NA NA
## 2 2 NA NA
## 3 3 NA NA
## 4 4 NA NA
## 5 5 1.40 0.837
## 6 6 1.39 1.13
## 7 7 1.03 0.991
## 8 8 1.26 0.950
## 9 9 1.76 1.16
## 10 10 1.31 1.20
## # … with 93 more rows
tbu <- full_join(tb_a, tb_b) %>%
select(-local, -silte, -cos)
## Joining, by = "local"
str(tbu)
## tibble [103 × 8] (S3: tbl_df/tbl/data.frame)
## $ fungo : num [1:103] 1 1 1 1 1 1 1 1 1 1 ...
## $ mos : num [1:103] NA 41 40.6 NA 38.9 ...
## $ areia : num [1:103] NA 58.1 53.7 NA 19.9 ...
## $ argila : num [1:103] NA 25.2 21.1 NA 32.8 ...
## $ ph : num [1:103] 7.27 6.7 6.18 NA 6.52 ...
## $ umidade : num [1:103] NA 0.222 0.321 NA 0.434 ...
## $ gtextura: num [1:103] NA NA NA NA 1.4 ...
## $ gumidade: num [1:103] NA NA NA NA 0.837 ...
Análise exploratória
#-----------------------------------------------------------------------
# Análise exploratória.
tbu %>%
gather(key = "var", value = "val", -fungo) %>%
drop_na() %>%
ggplot(mapping = aes(x = val, y = fungo)) +
facet_wrap(facets = ~var, scale = "free_x") +
geom_point() +
geom_smooth(se = FALSE)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
# Percentual de preenchimento das variáveis.
tbu %>%
summarise_all(~sum(!is.na(.))/length(.)) %>%
ungroup() %>%
gather("var", "val", mos:gumidade) %>%
ggplot(mapping = aes(x = reorder(var, val), y = val)) +
geom_col() +
labs(x = "Variável", y = "Preenchimento")
lattice::splom(tbu)
# Versão sem valores ausentes.
tbuc <- drop_na(tbu)
# Diferença na distribuição com e sem casos completos.
list(dirty = tbu, clean = tbuc) %>%
bind_rows(.id = "condition") %>%
gather("var", "val", -condition) %>%
ggplot(mapping = aes(x = val, color = condition)) +
facet_wrap(facets = ~var, scale = "free") +
geom_density()
Regressão logística
#-----------------------------------------------------------------------
# Regressão logística.
# Ajuste do modelo.
m0 <- glm(fungo ~ ., data = tbuc, family = quasibinomial)
# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Nada é significativo.
summary(m0)
##
## Call:
## glm(formula = fungo ~ ., family = quasibinomial, data = tbuc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7318 -1.1388 0.7280 0.8634 1.6500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.042138 10.068446 -0.501 0.619
## mos 0.071744 0.047091 1.524 0.134
## areia -0.001934 0.059950 -0.032 0.974
## argila -0.068304 0.083765 -0.815 0.418
## ph 0.066258 0.460778 0.144 0.886
## umidade 3.626580 10.785997 0.336 0.738
## gtextura -0.003947 1.574043 -0.003 0.998
## gumidade 3.002416 2.561026 1.172 0.246
##
## (Dispersion parameter for quasibinomial family taken to be 1.080562)
##
## Null deviance: 80.837 on 60 degrees of freedom
## Residual deviance: 71.375 on 53 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 6
summary(m0)$coefficients
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.042138248 10.06844551 -0.500786168 0.6185951
## mos 0.071743525 0.04709085 1.523513093 0.1335755
## areia -0.001933691 0.05995023 -0.032254944 0.9743899
## argila -0.068303537 0.08376452 -0.815423256 0.4184795
## ph 0.066258300 0.46077810 0.143796550 0.8862065
## umidade 3.626580107 10.78599676 0.336230409 0.7380247
## gtextura -0.003947107 1.57404267 -0.002507624 0.9980086
## gumidade 3.002415991 2.56102614 1.172348828 0.2462998
drop1(m0, test = "F")
## Single term deletions
##
## Model:
## fungo ~ mos + areia + argila + ph + umidade + gtextura + gumidade
## Df Deviance F value Pr(>F)
## <none> 71.375
## mos 1 74.094 2.0189 0.1612
## areia 1 71.376 0.0008 0.9771
## argila 1 72.112 0.5472 0.4627
## ph 1 71.397 0.0166 0.8981
## umidade 1 71.497 0.0906 0.7646
## gtextura 1 71.375 0.0000 0.9982
## gumidade 1 73.417 1.5159 0.2237
# Modelo nulo.
m00 <- update(m0, formula = . ~ 1, data = model.frame(m0))
# Razão entre modelos encaixados.
anova(m00, m0, test = "F")
## Analysis of Deviance Table
##
## Model 1: fungo ~ 1
## Model 2: fungo ~ mos + areia + argila + ph + umidade + gtextura + gumidade
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 60 80.837
## 2 53 71.375 7 9.4622 1.251 0.2925
# COMMENT: regressão logística não identificou variáveis relevantes para
# explicar a resposta.
Florestas aleatórias
#-----------------------------------------------------------------------
# Randon forest.
# Cria fator.
tbuc$Fungo <- factor(tbuc$fungo)
# 10 repetições independentes de 5 dobras.
control <- trainControl(method = "repeatedcv",
number = 5,
repeats = 10,
search = "grid")
# Grid do hiperparâmetro.
tunegrid <- expand.grid(.mtry = (1:(ncol(tbu) - 1)))
# Ajusta todos os modelos para selecionar o melhor valor do mtry.
set.seed(2020)
rf_gridsearch <- train(Fungo ~ .,
data = tbuc[, -1],
method = "rf",
metric = "Accuracy",
tuneGrid = tunegrid,
trControl = control)
# Resultados.
print(rf_gridsearch)
## Random Forest
##
## 61 samples
## 7 predictor
## 2 classes: '0', '1'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 10 times)
## Summary of sample sizes: 49, 50, 48, 48, 49, 49, ...
## Resampling results across tuning parameters:
##
## mtry Accuracy Kappa
## 1 0.5076457 -0.14715218
## 2 0.5193124 -0.09621316
## 3 0.5159557 -0.09190607
## 4 0.5123893 -0.09813724
## 5 0.5140326 -0.08794582
## 6 0.5205944 -0.07684311
## 7 0.5153147 -0.08777485
##
## Accuracy was used to select the optimal model using the
## largest value.
## The final value used for the model was mtry = 6.
plot(rf_gridsearch)
# Incluindo o erro padrão.
rf_gridsearch$results %>%
ggplot(mapping = aes(x = mtry, y = Accuracy)) +
geom_point() +
geom_line() +
geom_errorbar(mapping = aes(ymin = Accuracy - AccuracySD,
ymax = Accuracy + AccuracySD),
width = 0.1)
# Usa a raíz do número de variáveis mesmo.
set.seed(3030)
tunegrid <- expand.grid(.mtry = round(sqrt(7), 0))
rf_final <- train(Fungo ~ .,
data = tbuc[, -1],
method = "rf",
metric = "Accuracy",
tuneGrid = tunegrid,
trControl = control)
rf_final
## Random Forest
##
## 61 samples
## 7 predictor
## 2 classes: '0', '1'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 10 times)
## Summary of sample sizes: 48, 49, 49, 49, 49, 48, ...
## Resampling results:
##
## Accuracy Kappa
## 0.4952331 -0.1470569
##
## Tuning parameter 'mtry' was held constant at a value of 3
# Gráfico de importância das variáveis.
# plot(varImp(rf_final, scale = TRUE))
plot(varImp(rf_final, scale = FALSE))
# Matriz de confução (média).
confusionMatrix(rf_final)
## Cross-Validated (5 fold, repeated 10 times) Confusion Matrix
##
## (entries are percentual average cell counts across resamples)
##
## Reference
## Prediction 0 1
## 0 5.7 18.5
## 1 32.0 43.8
##
## Accuracy (average) : 0.4951
# Lista com os indices de treino e teste.
out <- rf_final$control$indexOut
int <- rf_final$control$index
# NOTE: a predição aqui está sendo feita como modelo final que usa todos
# os dados e não uma particação. Dessa forma, a avaliação da performance
# é otimista.
results <-
map_dfr(out,
function(i) {
cm <- confusionMatrix(
data = predict(rf_final,
newdata = tbuc[i, ]),
reference = tbuc[i, ]$Fungo)
as.data.frame(as.list(cm$overall))
})
str(results)
## 'data.frame': 50 obs. of 7 variables:
## $ Accuracy : num 1 1 1 1 1 1 1 1 1 1 ...
## $ Kappa : num 1 1 1 1 1 1 1 1 1 1 ...
## $ AccuracyLower : num 0.753 0.735 0.735 0.735 0.735 ...
## $ AccuracyUpper : num 1 1 1 1 1 1 1 1 1 1 ...
## $ AccuracyNull : num 0.615 0.667 0.583 0.667 0.583 ...
## $ AccuracyPValue: num 0.00182 0.00771 0.00155 0.00771 0.00155 ...
## $ McnemarPValue : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
# Junta com a acurácia feita nos dados de teste.
results <- rf_final$resample %>%
rename_all(~paste0(., "_")) %>%
bind_cols(results) %>%
arrange(Accuracy_) %>%
mutate(index = seq_along(AccuracyNull))
# Acurácia do modelo final contra o modelo nulo.
ggplot(data = results,
mapping = aes(x = index,
xend = index,
y = AccuracyNull,
yend = Accuracy_)) +
geom_segment() +
geom_point(mapping = aes(shape = "Mod. Nulo")) +
geom_point(mapping = aes(y = Accuracy_, shape = "Mod. Ajustado")) +
labs(x = "Reamostragem", y = "Acurácia") +
scale_shape_manual(name = "Performance",
values = c(19, 1))
# Resultado médio.
rf_final$resample %>%
summarise_at("Accuracy", "mean")
## Accuracy
## 1 0.4952331
confusionMatrix(rf_final)
## Cross-Validated (5 fold, repeated 10 times) Confusion Matrix
##
## (entries are percentual average cell counts across resamples)
##
## Reference
## Prediction 0 1
## 0 5.7 18.5
## 1 32.0 43.8
##
## Accuracy (average) : 0.4951
#-----------------------------------------------------------------------