ATENÇÃO: o código dessa vinheta não é mais avaliado para reduzir tempo de compilação do pacote. Portanto, ao usar esse código esteja ciente de que ele pode estar defasado.

opts_chunk$set(eval = FALSE)

Análise Exploratória

#-----------------------------------------------------------------------
# Carrega os pacotes necessários.

# rm(list = objects())
library(lattice)
library(latticeExtra)
library(corrplot)
library(EACS)
#-----------------------------------------------------------------------
# Análise exploratório dos dados.

str(teca_gran) # Variáveis granulométricas.
str(teca_qui)  # Variáveis químicas e físicas.
str(teca_arv)  # Contém a produção.

# Junção das informações.
db <- merge(subset(teca_qui, select = c("loc", "cam", "arg", "cas", "acc")),
            # subset(teca_gran, select = -are),
            subset(teca_gran, select = -c(afina, amedi, agros)),
            all = TRUE)
str(db)

# Renomeia níveis.
db$cam <- factor(db$cam, labels = 1:3)

# Empilha.
db_long <- reshape2::melt(data = subset(db, cam != "1"),
                          id.vars = c("loc", "cam"),
                          variable.name = "x",
                          value.name = "y")
db_long <- transform(db_long,
                     resp_cam = paste(x, cam, sep = "_"))
db_long <- droplevels(db_long)
str(db_long)

# Distribuição de densidade.
densityplot(~y | x,
            groups = cam,
            data = db_long,
            scales = "free",
            auto.key = TRUE,
            as.table = TRUE)

ecdfplot(~y | x,
         groups = cam,
         data = db_long,
         scales = "free",
         auto.key = TRUE,
         as.table = TRUE)

# Desempilha.
db_wide <- reshape2::dcast(data = db_long,
                           formula = loc ~ resp_cam,
                           value.var = "y")
# str(db_wide)

# Adiciona a variável resposta.
db_prod <- merge(subset(teca_arv, select = c("loc", "prod")),
                 db_wide,
                 by = "loc",
                 all = TRUE)
str(db_prod)

parallelplot(~db_prod[-1], data = db_prod,
             horizontal.axis = FALSE,
             scales = list(x = list(rot = 90)))

# splom(db_wide[, -1])
r <- cor(db_prod[, -(1)])

corrplot::corrplot(r,
                   type = "upper",
                   tl.pos = "d",
                   outline = TRUE,
                   method = "ellipse")

# TODO: fazer uma tabela de estatísticas descritivas da variável
# resposta e de cada preditora: média, mediana, quartis, desvio-padrão,
# amplitude, CV.
my_stats <- function(x) {
    m <- mean(x)
    fn <- fivenum(x)
    s <- sd(x)
    n <- sum(!is.na(x))
    c("n" = n,
      "Média" = m,
      "Mediana" = fn[3],
      "Mín." = fn[1],
      "Máx." = fn[5],
      "Q1" = fn[2],
      "Q3" = fn[3],
      "Desv. pad." = s,
      "CV%" = 100 * s/m)
}

resum <- do.call(rbind, lapply(db_prod[, -1], FUN = my_stats))
knitr::kable(resum, digits = 2)

O conjunto de dados obtido até aqui db_prod contém as variáveis

  • prod: produção de madeira.
  • arg: o conteúdo de argila do solo.
  • are: o conteúdo de areia do solo.
  • afina, amedi e agros: fracionamento da areia em fina, média e grossa.
  • cas: o conteúdo de cascalho do solo.
  • acc: soma de areia, cascalho e calhau (particulas grandes).

As variáveis granulométricas do solo foram obtidas em 3 camadas do solo. No entando, considerando que a primeira camada é muito fina e mais sujeita às ações de manejo, não será considerada nas análises.

Regressão múltipla com seleção de variáveis por stepwise

A seleção de variáveis por stepwise irá remover/adicionar termos com a finalidade de minimizar a medida de ajuste BIC.

m_full <- lm(prod ~ ., data = db_prod)
summary(m_full)

# Critério de BIC.
m_sel <- step(m_full,
              direction = "both",
              k = log(nrow(db_prod)),
              trace = FALSE)

# Resumo do modelo selecionado.
summary(m_sel)

# par(mfrow = c(2, 2))
# plot(m_sel)
# layout(1)

Controle do processo de particação dos dados

#-----------------------------------------------------------------------

library(caret)

y <- db_prod$prod
X <- as.matrix(db_prod[, 3:ncol(db_prod)])
dim(X)

# Controle da validação cruzada.
set.seed(123)
ctr <- trainControl(method = "repeatedcv",
                    number = 10,
                    repeats = 3,
                    returnResamp = "all",
                    savePredictions = "all")

Regressão com penalização LASSO

