#-----------------------------------------------------------------------
# Load packages.
library(lattice)
library(latticeExtra)
library(asbio)
library(car)
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)")
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\)).
## 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