Modelos de Regressão e aplicações no ambiente R

13 a 17 de Abril de 2015 - Manaus - AM
Prof. Dr. Walmes M. Zeviani
Fundação Oswaldo Cruz - FIOCRUZ
Lab. de Estatística e Geoinformação - LEG
Departamento de Estatística - UFPR


Inferência em modelos de regressão linear


Teste de hipótese sobre os parâmetros

##=============================================================================
## Modelos de Regressão e aplicações no ambiente R
##
##   13 a 17 de Abril de 2015 - Manaus/AM
##   Fundação Oswaldo Cruz - FIOCRUZ
## 
##                                                  Prof. Dr. Walmes M. Zeviani
##                                                                LEG/DEST/UFPR
##=============================================================================

##-----------------------------------------------------------------------------
## Definições da sessão.

pkg <- c("lattice", "latticeExtra", "gridExtra", "car", "alr3", "asbio",
         "plyr", "aod", "multcomp", "ellipse", "boot", "wzRfun")

sapply(pkg, require, character.only=TRUE)
##      lattice latticeExtra    gridExtra          car         alr3        asbio 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE 
##         plyr          aod     multcomp      ellipse         boot       wzRfun 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE
trellis.device(color=FALSE)

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

Seja o modelo \[ \begin{aligned} \text{E}(Y|x) = \beta_0+\beta_1 x_1+\beta_2 x_2 +\ldots+ \beta_k x_k. \end{aligned} \] O teste a hipótese sobre os parâmetros, \(\text{H}_0: \beta_i = 0\), permite verificar quais as covariáveis que têm influência sobre o valor esperado da variável observada. O teste de hipótese pode ser de várias formas:

  • Teste sobre um único parâmetro. É um teste som um grau de liberdade. Considera a distribuição t.
  • Teste sobre mais de um parâmetro, ou teste conjunto. É um teste com mais de um grau de liberdade. Considera a distribuição F.
  • Teste sobre uma função linear dos parâmetros. Aqui se encontram os testes para valores preditos e contrastes entre parâmentros ou valores preditos. Considera F se for uma hipótese conjunta e t se for uma hipótese simples.
##-----------------------------------------------------------------------------
## Journal of Applied Ecology, vol. 32, 1995.

url <- "http://www.leg.ufpr.br/~walmes/data/
business_economics_dataset/EXERCISE/SNOWGEES.DAT"

da <-
    read.table(paste0(strwrap(url), collapse=""), header=FALSE)
da[,1] <- NULL
names(da) <- c("diet", "wc", "de", "adf")
str(da)
## 'data.frame':    42 obs. of  4 variables:
##  $ diet: Factor w/ 2 levels "Chow","Plants": 2 2 2 2 2 2 2 2 2 2 ...
##  $ wc  : num  -6 -5 -4.5 0 2 3.5 -2 -2.5 -3.5 -2.5 ...
##  $ de  : num  0 2.5 5 0 0 1 2.5 10 20 12.5 ...
##  $ adf : num  28.5 27.5 27.5 32.5 32 30 34 36.5 28.5 29 ...
## Data on gosling weight change (wc), digestion efficiency (de),
## acid-detergent fiber (adf) (all in %) and diet (plants or duck show)
## for 42 feeding trials. The botanists were interested in predicting
## weight change as a function of the other variables.

pairs(da[,-1])

##-----------------------------------------------------------------------------
## Ajuste do modelo.

## m0 <- lm(wc~de+adf+de:adf, data=da)
m0 <- lm(wc~de+adf, data=da)
par(mfrow=c(2,2)); plot(m0); layout(1)

##-----------------------------------------------------------------------------
## Inferência sobre cada beta.

## summary(m0)
summary(m0)$coefficients
##                Estimate Std. Error    t value     Pr(>|t|)
## (Intercept) 12.18043771  4.4023590  2.7667979 0.0086100803
## de          -0.02653926  0.0534876 -0.4961759 0.6225553001
## adf         -0.45783404  0.1282758 -3.5691367 0.0009688975
## Como chegar a tais valores?
X <- model.matrix(m0); head(X)
##   (Intercept)  de  adf
## 1           1 0.0 28.5
## 2           1 2.5 27.5
## 3           1 5.0 27.5
## 4           1 0.0 32.5
## 5           1 0.0 32.0
## 6           1 1.0 30.0
vcov <- solve(crossprod(X))*summary(m0)$sigma^2; vcov
##             (Intercept)           de          adf
## (Intercept)  19.3807647 -0.221914245 -0.551723317
## de           -0.2219142  0.002860923  0.006038197
## adf          -0.5517233  0.006038197  0.016454693
vcov(m0) ## É a matriz de variância e covariância dos betas.
##             (Intercept)           de          adf
## (Intercept)  19.3807647 -0.221914245 -0.551723317
## de           -0.2219142  0.002860923  0.006038197
## adf          -0.5517233  0.006038197  0.016454693
sqrt(diag(vcov))
## (Intercept)          de         adf 
##   4.4023590   0.0534876   0.1282758
## Valor da estatística t.
coef(m0)/sqrt(diag(vcov(m0)))
## (Intercept)          de         adf 
##   2.7667979  -0.4961759  -3.5691367
2*pt(2.7667979, df=df.residual(m0), lower.tail=FALSE)
## [1] 0.008610081
##-----------------------------------------------------------------------------
## Teste por análise de variância.