A regressão com penalização LASSO será usada para fazer seleção de variáveis. Aquelas variáveis cujos coeficientes forem anulados pela penalização imposta serão removidas e então um modelo linear será ajustado apenas com as variáveis restantes. Será usado validação cruzada para determinar o valor do parâmetro de penalidade.

# str(getModelInfo(model = "glmnet")[[1]])

# Regressão com penalização LASSO.
set.seed(123)
fit_las <- train(y = y,
                 x = X,
                 method = "glmnet",
                 trControl = ctr,
                 tuneGrid = expand.grid(alpha = 1,
                                        lambda = 2^seq(from = -5,
                                                       to = 5,
                                                       by = 0.2)))
fit_las

# Gráfico para escolha do `lambda`.
sc <- list(x = list(log = 2))
c("Traço médio" = plot(fit_las, scale = sc),
  "Desempenhos individuais na CV" =
      xyplot(fit_las, , scale = sc) + latticeExtra::as.layer(plot(fit_las, scale = sc)))

# Resultado do modelo final.
# fit_las$finalModel

# Importância das variáveis.
plot(varImp(fit_las))

# names(fit_las$finalModel)

# Lambda ótimo.
fit_las$finalModel$lambdaOpt

# Variáveis não anuladas pela penalidade.
b <- coef(fit_las$finalModel, s = fit_las$finalModel$lambdaOpt)
b

b <- as.matrix(b)
v <- tail(names(b[abs(b) > 0, ]), n = -1)

# Resumo do modelo com as variáveis selecionadas.
m_las <- lm(prod ~ ., data = db_prod[, c("prod", v)])
summary(m_las)

#--------------------------------------------
# Para fazer o traço das estimativas com a penalização.

# library(glmnet)
#
# # Quando alpha = 0, então é Ridge. Se alpha = 1, então é Lasso.
# m_las <- glmnet(x = X, y = y, alpha = 1, nlambda = 100)
#
# # Traço das estimativas em função do hiperparâmetro de penalização.
# plot(m_las, xvar = "lambda", label = TRUE)
# abline(h = 0, lty = 2)
# abline(v = log(fit_las$finalModel$lambdaOpt), lty = 2, col = 2)
# abline(h = b[v, ], col = "orange", lty = 2)
# points(x = rep(log(fit_las$finalModel$lambdaOpt), length(v)),
#        y = b[v, ])

Melhor subconjunto de preditoras

O melhor subconjunto de variáveis é uma abordagem exaustiva que ajusta todos os modelos com \(k\) variáveis de cada vez. Com cada valor de \(k\), armazena-se o melhor modelo obtido. Depois verifica-se em qual \(k\) está o melhor modelo final. Para a determinação do valor de \(k\) será usada validação cruzada.

# str(getModelInfo(model = "leapForward")[[1]])
# help(rpart, h = "html")

library(leaps)

# Best subset regression.
set.seed(123)
fit_bsr <- train(y = y,
                 x = X,
                 method = "leapForward",
                 trControl = ctr,
                 tuneGrid = expand.grid(nvmax = seq(from = 1,
                                                    to = 8,
                                                    by = 1)))
fit_bsr

# Gráfico para escolha do `k`.
c("Traço médio" = plot(fit_bsr),
  "Desempenhos individuais na CV" =
      xyplot(fit_bsr) + latticeExtra::as.layer(plot(fit_bsr)))

# Resultado do modelo final.
summary(fit_bsr$finalModel)

# Ajustando com k = 1 até 4.
b0 <- regsubsets(prod ~ . - loc,
                 data = db_prod,
                 method = "exhaustive",
                 nvmax = 4)
sel <- summary(b0)

# Variáveis que são mantidas no modelo.
sel <- sel$which[, -1]
sel <- apply(sel,
             MARGIN = 1,
             FUN = function(x) {
                 colnames(sel)[x]
             })
sel

m_bsr <- lm(prod ~ ., data = db_prod[, c("prod", sel[[1]])])
summary(m_bsr)

Árvore de regressão

A árvore de regressão faz partições no espaço das variáveis preditodas gerando uma árvore semelhande as de classificação taxonomica. Cada ramificação da árvore faz decisão com base no valor de corte de uma variável. Os nós terminais contém os valores médios e tamanho de amostra para cada caminho percorrido pela árvore do nó raiz até cada nó terminal. A árvore de regressão possui um parâmetro de complexidade cp será definido por validação cruzada. Quando menor o valor de cp, mais ramificações a árvore terá, levando a um modelo mais flexível sujeito a superajuste. Quanto menor o valor de cp, mais simples é o modelo (com menos partições e nós terminais). Será usado um tamanho mínimo de nó de 5 observações (minsplit, 5 é 10% do número de registros) para ocorrência de partição.

# str(getModelInfo(model = "rpart")[[1]])
# help(rpart, h = "html")
library(rpart)

