Exemplos de Regressão - Forest Fires

Os dados utilizados no exemplo são referentes a incêndios em florestas na região Nordeste de Portugal. A tarefa é prever os tamanhos das áreas incendiadas (em hectares). Os dados podem ser encontrados neste link

Informações sobre os dados

O dataset é composto de 517 observações com 12 variáveis preditoras e uma variável resposta.

Informações sobre atributos:

* X -  coordenada espacial no eixo x
* Y - coordenada espacial no eixo y
* month - mês do ano
* day - dia da semana
* Índices do sistema FWI:
  + FFMC
  + DMC
  + DC
  + ISI
* temp - temperatura em graus celsius
* RH - humidade relativa
* wind - velocidade do vento em km/h
* rain - nível de chuva (em milímetro) por área (em metros quadrados)
* area - área incendiada na floresta (em hectares)

A seguir temos um breve sumário dos dados:

forestfires <- read.csv("~/Documents/ad2/problema3-regressao/forestfire-exemploCaret/forestfires.csv")

summary(forestfires)
##        X               Y           month      day          FFMC      
##  Min.   :1.000   Min.   :2.0   aug    :184   fri:85   Min.   :18.70  
##  1st Qu.:3.000   1st Qu.:4.0   sep    :172   mon:74   1st Qu.:90.20  
##  Median :4.000   Median :4.0   mar    : 54   sat:84   Median :91.60  
##  Mean   :4.669   Mean   :4.3   jul    : 32   sun:95   Mean   :90.64  
##  3rd Qu.:7.000   3rd Qu.:5.0   feb    : 20   thu:61   3rd Qu.:92.90  
##  Max.   :9.000   Max.   :9.0   jun    : 17   tue:64   Max.   :96.20  
##                                (Other): 38   wed:54                  
##       DMC              DC             ISI              temp      
##  Min.   :  1.1   Min.   :  7.9   Min.   : 0.000   Min.   : 2.20  
##  1st Qu.: 68.6   1st Qu.:437.7   1st Qu.: 6.500   1st Qu.:15.50  
##  Median :108.3   Median :664.2   Median : 8.400   Median :19.30  
##  Mean   :110.9   Mean   :547.9   Mean   : 9.022   Mean   :18.89  
##  3rd Qu.:142.4   3rd Qu.:713.9   3rd Qu.:10.800   3rd Qu.:22.80  
##  Max.   :291.3   Max.   :860.6   Max.   :56.100   Max.   :33.30  
##                                                                  
##        RH              wind            rain              area        
##  Min.   : 15.00   Min.   :0.400   Min.   :0.00000   Min.   :   0.00  
##  1st Qu.: 33.00   1st Qu.:2.700   1st Qu.:0.00000   1st Qu.:   0.00  
##  Median : 42.00   Median :4.000   Median :0.00000   Median :   0.52  
##  Mean   : 44.29   Mean   :4.018   Mean   :0.02166   Mean   :  12.85  
##  3rd Qu.: 53.00   3rd Qu.:4.900   3rd Qu.:0.00000   3rd Qu.:   6.57  
##  Max.   :100.00   Max.   :9.400   Max.   :6.40000   Max.   :1090.84  
## 
Com a função str, podemos ver o tipo das variáveis:
str(forestfires)
## 'data.frame':    517 obs. of  13 variables:
##  $ X    : int  7 7 7 8 8 8 8 8 8 7 ...
##  $ Y    : int  5 4 4 6 6 6 6 6 6 5 ...
##  $ month: Factor w/ 12 levels "apr","aug","dec",..: 8 11 11 8 8 2 2 2 12 12 ...
##  $ day  : Factor w/ 7 levels "fri","mon","sat",..: 1 6 3 1 4 4 2 2 6 3 ...
##  $ FFMC : num  86.2 90.6 90.6 91.7 89.3 92.3 92.3 91.5 91 92.5 ...
##  $ DMC  : num  26.2 35.4 43.7 33.3 51.3 ...
##  $ DC   : num  94.3 669.1 686.9 77.5 102.2 ...
##  $ ISI  : num  5.1 6.7 6.7 9 9.6 14.7 8.5 10.7 7 7.1 ...
##  $ temp : num  8.2 18 14.6 8.3 11.4 22.2 24.1 8 13.1 22.8 ...
##  $ RH   : int  51 33 33 97 99 29 27 86 63 40 ...
##  $ wind : num  6.7 0.9 1.3 4 1.8 5.4 3.1 2.2 5.4 4 ...
##  $ rain : num  0 0 0 0.2 0 0 0 0 0 0 ...
##  $ area : num  0 0 0 0 0 0 0 0 0 0 ...
Como podemos ver, existem duas variáveis categóricas (month e day). Iremos removê-las para facilitar nossa vida.
require(dplyr)
forestfires <- forestfires %>% select(-day, -month)

