Conjunto de dados: trees

Exploração preliminar do conjunto de dados

Detalhes do conjunto de dados:

help(trees)

Conhecendo o conjunto de dados carregado:

dim(trees)
## [1] 31  3
str(trees)
## 'data.frame':    31 obs. of  3 variables:
##  $ Girth : num  8.3 8.6 8.8 10.5 10.7 10.8 11 11 11.1 11.2 ...
##  $ Height: num  70 65 63 72 81 83 66 75 80 75 ...
##  $ Volume: num  10.3 10.3 10.2 16.4 18.8 19.7 15.6 18.2 22.6 19.9 ...
head(trees)
##   Girth Height Volume
## 1   8.3     70   10.3
## 2   8.6     65   10.3
## 3   8.8     63   10.2
## 4  10.5     72   16.4
## 5  10.7     81   18.8
## 6  10.8     83   19.7

Análise gráfica univariada

Relação entre a variável Volume (resposta) com Height e Girth (explicativas) , fig.align = ‘center’

  # relação entre volume e height
plot(Volume ~ Height, data = trees, pch = 16, col = "blue", 
    main = "Árvores de Black Cherry", 
    xlab = "Altura, em pés", ylab = "Volume, em pés cúbicos")
erre <- round(with(trees, cor(Volume, Height)), digits = 3)
text(70, 70, paste0("r = ", erre), cex = 0.8)
#
  # relação entre volume e girth
plot(Volume ~ Girth, data = trees, pch = 16, col = "blue", 
    main = "Árvores de Black Cherry", 
    xlab = "Diâmetro, em polegadas", ylab = "Volume, em pés cúbicos")
erre <- round(with(trees, cor(Volume, Girth)), digits = 3)
text(11, 70, paste0("r = ", erre), cex = 0.8)

A variável Girth está fortemente correlacionada com a resposta

Estimação Bootstrap

Objetivo: Escolher modelo de regressão, tendo como critério o coeficiente de determinação (\(R^2\)).

Tabela com os modelos possíveis, incluindo interação entre explicativas.

Modelos de regressão linear
# Modelo Explicativas Expressão
m1 Height Volume ~ Height
m2 Girth Volume ~ Girth
m3 razao = Girth/Height Volume ~ razao
m4 Height e Girth Volume ~ Height + Girth
m5 Height, Girth e razao Volume ~ Height + Girth + razao
Função para estimação dos valores de interesse na amostra bootstrap
volume.estimado <- function(dados, indices){
    d <- dados[indices, ]
    d$razao <- d$Girth / d$Height
      # Height
    m1 <- lm(Volume ~ Height, data = d)
    m1.R2 <- summary(m1)$r.square
      # Girth
    m2 <- lm(Volume ~ Girth, data = d)
    m2.R2 <- summary(m2)$r.square
      # Razão Girth/Height
    m3 <- lm(Volume ~ razao, data = d)
    m3.R2 <- summary(m3)$r.square
      # Height e Girth
    m4 <- lm(Volume ~ Height + Girth, data = d)
    m4.R2 <- summary(m4)$r.square
      # Height, Girth e razão
    m5 <- lm(Volume ~ Height + Girth + razao, data = d)
    m5.R2 <- summary(m5)$r.square
      # vetor com estimativas de R2
    relacoes <- c(m1.R2, m2.R2, m3.R2, m4.R2, m5.R2)
    return(relacoes)
}
Estimação bootstrap

Estimação usando o pacote boot, com semente para permitir replicação dos resultados

library(boot)
set.seed(666)

Cálculo do \(R^2\) de cada modelo usando B = 5.000 reamostras. Os cálculos demandarão algum tempo.

  # Estimativa bootstrap - 5.000 iterações
resultados <- boot(data = trees, statistic = volume.estimado, R = 5000)

resultados é objeto de classe boot, tendo a seguinte estrutura

