Medidas de Influência e Bandas de Confiança
Carregando Pacotes
#-----------------------------------------------------------------------
# Load packages.
library(lattice)
library(latticeExtra)
library(asbio)
library(car)Carregando os Dados
Os dados referem-se a valores de altura (h, m) e diâmetro a altura do peito (d, DAP, cm) de uma amostra de árvores medidas para estimação do volume de madeira de uma floresta. As análises aqui mostradas vão ajustar um modelo para altura como função do DAP.
#-----------------------------------------------------------------------
# Dados de inventário florestal.
dap <- read.table("http://www.leg.ufpr.br/~walmes/data/dap.txt",
header = TRUE, sep = "\t")
names(dap) <- c("d", "h")
str(dap)## 'data.frame': 991 obs. of 2 variables:
## $ d: num 4.87 7.38 5.95 5.73 6.4 ...
## $ h: num 9.5 9.8 13 13.5 13.5 13.5 13.5 14.3 14.8 14.8 ...# Ordenar para evitar efeito espaguete quando plotar.
dap <- dap[order(dap$d), ]
# Nova base que contém d e h observados.
dapcc <- dap[complete.cases(dap), ]
# Renomeia as linhas
rownames(dapcc) <- NULL
head(dapcc)## d h
## 1 4.8701 9.5
## 2 5.7296 13.5
## 3 5.9524 13.0
## 4 6.3344 15.5
## 5 6.3980 13.5
## 6 6.7482 13.5#-----------------------------------------------------------------------
# Visualização.
xyplot(h ~ d,
data = dap,
type = c("p", "smooth", "g"),
xlab = "Diâmetro à altura do peito (cm)",
ylab = "Altura (m)")
Ajuste do Modelo Hipsométrico
Um modelo hipsométrico utilizado, às vezes por sua simplicidade, é o modelo linear com o termo raíz quadrada. Após ajustar o modelo, seão obtidas as medidas de influência das observações.
#-----------------------------------------------------------------------
# Linear com raíz.
m0 <- lm(h ~ d + sqrt(d), data = dapcc)
summary(m0)##
## Call:
## lm(formula = h ~ d + sqrt(d), data = dapcc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.3596 -1.1310 0.0444 1.1666 10.8313
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.0753 4.2832 -3.987 9.12e-05 ***
## d -1.2162 0.2889 -4.210 3.73e-05 ***
## sqrt(d) 15.3125 2.2487 6.809 9.22e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.099 on 220 degrees of freedom
## Multiple R-squared: 0.7843, Adjusted R-squared: 0.7823
## F-statistic: 399.9 on 2 and 220 DF, p-value: < 2.2e-16# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
#-----------------------------------------------------------------------
# Medidas de influência.
inf <- influence.measures(m0)
# Sinalizando os pontos influentes.
summary(inf)## Potentially influential observations of
## lm(formula = h ~ d + sqrt(d), data = dapcc) :
##
## dfb.1_ dfb.d dfb.sq() dffit cov.r cook.d hat
## 1 -0.25 -0.22 0.24 -0.27 1.17_* 0.02 0.14_*
## 2 0.12 0.11 -0.12 0.14 1.11_* 0.01 0.09_*
## 3 -0.01 0.00 0.01 -0.01 1.10_* 0.00 0.08_*
## 4 0.19 0.16 -0.18 0.22 1.07_* 0.02 0.06_*
## 5 -0.04 -0.03 0.04 -0.05 1.08_* 0.00 0.06_*
## 6 -0.09 -0.08 0.08 -0.11 1.07_* 0.00 0.05_*
## 7 -0.03 -0.03 0.03 -0.04 1.07_* 0.00 0.05_*
## 8 -0.12 -0.10 0.11 -0.15 1.06_* 0.01 0.05_*
## 9 0.10 0.08 -0.09 0.12 1.06_* 0.00 0.05_*
## 10 -0.41 -0.33 0.37 -0.56_* 0.94_* 0.10 0.04
## 11 -0.04 -0.03 0.04 -0.05 1.05_* 0.00 0.04
## 15 0.38 0.25 -0.31 0.78_* 0.71_* 0.18 0.02
## 29 -0.02 0.03 -0.01 -0.25 0.96_* 0.02 0.01
## 41 -0.17 -0.26 0.23 0.58_* 0.69_* 0.10 0.01
## 160 0.02 0.00 -0.01 -0.20 0.94_* 0.01 0.01
## 194 -0.08 -0.11 0.09 -0.27 0.93_* 0.02 0.01
## 209 -0.25 -0.31 0.28 -0.50_* 0.87_* 0.08 0.02
## 221 0.03 0.04 -0.04 0.05 1.05_* 0.00 0.03
## 222 0.13 0.15 -0.14 0.19 1.04_* 0.01 0.04
## 223 -0.19 -0.22 0.20 -0.26 1.05_* 0.02 0.05_*# Influentes pelo DFFITS.
dfits <- inf$is.inf[, 4]
er <- extendrange(dapcc$d, f = 0.05)
pred <- data.frame(d = seq(er[1], er[2], length.out = 30))
pred$y <- predict(m0, newdata = pred)
plot(h ~ d, dapcc)
lines(y ~ d, data = pred, col = 2)
with(dapcc, points(d[dfits], h[dfits], col = 2, pch = 19))
Dentre os vários critérios de influência calculados, desenvolvi simpatia pelo DFFits por ele indicar observações influêntes que saltam aos olhos nos gráficos. Então essa correspondência visual me parece interessante. O DFFits expressa a influência de uma observação sobre o seu valor predito, o que é uma interpretação interessante também.
As observações influentes serão removidas e o ajuste será refeito.
#-----------------------------------------------------------------------
id <- which(dfits)
# Refazer a análise com os pontos removidos.
m1 <- update(m0, data = dapcc[-id, ])
summary(m1)##
## Call:
## lm(formula = h ~ d + sqrt(d), data = dapcc[-id, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.6150 -1.0316 0.1181 1.2315 3.9940
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -15.210 3.585 -4.243 3.28e-05 ***
## d -1.034 0.241 -4.291 2.68e-05 ***
## sqrt(d) 14.081 1.878 7.497 1.66e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.729 on 216 degrees of freedom
## Multiple R-squared: 0.8459, Adjusted R-squared: 0.8444
## F-statistic: 592.7 on 2 and 216 DF, p-value: < 2.2e-16# Fazendo predição com intervalo de confiança e predição futura.
Yp <- predict(m1, newdata = pred, interval = "confidence")
Yf <- predict(m1, newdata = pred, interval = "prediction")
colnames(Yf) <- toupper(colnames(Yf))
pred <- cbind(pred, Yp, Yf)
str(pred)## 'data.frame': 30 obs. of 8 variables:
## $ d : num 3.74 4.6 5.46 6.31 7.17 ...
## $ y : num 7.98 10.16 12.05 13.72 15.21 ...
## $ fit: num 8.15 10.23 12.04 13.64 15.08 ...
## $ lwr: num 6.39 8.83 10.92 12.75 14.37 ...
## $ upr: num 9.91 11.63 13.15 14.53 15.8 ...
## $ FIT: num 8.15 10.23 12.04 13.64 15.08 ...
## $ LWR: num 4.32 6.54 8.45 10.12 11.6 ...
## $ UPR: num 12 13.9 15.6 17.2 18.6 ...#-----------------------------------------------------------------------
# Inclusão da expressão do modelo.
c0 <- coef(m1)
co <- format(c(abs(c0), summary(m1)$r.squared), digits = 3, trim = TRUE)
sinal <- ifelse(c0 < 0, "-", "+")
co[seq_along(sinal)] <- paste(sinal, co[seq_along(sinal)], sep = "")
l <- c("Predito", "IC predito", "IC obs. futura")
key <- list(lines = list(lty = 1),
text = list(l[1]),
rect = list(col = c("red"), alpha = 0.1, lty = 3),
text = list(l[2]),
rect = list(col = c("blue"), alpha = 0.1, lty = 3),
text = list(l[3]))
source(paste0("https://raw.githubusercontent.com/",
"walmes/wzRfun/master/R/panel.cbH.R"))
p1 <- xyplot(h ~ d, data = dapcc[-id, ], xlim = er, xlab = "DAP (cm)",
ylab = "Altura (m)", key = key)
p2 <- xyplot(fit ~ d, data = pred, type = "l", ly = pred$lwr,
uy = pred$upr, cty = "bands", fill = "red",
prepanel = prepanel.cbH, panel = panel.cbH)
p3 <- xyplot(FIT ~ d, data = pred, type = "l", ly = pred$LWR,
uy = pred$UPR, cty = "bands", fill = "blue",
prepanel = prepanel.cbH, panel = panel.cbH)
eqn <- substitute(hat(h) == b0 ~ b1 * d ~ b2 * sqrt(d) ~~~~ (R^2 == r2),
list(b0 = co[1], b1 = co[2], b2 = co[3], r2 = co[4]))
p1 +
as.layer(p2) +
as.layer(p3) +
layer(panel.text(x = 20, y = 10, label = eqn, bty = "n"))
No gráfico foram adicionadas a curva predita, a banda de confiança para o valor predito (\(\hat{y}\)) e a banda de confiança para a observação futura (\(Y\)).
Informações da Sessão
## Updated in 2016-10-10.
##
## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.1 LTS
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=pt_BR.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] tcltk stats graphics grDevices utils datasets
## [7] base
##
## other attached packages:
## [1] car_2.1-2 asbio_1.3-4 latticeExtra_0.6-28
## [4] RColorBrewer_1.1-2 lattice_0.20-33 rmarkdown_1.0
## [7] knitr_1.14
##
## loaded via a namespace (and not attached):
## [1] pixmap_0.4-11 Rcpp_0.12.7 magrittr_1.5
## [4] splines_3.3.1 MASS_7.3-45 scatterplot3d_0.3-37
## [7] minqa_1.2.4 stringr_1.0.0 tools_3.3.1
## [10] parallel_3.3.1 nnet_7.3-12 pbkrtest_0.4-6
## [13] grid_3.3.1 nlme_3.1-128 mgcv_1.8-13
## [16] quantreg_5.26 multcompView_0.1-7 plotrix_3.6-2
## [19] MatrixModels_0.4-1 htmltools_0.3.5 lme4_1.1-12
## [22] yaml_2.1.13 digest_0.6.10 Matrix_1.2-6
## [25] nloptr_1.0.4 formatR_1.4 deSolve_1.13
## [28] evaluate_0.9 stringi_1.1.1 methods_3.3.1
## [31] SparseM_1.7 mvtnorm_1.0-5