Examinando histogramas das variáveis

require(reshape2)
require(ggplot2)


df <- melt(forestfires)
ggplot(df,aes(x = value)) + 
    facet_wrap(~variable, scales = "free_x") + 
    geom_histogram()

Transformação nos Dados

Ao analisar os histogramas, é interessante observar em quais gráficos existem muita assimetria. Recomenda-se realizar um ajuste nessas variáveis aplicando alguma função que gere dados mais próximos de uma distribuição normal.

forestfires$log.Y <- log(forestfires$Y)
forestfires$log.DC <- log(forestfires$DC + 1)
forestfires$log.RH <- log(forestfires$RH)
forestfires$log.wind <- log(forestfires$wind)
forestfires$log.area <- log(forestfires$area + 1)
df <- melt(forestfires)
ggplot(df,aes(x = value)) + 
    facet_wrap(~variable, scales = "free_x") + 
    geom_histogram()

Particionamento dos dados: Treino e Teste

require(caret)

trainIndex <- createDataPartition(forestfires$log.area, p = .75, list = FALSE)

train <- forestfires[trainIndex,]
test <- forestfires[-trainIndex,]

Treinamento de Regressão Linear

Para realizar o treinamento, deve-se assumir as seguintes proposições:

  - A relação entre preditores e variável resposta é aproximadamente linear.
  - Um bom preditor está consideravelmente bem correlacionado com a variável resposta.
  - Preditores com alta colinearidade geram instabilidade no modelo.
  - Um preditor não deve ser derivado da combinação linear dos demais.
  - Há um número igual ou maior de observações do que preditores (isso tenta evitar a maldição da dimensionalidade).