## List of 11
##  $ t0       : num [1:5] 0.358 0.935 0.731 0.948 0.973
##  $ t        : num [1:5000, 1:5] 0.419 0.434 0.484 0.427 0.542 ...
##  $ R        : num 5000
##  $ data     :'data.frame':   31 obs. of  3 variables:
##   ..$ Girth : num [1:31] 8.3 8.6 8.8 10.5 10.7 10.8 11 11 11.1 11.2 ...
##   ..$ Height: num [1:31] 70 65 63 72 81 83 66 75 80 75 ...
##   ..$ Volume: num [1:31] 10.3 10.3 10.2 16.4 18.8 19.7 15.6 18.2 22.6 19.9 ...
##  $ seed     : int [1:626] 403 624 1904417942 -1425571297 -658651692 -120115515 1618310338 -1447143077 -246272800 -1745531039 ...
##  $ statistic:function (dados, indices)  
##   ..- attr(*, "srcref")= 'srcref' int [1:8] 1 20 22 1 20 1 1 22
##   .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x000000000ac7c7d0> 
##  $ sim      : chr "ordinary"
##  $ call     : language boot(data = trees, statistic = volume.estimado, R = 5000)
##  $ stype    : chr "i"
##  $ strata   : num [1:31] 1 1 1 1 1 1 1 1 1 1 ...
##  $ weights  : num [1:31] 0.0323 0.0323 0.0323 0.0323 0.0323 ...
##  - attr(*, "class")= chr "boot"
##  - attr(*, "boot_type")= chr "boot"

Saída com os resultados das estimativas para cada um dos cinco modelos, com informações sobre a estimativas bootstrap de \(R^2\), do viés da estimativa e de seu erro padrão.

print(resultados)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = trees, statistic = volume.estimado, R = 5000)
## 
## 
## Bootstrap Statistics :
##      original       bias    std. error
## t1* 0.3579026 0.0027022293  0.11693873
## t2* 0.9353199 0.0003293783  0.01753864
## t3* 0.7309204 0.0029201793  0.08064173
## t4* 0.9479500 0.0031197047  0.01193377
## t5* 0.9732894 0.0004290413  0.01032291

Plot do modelo 1: Volume ~ Height

plot(resultados, index = 1)

Obtenção do intervalo bootstrap com 95% de confiança para o \(R^2\) do modelo m1 (Volume ~ Height)

  # Intervalos bootstrap
IC.m1.obj <- boot.ci(resultados, index = 1, conf = 0.95, type = "bca")
print(IC.m1.obj)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 5000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = resultados, conf = 0.95, type = "bca", index = 1)
## 
## Intervals : 
## Level       BCa          
## 95%   ( 0.1333,  0.5929 )  
## Calculations and Intervals on Original Scale
ic.m1 <- IC.m1.obj$bca[ , c(4, 5)]
print(ic.m1)
##                     
## 0.1333111 0.5929402

Obtenção dos intervalos bootstrap com 95% para o \(R^2\) para todos os modelos:

(intervalos <- t(
            sapply(1:5, FUN = function(x) boot.ci(resultados, index = x, 
                conf = 0.95, type = "bca")$bca[ , c(4, 5)])
            )
)
##                         
## [1,] 0.1333111 0.5929402
## [2,] 0.8762088 0.9581551
## [3,] 0.4776446 0.8434214
## [4,] 0.9033394 0.9642847
## [5,] 0.9400533 0.9867468

Histograma das estimativas bootstrap do modelo 1, com densidade estimada (azul) e limites do intervalo de confiança bootstrao

  # histograma, suavização densidade (blue), IC (red)
hist(resultados$t[, 1], main = "Coeficiente de Determinação: Height", 
    xlab = expression(R^2), ylab = "Densidade", col = "grey", prob = T)
lines(density(resultados$t[, 1]), col = "blue")
abline(v = ic.m1, col = "red")

Recuperação das estimativas boostrap de \(R^2\) do modelo m1 (Volume ~ Height):

  # vetor com valores de estatistica das amostras bootstrap
head(resultados$t[, 1])
## [1] 0.4190074 0.4337527 0.4840151 0.4269644 0.5415744 0.1198099

Tabela com as estimativas bootstrap e intervalares de \(R^2\) de todos os modelos:

tabela <- cbind(resultados$t0, intervalos)
colnames(tabela) <- c("R2", "IC.inf", "IC.sup")
rownames(tabela) <- paste0("m", 1:5)
tabela
##           R2    IC.inf    IC.sup
## m1 0.3579026 0.1333111 0.5929402
## m2 0.9353199 0.8762088 0.9581551
## m3 0.7309204 0.4776446 0.8434214
## m4 0.9479500 0.9033394 0.9642847
## m5 0.9732894 0.9400533 0.9867468

Qual modelo escolher? Apenas coeficiente de determinação é um bom critério de escolha do modelo? Repetir para o cirtério AIC (ou BIC) e verificar se a conlusão é a mesma.