## Sequêncial.
anova(m0)
## Analysis of Variance Table
## 
## Response: wc
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## de         1 384.24  384.24  31.020  2.05e-06 ***
## adf        1 157.79  157.79  12.739 0.0009689 ***
## Residuals 39 483.08   12.39                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Marginal.
## drop1(m0, test="F")
drop1(m0, scope=.~., test="F")
## Single term deletions
## 
## Model:
## wc ~ de + adf
##        Df Sum of Sq    RSS    AIC F value    Pr(>F)    
## <none>              483.08 108.59                      
## de      1      3.05 486.13 106.85  0.2462 0.6225553    
## adf     1    157.79 640.88 118.46 12.7387 0.0009689 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(m0, type="II")
## Anova Table (Type II tests)
## 
## Response: wc
##           Sum Sq Df F value    Pr(>F)    
## de          3.05  1  0.2462 0.6225553    
## adf       157.79  1 12.7387 0.0009689 ***
## Residuals 483.08 39                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(m0, type="III")
## Anova Table (Type III tests)
## 
## Response: wc
##             Sum Sq Df F value    Pr(>F)    
## (Intercept)  94.82  1  7.6552 0.0086101 ** 
## de            3.05  1  0.2462 0.6225553    
## adf         157.79  1 12.7387 0.0009689 ***
## Residuals   483.08 39                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## m1 <- lm(cbind(wc, wc^2)~de+adf, data=da)
## anova(m1)
## summary(m1)

Os testes t retornados no summary são marginais, ou seja, são fazem restrição aos demais parâmetros. Sendo assim, o teste para verificar se de tem efeito nulo é um teste feito considera adf presente no modelo (e não considerando-o igual a zero ou outro particular valor). Por outro lado, o teste retornado pela função anova() é uma lista de testes sequênciais. No caso, testa df considerando como referência o modelo com intercepto (sem adf, efeito zero). Na sequência, o teste para adf é feito com relação ao modelo de referência que contém de. A função drop1 também faz testes marginais, mas considerando a distribuição F.

##-----------------------------------------------------------------------------
## Testes conjuntos.

## Testar se simultaneamente o efeito de de e adf são zero. Ao
## considerar isso, tem-se que o modelo sem covariáveis (só com
## intercepto) é o modelo de referência. Esse teste pode ser feito de
## várias formas.

summary(m0)
## 
## Call:
## lm(formula = wc ~ de + adf, data = da)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.0649 -2.0241  0.5645  2.4590  6.8556 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12.18044    4.40236   2.767 0.008610 ** 
## de          -0.02654    0.05349  -0.496 0.622555    
## adf         -0.45783    0.12828  -3.569 0.000969 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.519 on 39 degrees of freedom
## Multiple R-squared:  0.5288, Adjusted R-squared:  0.5046 
## F-statistic: 21.88 on 2 and 39 DF,  p-value: 4.25e-07
m00 <- update(m0, formula=.~1)
anova(m00, m0)
## Analysis of Variance Table
## 
## Model 1: wc ~ 1
## Model 2: wc ~ de + adf
##   Res.Df     RSS Df Sum of Sq     F   Pr(>F)    
## 1     41 1025.12                                
## 2     39  483.08  2    542.03 21.88 4.25e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Teste de Wald.
coef(m0)
## (Intercept)          de         adf 
## 12.18043771 -0.02653926 -0.45783404
C <- cbind(0, diag(2))
linearHypothesis(m0, hypothesis.matrix=C, rhs=c(0,0), test="F")
## Linear hypothesis test
## 
## Hypothesis:
## de = 0
## adf = 0
## 
## Model 1: restricted model
## Model 2: wc ~ de + adf
## 
##   Res.Df     RSS Df Sum of Sq     F   Pr(>F)    
## 1     41 1025.12                                
## 2     39  483.08  2    542.03 21.88 4.25e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## linearHypothesis(m0, hypothesis.matrix=C, rhs=c(1,1), test="F")