correlationMatrix <- cor(forestfires)
print(correlationMatrix)
##                     X            Y        FFMC          DMC          DC
## X         1.000000000  0.539548171 -0.02103927 -0.048384178 -0.08591612
## Y         0.539548171  1.000000000 -0.04630755  0.007781561 -0.10117777
## FFMC     -0.021039272 -0.046307546  1.00000000  0.382618800  0.33051180
## DMC      -0.048384178  0.007781561  0.38261880  1.000000000  0.68219161
## DC       -0.085916123 -0.101177767  0.33051180  0.682191612  1.00000000
## ISI       0.006209941 -0.024487992  0.53180493  0.305127835  0.22915417
## temp     -0.051258262 -0.024103084  0.43153226  0.469593844  0.49620805
## RH        0.085223194  0.062220731 -0.30099542  0.073794941 -0.03919165
## wind      0.018797818 -0.020340852 -0.02848481 -0.105342253 -0.20346569
## rain      0.065387168  0.033234103  0.05670153  0.074789982  0.03586086
## area      0.063385299  0.044873225  0.04012200  0.072994296  0.04938323
## log.Y     0.531494601  0.973352433 -0.05172036 -0.000110987 -0.09158812
## log.DC   -0.063975630 -0.070535390  0.34530923  0.653334918  0.94651532
## log.RH    0.080301614  0.055065326 -0.25946727  0.084086452 -0.01430512
## log.wind  0.033862600  0.003378939  0.01494708 -0.069256421 -0.15982208
## log.area  0.061994908  0.038838213  0.04679856  0.067152740  0.06635976
##                   ISI        temp          RH         wind         rain
## X         0.006209941 -0.05125826  0.08522319  0.018797818  0.065387168
## Y        -0.024487992 -0.02410308  0.06222073 -0.020340852  0.033234103
## FFMC      0.531804931  0.43153226 -0.30099542 -0.028484809  0.056701533
## DMC       0.305127835  0.46959384  0.07379494 -0.105342253  0.074789982
## DC        0.229154169  0.49620805 -0.03919165 -0.203465691  0.035860862
## ISI       1.000000000  0.39428710 -0.13251718  0.106825888  0.067668190
## temp      0.394287104  1.00000000 -0.52739034 -0.227116220  0.069490547
## RH       -0.132517177 -0.52739034  1.00000000  0.069410067  0.099751223
## wind      0.106825888 -0.22711622  0.06941007  1.000000000  0.061118880
## rain      0.067668190  0.06949055  0.09975122  0.061118880  1.000000000
## area      0.008257688  0.09784411 -0.07551856  0.012317277 -0.007365729
## log.Y    -0.026017226 -0.03807040  0.06945431 -0.007171111  0.037843476
## log.DC    0.286186198  0.53133294 -0.06596066 -0.171909763  0.033787859
## log.RH   -0.128587013 -0.50179002  0.97586364  0.044159611  0.094891378
## log.wind  0.116384763 -0.15355750  0.04786933  0.949644419  0.059946913
## log.area -0.010346879  0.05348655 -0.05366216  0.066973489  0.023311313
##                  area        log.Y      log.DC      log.RH     log.wind
## X         0.063385299  0.531494601 -0.06397563  0.08030161  0.033862600
## Y         0.044873225  0.973352433 -0.07053539  0.05506533  0.003378939
## FFMC      0.040122004 -0.051720356  0.34530923 -0.25946727  0.014947084
## DMC       0.072994296 -0.000110987  0.65333492  0.08408645 -0.069256421
## DC        0.049383225 -0.091588119  0.94651532 -0.01430512 -0.159822076
## ISI       0.008257688 -0.026017226  0.28618620 -0.12858701  0.116384763
## temp      0.097844107 -0.038070402  0.53133294 -0.50179002 -0.153557500
## RH       -0.075518563  0.069454306 -0.06596066  0.97586364  0.047869333
## wind      0.012317277 -0.007171111 -0.17190976  0.04415961  0.949644419
## rain     -0.007365729  0.037843476  0.03378786  0.09489138  0.059946913
## area      1.000000000  0.035214168  0.05203853 -0.07895492  0.027440629
## log.Y     0.035214168  1.000000000 -0.06798841  0.06262184  0.013881816
## log.DC    0.052038530 -0.067988410  1.00000000 -0.04187272 -0.131526348
## log.RH   -0.078954922  0.062621841 -0.04187272  1.00000000  0.033396499
## log.wind  0.027440629  0.013881816 -0.13152635  0.03339650  1.000000000
## log.area  0.524134036  0.047660541  0.06900635 -0.04964497  0.066917694
##             log.area
## X         0.06199491
## Y         0.03883821
## FFMC      0.04679856
## DMC       0.06715274
## DC        0.06635976
## ISI      -0.01034688
## temp      0.05348655
## RH       -0.05366216
## wind      0.06697349
## rain      0.02331131
## area      0.52413404
## log.Y     0.04766054
## log.DC    0.06900635
## log.RH   -0.04964497
## log.wind  0.06691769
## log.area  1.00000000

Analisando graficamente a matriz de correlação

library(corrplot)
corrplot(correlationMatrix, method="circle", type="lower", order="hclust")

Observe este gráfico atentamente. Variáveis correlacionadas com seus respectivos logs estão altamente relacionadas. Seria algo redundante utilizar ambas para gerar no modelo. Variáveis bem correlacionadas com a variável resposta produzirão modelos com resultados mais expressivos.
É possível perceber que as variáveis preditoras estão fracamente correlacionadas com a variável resposta log.area. Isso não é um bom sinal.

Treinamento de modelo

#parâmetro de ajuste
fitControl <- trainControl(## 10-fold CV
                           method = "cv",
                           number = 10)