# Árvore regressão.
set.seed(123)
fit_arv <- train(y = y,
                 x = X,
                 method = "rpart",
                 trControl = ctr,
                 control = rpart.control(minsplit = 5),
                 tuneGrid = expand.grid(cp = seq(from = 0.01,
                                                 to = 0.3,
                                                 by = 0.01)))
fit_arv

# Gráfico para escolha do `cp`.
c("Traço médio" = plot(fit_arv),
  "Desempenhos individuais na CV" =
      xyplot(fit_arv) + latticeExtra::as.layer(plot(fit_arv)))

m_arv <- rpart(prod ~ .,
               data = db_prod,
               control = rpart.control(cp = 0.23,
                                       minsplit = 5))
m_arv

par(mfrow = c(1, 2), xpd = FALSE)
plot(m_arv)
text(m_arv, use.n = TRUE)
plot(prod ~ arg_2, data = db_prod)
points(predict(m_arv) ~ arg_2, data = db_prod, col = 2)
abline(v = 331.6, lty = 2, col = 2)
layout(1)

Regressão de componentes principais

A regressão com componentes principais irá usar como variável regressora um conjunto de índices produzidos com combinações lineares das variáveis preditoras. Os índices sintetizam a informação das regressoras em uma dimensão menor porém prevervando ao máximo a contribuição das preditoras na construção de tais índices. O número de componentes principais (ncomp) será determinado por validação cruzada.

# str(getModelInfo(model = "pcr")[[1]])

# Principal components regression.
set.seed(123)
fit_pcr <- train(y = y,
                 x = X,
                 method = "pcr",
                 trControl = ctr,
                 tuneGrid = expand.grid(ncomp = seq(from = 1,
                                                    to = 8,
                                                    by = 1)))
fit_pcr

# Gráfico para escolha do `ncomp`.
c("Traço médio" = plot(fit_pcr),
  "Desempenhos individuais na CV" =
      xyplot(fit_pcr) + latticeExtra::as.layer(plot(fit_pcr)))

# Resultado do modelo final.
summary(fit_pcr$finalModel)

# Carregamentos para contrução dos índices.
fit_pcr$finalModel$loadings

# Coeficientes de regressão.
fit_pcr$finalModel$coefficients

# biplot(fit_pcr$finalModel)
loa <- fit_pcr$finalModel$loadings[, ]
cbind(loa)

plot(loa, asp = 1, type = "n")
abline(v = 0, h = 0, lty = 3)
arrows(loa[, 1], loa[, 2], 0, 0, length = 0.05, code = 1, col = "orange")
text(loa[, 1], loa[, 2], rownames(loa))

Floresta aleatória

# str(getModelInfo(model = "rf")[[1]])

# Random forest.
set.seed(123)
fit_rfr <- train(y = y,
                 x = X,
                 method = "rf",
                 importance = TRUE,
                 ntree = 100,
                 nodesize = 8,
                 trControl = ctr,
                 tuneGrid = expand.grid(mtry = seq(from = 1,
                                                   to = 5,
                                                   by = 1)))
fit_rfr

# Gráfico para escolha do `mtry`.
c("Traço médio" = plot(fit_rfr),
  "Desempenhos individuais na CV" =
      xyplot(fit_rfr) + latticeExtra::as.layer(plot(fit_rfr)))

# Resultado do modelo final.
fit_rfr$finalModel

# Importância das variáveis.
plot(varImp(fit_rfr, scale = TRUE))
# plot(varImp(fit_rfr, scale = FALSE))

# ip <- importance(fit_rfr$finalModel)
# head(ip[order(ip[, 1], decreasing = TRUE), ], n = 10)

Integração dos resultados

preds <- cbind(resp = db_prod$prod,
               m_sel = predict(m_sel),
               m_bsr = predict(m_bsr),
               m_arv = predict(m_arv),
               pcr = predict(fit_pcr),
               rfr = predict(fit_rfr))

splom(preds, as.matrix = TRUE, type = c("p", "r"))

# Correlação entre preditos e observados.
cbind(sort(cor(preds)[, 1]))

# Importância relativa das variáveis por cada modelo.
vimp <- cbind("Best subset regression" = unname(varImp(fit_bsr)$importance),
              "Lasso penalized regression" = unname(varImp(fit_las)$importance),
              "Random forest regression" = unname(varImp(fit_rfr)$importance),
              "Principal components regression" = unname(varImp(fit_pcr)$importance))
vimp

vimp <- cbind(variable = rownames(vimp), stack(vimp))
# names(vimp)

barchart(reorder(variable, values) ~ values,
         groups = ind,
         data = vimp,
         horizon = TRUE,
         auto.key = TRUE) +
    latticeExtra::layer({
        panel.abline(v = c(0, 100), lty = 2)
    })

barchart(reorder(variable, values) ~ values | ind,
         data = vimp,
         horizon = TRUE,
         layout = c(NA, 2)) +
    latticeExtra::layer({
        panel.abline(v = c(0, 100), lty = 2)
    })

Informações da Sessão