wald.test(Sigma=vcov(m0), b=coef(m0),
          Terms=2:3, H0=c(0,0), df=df.residual(m0))
## Wald test:
## ----------
## 
## Chi-squared test:
## X2 = 43.8, df = 2, P(> X2) = 3.1e-10
## 
## F test:
## W = 21.9, df1 = 2, df2 = 39, P(> W) = 4.2e-07
## wald.test(Sigma=vcov(m0), b=coef(m0),
##           L=C, H0=c(0,0), df=df.residual(m0))

##-----------------------------------------------------------------------------
## Testes de hipótese simples com correção para multiplicidade.

## Sem correção. Equivalente ao summary().
cftest(m0)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = wc ~ de + adf, data = da)
## 
## Linear Hypotheses:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept) == 0 12.18044    4.40236   2.767 0.008610 ** 
## de == 0          -0.02654    0.05349  -0.496 0.622555    
## adf == 0         -0.45783    0.12828  -3.569 0.000969 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
## Com correção.
g <- glht(m0, linfct=C)
summary(g)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = wc ~ de + adf, data = da)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)   
## 1 == 0 -0.02654    0.05349  -0.496  0.75768   
## 2 == 0 -0.45783    0.12828  -3.569  0.00155 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
summary(g, test=adjusted(type="none"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = wc ~ de + adf, data = da)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)    
## 1 == 0 -0.02654    0.05349  -0.496 0.622555    
## 2 == 0 -0.45783    0.12828  -3.569 0.000969 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)
summary(g, test=adjusted(type="bonferroni"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = wc ~ de + adf, data = da)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)   
## 1 == 0 -0.02654    0.05349  -0.496  1.00000   
## 2 == 0 -0.45783    0.12828  -3.569  0.00194 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)
summary(g, test=adjusted(type="fdr"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = wc ~ de + adf, data = da)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)   
## 1 == 0 -0.02654    0.05349  -0.496  0.62256   
## 2 == 0 -0.45783    0.12828  -3.569  0.00194 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
summary(m0)$coef
##                Estimate Std. Error    t value     Pr(>|t|)
## (Intercept) 12.18043771  4.4023590  2.7667979 0.0086100803
## de          -0.02653926  0.0534876 -0.4961759 0.6225553001
## adf         -0.45783404  0.1282758 -3.5691367 0.0009688975
## Teste conjunto.
summary(g, test=Ftest())
## 
##   General Linear Hypotheses
## 
## Linear Hypotheses:
##        Estimate
## 1 == 0 -0.02654
## 2 == 0 -0.45783
## 
## Global Test:
##       F DF1 DF2   Pr(>F)
## 1 21.88   2  39 4.25e-07
##-----------------------------------------------------------------------------
## Qual a diferença entre testes marginais e conjuntos?

ci <- cbind(coef(m0), confint(m0))

plot(ellipse(m0, which=c("de", "adf"), level=0.95), type="l")
abline(v=ci["de",1], h=ci["adf",1], lty=2)
abline(v=ci["de",-1], h=ci["adf",-1], col=2, lty=2)

\[ \begin{aligned} \text{H}_0 &: \begin{pmatrix} 0 & 1 & 0\\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \beta_0\\ \beta_1\\ \beta_2 \end{pmatrix} = \begin{pmatrix} m_1\\ m_2 \end{pmatrix}\\ &: C\beta = M \\ \hat{V}_{\beta} &= \hat{\sigma}^2 (X^\top X)^{-1}\\ F &= \frac{(C\hat{\beta}-M)^\top (C\hat{V}_{\beta}C^\top)^{-1} (C\hat{\beta}-M)}{\text{rank}(C)} \stackrel{H_0}{\sim} F_{\text{rank}(C), n-p} \end{aligned} \]

##-----------------------------------------------------------------------------
## Como fazer esse teste matricialmente? Chegar ao valor de F?

b <- coef(m0); V <- vcov(m0); M <- c(0,0)

CB <- (C%*%b-M)
CVi <- solve(C%*%V%*%t(C))
F <- t(CB)%*%CVi%*%CB/nrow(C); F
##          [,1]
## [1,] 21.87958
## Aplicativos que fazem essas contas certamente empregam métodos
## númericos, como decomposições de matrizes, ao invés de fazer as
## contas matriciais como colocadas acima.

Intervalos de confiança para os parâmetros

http://socserv.socsci.mcmaster.ca/jfox/Courses/SPIDA/linear-models-notes.pdf