lm <- lm(log.area ~ FFMC + temp + DMC + log.DC + log.wind + X + log.Y + RH ,
               data = forestfires)
lm
## 
## Call:
## lm(formula = log.area ~ FFMC + temp + DMC + log.DC + log.wind + 
##     X + log.Y + RH, data = forestfires)
## 
## Coefficients:
## (Intercept)         FFMC         temp          DMC       log.DC  
##    0.472996    -0.001689    -0.007535     0.001324     0.084991  
##    log.wind            X        log.Y           RH  
##    0.212060     0.035626     0.106713    -0.007121
summary(lm)
## 
## Call:
## lm(formula = log.area ~ FFMC + temp + DMC + log.DC + log.wind + 
##     X + log.Y + RH, data = forestfires)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5929 -1.0951 -0.5969  0.9130  5.6667 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.472996   1.309338   0.361   0.7181  
## FFMC        -0.001689   0.013070  -0.129   0.8972  
## temp        -0.007535   0.016585  -0.454   0.6498  
## DMC          0.001324   0.001415   0.935   0.3501  
## log.DC       0.084991   0.091025   0.934   0.3509  
## log.wind     0.212060   0.125435   1.691   0.0915 .
## X            0.035626   0.031521   1.130   0.2589  
## log.Y        0.106713   0.236857   0.451   0.6525  
## RH          -0.007121   0.005073  -1.404   0.1611  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.395 on 508 degrees of freedom
## Multiple R-squared:  0.0204, Adjusted R-squared:  0.004971 
## F-statistic: 1.322 on 8 and 508 DF,  p-value: 0.2297
Residuals: fornecem o valor mínimo, primeiro quartil, mediana
terceiro quartil e valor máximo dos resíduos da regressão
A distribuição dos resíduos, tal como estimada com esses números,
deve ser simétrica, a mediana deveria ser perto de 0, 
Os valores de 1Q e 3Q deveriam ter valor absoluto semelhantes.
Coefficients: fornecem os valores estimados de beta_0 hat e beta_1 hat
(Estimate), o desvio padrão dos valores (Std. Error), 
o valor t (t value) e o p-valor associado Pr(>|t|)
O desvio padrão fornece uma estimativa da incerteza associada aos 
coeficientes do modelo.
A estatística t é igual ao coeficiente dividido pelo desvio padrão.
Neste caso, com p-valor muito baixo, os dois coeficientes são estatisticamente significativos, isto é, diferentes de zero.
Residual standard error: valor estimado do desvio padrão dos resíduos.
R-squared: a fração da variância em log.psa explicada por log.vol
Adjusted R-squared: igual a R-squared com ajuste pela complexidade 
do modelo (o número de parâmetros).
Não tem importância para um modelo tão simples com 1 preditor.
F-statistic: a razão entre a variância explicada pelo modelo e a variância residual ou não explicada.
Se este valor for alto (com p-valor baixo), significa que o modelo explica
mais variância do que a não explicada.

Realizando previsão

prediction <- predict(lm, select(test,FFMC,temp ,DMC ,log.DC,log.wind,X,log.Y,RH))

exp.prediction <- exp(prediction) - 1 #transformando o vetor de predição para a unidade original do problema (em ha).

lm_prediction <- data.frame(pred = exp.prediction, obs = test$area)
ggplot(lm_prediction, aes(x = pred, y = obs)) + geom_point(alpha = 0.5, position = position_jitter(width=0.2)) + geom_abline(colour = "blue") + ggtitle("Previsão x Observado (validação)")

round(defaultSummary(lm_prediction), digits = 3)
##     RMSE Rsquared 
##   21.949    0.002
O RMSE alto indica que o erro na predição foi alto. O que acontece é que, por conta dos outliers, o modelo se ajusta bem a boa parte dos dados tendo como consequência grandes erros na previsão alguns valores. Dependendo do caso, isso pode indicar que a regressão não seria a melhor técnica para resolver o problema.

Resíduos na validação

ggplot(lm_prediction, aes(y = obs - pred, x = pred)) + 
  geom_point(alpha = 0.5, position = position_jitter(width=0.1)) + 
  geom_abline(slope = 0, intercept = 0, colour = "darkred")