O intervalo de confiança pára um parâmetro é dado por \[ \begin{aligned} \frac{\hat{\beta}_i-\beta_i}{\sqrt{\text{diag}(\hat{V}_{\beta})_i}} &\stackrel{H_0}{\sim} t_{n-p} \\ \text{IC}(\beta_i) &= \left\{\beta_i : |\beta_i-\hat{\beta_i}|\leq t_{\alpha/2}\sqrt{\text{diag}(\hat{V}_{\beta})_i}\right\} \end{aligned} \]

A região de confiança é dada por \[ \begin{aligned} F &= \frac{(\hat{\beta}_{1:k}-\beta_{1:k})^\top (\hat{V}_{1:k})^{-1} (\hat{\beta}_{1:k}-\beta_{1:k})}{k} \stackrel{H_0}{\sim} F_{k, n-p}\\ RC(\beta_{1:k}) &= \{\beta_{1:k} : F\leq F_{k,n-p}\}. \end{aligned} \]

O intervalo de confiança para uma função linear dos parâmetros é obtido por \[ \begin{aligned} F &= (L\hat{\beta})^\top (L\hat{V}_{\beta}L^\top)^{-1} (L\hat{\beta}) \stackrel{H_0}{\sim} F_{1, n-p} = t^2_{n-p}\\ \text{IC}(L\beta) &= L\hat{\beta} \pm t_{\alpha/2, n-p}\sqrt{L\hat{V}_{\beta}L^\top}. \end{aligned} \] O caso mais frequente de função linear dos parâmetros é a estimação dos valores preditos pelo modelo. Nesse caso \(L = x_i^{+}\) representa um vetor com valores das covariáveis para o qual se deseja obter o valor esperado. O intervalo de confiança para o valor futuro é dado por \[ \text{IC}(L\beta) = L\hat{\beta} \pm t_{\alpha/2, n-p}\sqrt{\hat{\sigma}^2+L\hat{V}_{\beta}L^\top} \]


Representando o ajuste com bandas de confiança e de predição

##-----------------------------------------------------------------------------
## Dados de inventário florestal.

rm(list=ls())

dap <- read.table("http://www.leg.ufpr.br/~walmes/data/dap.txt",
                  header=TRUE, sep="\t")
str(dap)
## 'data.frame':    991 obs. of  2 variables:
##  $ DAP: num  4.87 7.38 5.95 5.73 6.4 ...
##  $ HT : num  9.5 9.8 13 13.5 13.5 13.5 13.5 14.3 14.8 14.8 ...
## Renomeia.
names(dap) <- c("d","h")

## 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),]

## reseta o nome das 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"))

##-----------------------------------------------------------------------------
## Ajuste de modelos.

## Linear simples.
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
## Resultado.
layout(matrix(c(1,1,2,3,4,5),2,3))
plot(h~d, dapcc)
lines(fitted(m0)~d, dapcc, col=2)
plot(m0)

layout(1)

##-----------------------------------------------------------------------------
## Medidas de influência.

inf <- influence.measures(m0)
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_*
## Sinalizando os pontos influentes.
## str(inf)

## Influentes pelo DFITS (influência sobre ajuste).
dfits <- inf$is.inf[,4]

plot(h~d, dapcc)
lines(fitted(m0)~d, dapcc, col=2)
with(dapcc, points(d[dfits], h[dfits], col=2, pch=19))

##-----------------------------------------------------------------------------
## Identificar/remover os pontos discrepantes/influentes manualmente
## (critério visual).

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
##-----------------------------------------------------------------------------
## Predição.

er <- extendrange(dapcc$d, f=0.1)
d.new <- seq(er[1], er[2], length=100)
c(head(d.new), tail(d.new))
##  [1]  2.606910  2.881236  3.155562  3.429888  3.704214  3.978540 28.393560 28.667886
##  [9] 28.942212 29.216538 29.490864 29.765190
pred <- data.frame(d=d.new)

## 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':    100 obs. of  7 variables:
##  $ d  : num  2.61 2.88 3.16 3.43 3.7 ...
##  $ fit: num  4.83 5.71 6.54 7.32 8.06 ...
##  $ lwr: num  2.45 3.5 4.49 5.41 6.29 ...
##  $ upr: num  7.21 7.92 8.59 9.23 9.83 ...
##  $ FIT: num  4.83 5.71 6.54 7.32 8.06 ...
##  $ LWR: num  0.672 1.651 2.563 3.416 4.219 ...
##  $ UPR: num  8.99 9.77 10.52 11.23 11.9 ...
## 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(#corner=c(0.05,0.95),
            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]))

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)

p1+as.layer(p2)+as.layer(p3)+
    layer(panel.text(
        x=20, y=10,
        label=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])), bty="n"))