Regressão com o pacote Caret: utilizando o modelo LASSO para a seleção de preditores

fitControl <- trainControl(method='cv', number = 10)

lasso.fit <- train(log.area ~ ., data=select(train, FFMC,ISI,temp,DMC,log.DC,log.wind,X,log.Y ,RH, log.area), 
                  method='lasso', 
                  metric="RMSE",
                  tuneLength = 10,
                  trControl=fitControl)

Validação Cruzada: Resultados

plot(lasso.fit)

lasso.fit
## The lasso 
## 
## 389 samples
##   9 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 350, 350, 349, 350, 351, 350, ... 
## Resampling results across tuning parameters:
## 
##   fraction   RMSE      Rsquared    RMSE SD    Rsquared SD
##   0.1000000  1.424426  0.02738409  0.1749171  0.02680433 
##   0.1888889  1.424673  0.02084919  0.1752045  0.01790298 
##   0.2777778  1.425126  0.02005852  0.1755858  0.02333097 
##   0.3666667  1.425356  0.02027669  0.1763439  0.02804242 
##   0.4555556  1.424778  0.02211763  0.1776085  0.03179565 
##   0.5444444  1.423886  0.02457893  0.1793339  0.03484948 
##   0.6333333  1.422838  0.02667298  0.1805959  0.03694805 
##   0.7222222  1.421879  0.02840167  0.1810235  0.03822739 
##   0.8111111  1.421340  0.02984201  0.1815390  0.03938656 
##   0.9000000  1.421776  0.02963991  0.1824419  0.03808917 
## 
## RMSE was used to select the optimal model using  the smallest value.
## The final value used for the model was fraction = 0.8111111.
A função train tentou 10 valores para fraction e achou valor ótimo (RMSE mais baixo).

Importância das Variáveis para o modelo

plot(varImp(lasso.fit))

Realizando previsão com modelo LASSO

lasso_prediction <- predict(lasso.fit, select(test, FFMC,ISI,temp,DMC,log.DC,log.wind,X,log.Y ,RH))

exp.lasso_prediction <- exp(lasso_prediction) - 1 #transformando o vetor de predição para a unidade original do problema (em ha).

lasso_prediction <- data.frame(pred = exp.lasso_prediction, obs = test$area)
round(defaultSummary(lasso_prediction), digits = 3)
##     RMSE Rsquared 
##   21.954    0.004

Comparando os modelos

compare <- lm_prediction
compare$model <- "RL"

lasso_prediction$model <- "LASSO"

compare <- rbind(compare, lasso_prediction)

ggplot(compare, aes(x = pred, y = obs)) + 
  geom_point(alpha = 0.5, position = position_jitter(width=0.2)) + 
  facet_grid(. ~ model) + 
  geom_abline() +
  ggtitle("Observado x Previsão (validação)")

round(defaultSummary(lasso_prediction), digits = 3)
##     RMSE Rsquared 
##   21.954    0.004
round(defaultSummary(lm_prediction), digits = 3)
##     RMSE Rsquared 
##   21.949    0.002

Conclusão

O melhor modelo será aquele que possuir o RMSE mais baixo. Como podemos ver, ambos os modelos produziram um RMSE alto, sugerindo que a regressão linear talvez não seja uma boa solução para este problema. 
Contudo, analisando os gráficos, podemos observar que a regressão se ajustou bem aos valores observados e parece encontrar resultados bem aproximados. Sendo assim, o que causa o alto valor do RMSE?
Podemos culpar os outliers. O RMSE é calculado através da raiz quadrada da média dos erros quadrados. Com isso, quanto maior os erros, maior o RMSE. Como os outliers possuem valores altíssimos (houve um incêndio em 1200 ha, por exemplo), isso faz com que a métrica aumente consideravelmente. 
Dependendo do problema, indica-se remover os outliers, pois estes são considerados comportamentos atípicos nos dados. Neste problema em específico, entendo que o objetivo seja prever os outliers, que seriam catástrofes de grandes áreas incendiadas. Com isso, o indicado seria buscar um modelo alternativo a regressão linear.