Capítulo 7: Modelo de Regressão Linear Multivariada

DES4033: Análise Multivariada em Biologia

Bastão de Asclépio & Distribuição Normal

invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))

options(warn=-1)
suppressMessages(library(apaTables, warn.conflicts=FALSE))
suppressMessages(library(aplpack, warn.conflicts=FALSE))
suppressMessages(library(calculus, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(cellWise, warn.conflicts=FALSE))
suppressMessages(library(cluster, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(diptest, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts=FALSE))
suppressMessages(library(ellipse, warn.conflicts=FALSE))
suppressMessages(library(far, warn.conflicts=FALSE))
suppressMessages(library(FlexReg, warn.conflicts=FALSE))
suppressMessages(library(fmsb, warn.conflicts=FALSE))
suppressMessages(library(GGally, warn.conflicts=FALSE))
suppressMessages(library(ggfortify, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(HSAUR2, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(lawstat, warn.conflicts=FALSE))
suppressMessages(library(lfda, warn.conflicts=FALSE))
suppressMessages(library(MASS, warn.conflicts=FALSE))
suppressMessages(library(MatchIt, warn.conflicts=FALSE)) # distancia robusta
suppressMessages(library(matlib, warn.conflicts=FALSE))
suppressMessages(library(matrixcalc, warn.conflicts=FALSE))
suppressMessages(library(MatrixModels, warn.conflicts=FALSE))
suppressMessages(library(matrixStats, warn.conflicts=FALSE))
suppressMessages(library(MomTrunc, warn.conflicts=FALSE))
suppressMessages(library(MVA, warn.conflicts=FALSE)) # An Introduction to Applied Multivariate Analysis with R
suppressMessages(library(MVar.pt, warn.conflicts=FALSE))
# suppressMessages(library(MVLM, warn.conflicts=FALSE))
suppressMessages(library(MVN, warn.conflicts=FALSE))
suppressMessages(library(mvtnorm, warn.conflicts=FALSE))
suppressMessages(library(olsrr, warn.conflicts=FALSE))
suppressMessages(library(plotly, warn.conflicts=FALSE))
suppressMessages(library(ppcor, warn.conflicts=FALSE))
suppressMessages(library(pracma, warn.conflicts=FALSE))
suppressMessages(library(rcompanion, warn.conflicts=FALSE))
suppressMessages(library(reticulate, warn.conflicts=FALSE))
suppressMessages(library(rgl, warn.conflicts=FALSE))
suppressMessages(library(scatterplot3d, warn.conflicts=FALSE))
source("summarySEwithin2.R")

Material

• HTML de R Markdown em RPubs
• Prof. José Siqueira: ResearchGate

Pensamento

“Não é paradoxo dizer que nos nossos momentos de inspiração mais teórica podemos estar o mais próximo possível de nossas aplicações mais práticas.”

WHITEHEAD, AN apud BOYER, CB (1974) História da matemática. São Paulo: Blücher/EDUSP, p. 419.

Sumário

1. Aspectos de análise multivariada
2. Álgebra matricial e vetor estocástico
3. Geometria amostral e amostragem aleatória
4. Distribuição normal multivariada
5. Inferência sobre vetor de média
6. Comparação de várias médias multivariadas
7. Modelo de regressão linear multivariada
8. Componentes principais
9. Análise de fatores e inferência sobre matriz de covariância
10. Análise de correlação canônica
11. Discriminação e classificação
12. Clusterização, métodos de distância e ordination

Incluir

• Mallow’s Cp: olsrr::ols_mallows_cp

full_model <- lm(mpg ~ ., data = mtcars)
model <- lm(mpg ~ disp + hp + wt + qsec, data = mtcars)
ols_mallows_cp(model, full_model)

• Sawa’s bayesian information criterion: olsrr::ols_sbic

• Classical Multivariate Regression

• Alexopoulos E. C. (2010). Introduction to multivariate regression analysis. Hippokratia, 14(Suppl 1), 23–28. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3049417/

• Notes for Predictive Modeling: 4.3 Multivariate multiple linear model: Eduardo García-Portugués, 2023: https://bookdown.org/egarpor/PM-UC3M/lm-iii-mult.html

• The Epidemiologist: R Handbook: 19 Univariate and multivariable regression: https://epirhandbook.com/en/univariate-and-multivariable-regression.html

Pré-requisitos do Capítulo 7

• Análise de Correlação em R

• Análise de Regressão Linear Simples & Regressão de Deming em R

Introdução

A análise de regressão é a metodologia estatística para prever valores de uma ou mais variáveis de resposta (dependentes) a partir de uma coleção de valores de variáveis preditoras (independentes). Também pode ser usada para avaliar os efeitos das variáveis preditoras nas respostas. Infelizmente, o nome “regressão”, derivado do título do primeiro artigo sobre o assunto por F. Galton, de forma alguma reflete a importância ou a amplitude de aplicação desta metodologia.

Neste capítulo, discutimos primeiramente o modelo de regressão múltipla para a previsão de uma única resposta. Esse modelo é então generalizado para lidar com a previsão de várias variáveis dependentes. Nosso tratamento deve ser um tanto conciso, pois existe uma vasta literatura sobre o assunto. (Se você estiver interessado em aprofundar-se na análise de regressão, consulte os seguintes livros, em ordem ascendente de dificuldade: Abraham e Ledolter, Bowerman e O’Connell, Kutner et al. (2005), Draper e Smith (1998), Cook e Weisberg, Seber e Goldberger.) Nosso tratamento abreviado destaca as suposições da regressão e suas consequências, formulações alternativas do modelo de regressão e a aplicabilidade geral das técnicas de regressão em situações aparentemente diferentes.

Modelo Clássico de Regressão Linear Múltipla (MLR)

Seja \(z_1, z_2, \ldots, z_r\) as variáveis preditoras que se acredita estarem relacionadas com uma variável de resposta \(Y\). Por exemplo, com \(r = 4\), poderíamos ter o seguinte modelo de regressão hedônica:

• \(Y\): valor de mercado atual do imóvel

e

• \(z_1\): metros quadrados de área habitável
• \(z_2\): localização (indicador para a zona da cidade central ou não central)
• \(z_3\): valor avaliado no ano passado
• \(z_4\): qualidade da construção (preço por metro quadrado)

O modelo de regressão linear clássico afirma que \(Y\) é composto por uma média, que depende de maneira contínua das variáveis \(z_i\), e um erro aleatório \(e\), que contabiliza o erro de medição e os efeitos de outras variáveis não explicitamente consideradas no modelo. Os valores das variáveis preditoras registrados no experimento ou definidos pelo pesquisador são tratados como fixos (sem erro de medição). O erro (e, portanto, a resposta) é visto como uma variável aleatória cujo comportamento é caracterizado por um conjunto de pressupostos distribucionais.

Erro de Medição (Kutner et al. 2005, p. 165)

Em modelo de regressão não consideramos explicitamente a presença de erros de medição nas observações, seja na variável de resposta \(Y\) ou na variável preditora \(X\). Agora examinaremos brevemente os efeitos dos erros de medição nas observações das variáveis de resposta e preditora.

Erro de Medição na VD

Quando erros de medição aleatórios estão presentes nas observações da variável de resposta \(Y\), não são criados novos problemas quando esses erros são não correlacionados e não tendenciosos (erros de medição positivos e negativos tendem a se anular). Considere, por exemplo, um estudo da relação entre o tempo necessário para completar uma tarefa (\(Y\)) e a complexidade da tarefa (\(X\)). O tempo para completar a tarefa pode não ser medido com precisão porque a pessoa operando o cronômetro pode não fazê-lo nos instantes precisos necessários. Contanto que esses erros de medição sejam de natureza aleatória, não correlacionados e não tendenciosos, esses erros de medição são simplesmente absorvidos no termo de erro do modelo \(\varepsilon\). O termo de erro do modelo sempre reflete os efeitos compostos de um grande número de fatores não considerados no modelo, um dos quais agora seria a variação aleatória devido à imprecisão no processo de medição de \(Y\).

Erro de Medição na VI

Infelizmente, uma situação diferente ocorre quando as observações na variável preditora \(X\) estão sujeitas a erros de medição. Frequentemente, é claro, as observações em \(X\) são precisas, sem erros de medição, como quando a variável preditora é o preço de um produto em diferentes lojas, o número de variáveis em diferentes problemas de otimização ou a taxa salarial para diferentes classes de empregados. Em outros momentos, no entanto, erros de medição podem entrar no valor observado para a variável preditora, por exemplo, quando a variável preditora é a pressão em um tanque, a temperatura em um forno, a velocidade de uma linha de produção ou a idade relatada de uma pessoa.

Regressão Linear Univariada Simples

• \(r\): número de variáveis explicativas/preditoras/indepedentes/covariáveis intervalares

• \(p=1, q=1, r=1\)

Tabela 1. Representações univariada e multivariada do GLM (Modelo Linear General). Chartier & Faulkner, 2008

Uma forma de definir com mais precisão o número de variáveis dependentes (DV) é \(pq\). Se \(pq=1\), o modelo é univariado. Se \(pq\ge2\), o modelo é multivariado.

Observar que Correlação Canônica não tem VD-VI, pois não é um modelo de regressão. A confusão da classificação começa ao misturar variável manifesta com variável latente. Variável latente pode ser nominal ou intervalar (continuous). Note que Regressão Multivariada em GLM multivariado.

Figura 1. Diagramas de Venn para a) dois preditores intervalares em regressão, b) análise de variância unifatorial, c) ANOVA bifatorial e d) medidas repetidas. Chartier & Faulkner, 2008

• Chartier, S & Faulkner, A (2008) General Linear Models: An integrated approach to statistics. Tutorial in Quantitative Methods for Psychology 4(2): 65‐78.

• JW6_Cap7_RLS_ANOVA_Testet.R

invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
alfa <- 0.05 
# The Effect of Vitamin C on Tooth Growth in Guinea Pigs
# len: numeric  Tooth length
# supp: factor  Supplement type (VC or OJ).
#   dose:   numeric Dose in milligrams/day
Dados <- datasets::ToothGrowth
print.data.frame(Dados)

    len supp dose
1   4.2   VC  0.5
2  11.5   VC  0.5
3   7.3   VC  0.5
4   5.8   VC  0.5
5   6.4   VC  0.5
6  10.0   VC  0.5
7  11.2   VC  0.5
8  11.2   VC  0.5
9   5.2   VC  0.5
10  7.0   VC  0.5
11 16.5   VC  1.0
12 16.5   VC  1.0
13 15.2   VC  1.0
14 17.3   VC  1.0
15 22.5   VC  1.0
16 17.3   VC  1.0
17 13.6   VC  1.0
18 14.5   VC  1.0
19 18.8   VC  1.0
20 15.5   VC  1.0
21 23.6   VC  2.0
22 18.5   VC  2.0
23 33.9   VC  2.0
24 25.5   VC  2.0
25 26.4   VC  2.0
26 32.5   VC  2.0
27 26.7   VC  2.0
28 21.5   VC  2.0
29 23.3   VC  2.0
30 29.5   VC  2.0
31 15.2   OJ  0.5
32 21.5   OJ  0.5
33 17.6   OJ  0.5
34  9.7   OJ  0.5
35 14.5   OJ  0.5
36 10.0   OJ  0.5
37  8.2   OJ  0.5
38  9.4   OJ  0.5
39 16.5   OJ  0.5
40  9.7   OJ  0.5
41 19.7   OJ  1.0
42 23.3   OJ  1.0
43 23.6   OJ  1.0
44 26.4   OJ  1.0
45 20.0   OJ  1.0
46 25.2   OJ  1.0
47 25.8   OJ  1.0
48 21.2   OJ  1.0
49 14.5   OJ  1.0
50 27.3   OJ  1.0
51 25.5   OJ  2.0
52 26.4   OJ  2.0
53 22.4   OJ  2.0
54 24.5   OJ  2.0
55 24.8   OJ  2.0
56 30.9   OJ  2.0
57 26.4   OJ  2.0
58 27.3   OJ  2.0
59 29.4   OJ  2.0
60 23.0   OJ  2.0

Dados$supp <- factor(Dados$supp,
                     levels=c("VC", "OJ"))
Dados$dose <- factor(Dados$dose)
boxplot(len~supp,
        data=Dados)

alfaBonf <- alfa/length(unique(Dados$supp))
gplots::plotmeans(len~supp,
                  connect=FALSE,
                  col="black",
                  barcol="black",
                  p=1-alfaBonf,
                  main="Porquinho-da-Índia",
                  data=Dados)

Registered S3 method overwritten by 'gplots':
  method         from     
  reorder.factor DescTools

# Solução 1: lm(len~supp, data=Dados), car::Anova e summary
fit <- lm(len~supp,
          data=Dados)
print(sai <- car::Anova(fit))

Anova Table (Type II tests)

Response: len
          Sum Sq Df F value  Pr(>F)  
supp       205.4  1  3.6683 0.06039 .
Residuals 3246.9 58                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(out <- summary(fit))

Call:
lm(formula = len ~ supp, data = Dados)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.7633  -5.7633   0.4367   5.5867  16.9367 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   16.963      1.366  12.418   <2e-16 ***
suppOJ         3.700      1.932   1.915   0.0604 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.482 on 58 degrees of freedom
Multiple R-squared:  0.05948,   Adjusted R-squared:  0.04327 
F-statistic: 3.668 on 1 and 58 DF,  p-value: 0.06039

R2 <- out$r.squared
df1 <- sai$Df[1]
df2 <- sai$Df[2]
F_omnibus <- (R2/df1)/((1-R2)/df2) # Pestana & Gageiro, 2005, p. 77
pv <- 1-pf(F_omnibus, df1, df2)
if(pv<2.2e-16) p_omnibus <- "< 2.2e-16"
if(pv>2.2e-16) p_omnibus <- paste0("= ", formatC(pv, format="e", digits=2))
Fcrit <- qf(1-alfa, df1, df2)
cat("Omnibus F(",df1, ",", df2, ") = ", round(Fcrit,4), 
    ", F = ", round(F_omnibus,4), 
    ", p ", p_omnibus, ", R^2 = eta^2 = ", round(R2,4), sep="")

Omnibus F(1,58) = 4.0069, F = 3.6683, p = 6.04e-02, R^2 = eta^2 = 0.0595

print(effectsize::eta_squared(fit,
                              partial = TRUE,
                              generalized = FALSE,
                              ci = 1-alfa,
                              alternative = "two.sided",
                              verbose = TRUE), 
      digits=2)

For one-way between subjects designs, partial eta squared is equivalent
  to eta squared. Returning eta squared.

# Effect Size for ANOVA

Parameter | Eta2 |       95% CI
-------------------------------
supp      | 0.06 | [0.00, 0.21]

# Solução 2: lm(Dados$len~supp.num) e model.matrix(~supp, data=Dados) 
supp.num <- as.numeric(Dados$supp)
xtabs(~Dados$supp)

Dados$supp
VC OJ 
30 30 

xtabs(~supp.num)

supp.num
 1  2 
30 30 

fit <- lm(Dados$len~supp.num)
print(car::Anova(fit))

Anova Table (Type II tests)

Response: Dados$len
          Sum Sq Df F value  Pr(>F)  
supp.num   205.3  1  3.6683 0.06039 .
Residuals 3246.9 58                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(summary(fit))

Call:
lm(formula = Dados$len ~ supp.num)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.7633  -5.7633   0.4367   5.5867  16.9367 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   13.263      3.055   4.342 5.73e-05 ***
supp.num       3.700      1.932   1.915   0.0604 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.482 on 58 degrees of freedom
Multiple R-squared:  0.05948,   Adjusted R-squared:  0.04327 
F-statistic: 3.668 on 1 and 58 DF,  p-value: 0.06039

r <- cor(Dados$len, supp.num)
r

[1] 0.2438927

R2 <- r^2
dm <- model.matrix(~supp, data=Dados)
print(dm)

   (Intercept) suppOJ
1            1      0
2            1      0
3            1      0
4            1      0
5            1      0
6            1      0
7            1      0
8            1      0
9            1      0
10           1      0
11           1      0
12           1      0
13           1      0
14           1      0
15           1      0
16           1      0
17           1      0
18           1      0
19           1      0
20           1      0
21           1      0
22           1      0
23           1      0
24           1      0
25           1      0
26           1      0
27           1      0
28           1      0
29           1      0
30           1      0
31           1      1
32           1      1
33           1      1
34           1      1
35           1      1
36           1      1
37           1      1
38           1      1
39           1      1
40           1      1
41           1      1
42           1      1
43           1      1
44           1      1
45           1      1
46           1      1
47           1      1
48           1      1
49           1      1
50           1      1
51           1      1
52           1      1
53           1      1
54           1      1
55           1      1
56           1      1
57           1      1
58           1      1
59           1      1
60           1      1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$supp
[1] "contr.treatment"

n <- dim(dm)[1]
df1 <- dim(dm)[2] - 1
df2 <- dim(dm)[1] - 2
# n <- nrow(Dados)
# df1 <- 1
# df2 <- n - 2
F_omnibus <- (R2/df1)/((1-R2)/df2) # Pestana & Gageiro, 2005, p. 77
pv <- 1-pf(F_omnibus, df1, df2)
if(pv<2.2e-16) p_omnibus <- "< 2.2e-16"
if(pv>2.2e-16) p_omnibus <- paste0("= ", formatC(pv, format="e", digits=2))
Fcrit <- qf(1-alfa, df1, df2)
cat("Omnibus F(",df1, ",", df2, ") = ", round(Fcrit,4), 
    ", F = ", round(F_omnibus,4), 
    ", p ", p_omnibus, ", R^2 = eta^2 = ", round(R2,4), "\n", sep="")

Omnibus F(1,58) = 4.0069, F = 3.6683, p = 6.04e-02, R^2 = eta^2 = 0.0595

print(effectsize::eta_squared(fit,
                              partial = TRUE,
                              generalized = FALSE,
                              ci = 1-alfa,
                              alternative = "two.sided",
                              verbose = TRUE), digits=2)

For one-way between subjects designs, partial eta squared is equivalent
  to eta squared. Returning eta squared.

# Effect Size for ANOVA

Parameter | Eta2 |       95% CI
-------------------------------
supp.num  | 0.06 | [0.00, 0.21]

confint(fit)

                 2.5 %    97.5 %
(Intercept)  7.1490593 19.377607
supp.num    -0.1670064  7.567006

len.media <- mean(Dados$len)
VC.media <- mean(supp.num)
len.sd <- sd(Dados$len)
VC.sd <- sd(supp.num)
beta1_hat <- r*len.sd/VC.sd
beta0.hat <- len.media - beta1_hat*VC.media
y.hat <- beta0.hat + beta1_hat*supp.num
MSE <- sum((Dados$len-y.hat)^2)/(df2)
sqrt(MSE)

[1] 7.482001

beta1.se <- sqrt(MSE/((n-1)*(VC.sd)^2))
t.obs <- beta1_hat/beta1.se
beta1.LI <- beta1_hat - qt(1-alfa/2, df2)*beta1.se
beta1.LS <- beta1_hat + qt(1-alfa/2, df2)*beta1.se
cat("\nIC95%(beta1) = [", round(beta1.LI,4), ", ", 
    round(beta1.LS,4), "]", "\n", sep="")

IC95%(beta1) = [-0.167, 7.567]

confint(fit)

                 2.5 %    97.5 %
(Intercept)  7.1490593 19.377607
supp.num    -0.1670064  7.567006

# Solução 3: DescTools::TTestA
len.mean <- aggregate(len~supp, FUN=mean, data=Dados)
len.sd <- aggregate(len~supp, FUN=sd, data=Dados)
len.n <- aggregate(len~supp, FUN=length, data=Dados)
DescTools::TTestA(mx=len.mean[2,2],
                  sx=len.sd[2,2],
                  nx=len.n[2,2],
                  my=len.mean[1,2],
                  sy=len.sd[1,2],
                  ny=len.n[1,2],
                  paired=FALSE,
                  var.equal=TRUE,
                  conf.level=1-alfa)

    Two Sample t-test

data:  x and y
t = 1.9153, df = 58, p-value = 0.06039
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1670064  7.5670064
sample estimates:
mean of x mean of y 
 20.66333  16.96333 

# Solução 4: Fórmula com n, média e desvio-padrão
nA <- len.n[2,2]
nB <- len.n[1,2]
mediaA <- len.mean[2,2]
mediaB <- len.mean[1,2]
dpA <- len.sd[2,2]
dpB <- len.sd[1,2]
dif <- mediaA - mediaB
dfA <- nA - 1
dfB <- nB - 1
df <- dfA + dfB
sdp <- sqrt((dfA*dpA^2+dfB*dpB^2)/df)
ep <- sdp*sqrt((1/nA+1/nB))
t <- dif/ep
p <- 2*pt(-abs(t),df)
eta2 <- t^2/(t^2 + df)
mag_eta2 <- effectsize::interpret_eta_squared(eta2)
cat("\nAnálise de significância estatística: valor p\n")

Análise de significância estatística: valor p

cat("\tt = ",t,"\n",sep="")

    t = 1.915268

cat("\tdf = ",df,"\n",sep="")

    df = 58

cat("\tValor-p = ",p,"\n",sep="")

    Valor-p = 0.06039337

cat("Análise de significância prática: tamanho de efeito\n")

Análise de significância prática: tamanho de efeito

cat("\tEta^2 = ",eta2, "\n", sep="")

    Eta^2 = 0.05948365

print(mag_eta2)

[1] "small"
(Rules: field2013)

\[\Diamond\]

Modelo de Regressão Linear Univariada Múltipla (UMLR)

Especificamente, o modelo de regressão linear univariada múltipla com uma única resposta tem a forma:

\[ \begin{align} Y &= \sum_{i=0}^{r}{\beta_iz_i} + \varepsilon\\ z_0&=1 \end{align} \]

\[ \text{Resposta} = \text{média (dependendo de } z_1, z_2, \ldots, z_r) + \text{erro} \]

O termo “linear” refere-se ao fato de que a média é uma função linear dos parâmetros desconhecidos \(\beta_0, \beta_1, \ldots, \beta_r\). As variáveis preditoras podem ou não entrar no modelo como termos de primeira ordem.

Com \(n\) observações independentes em \(Y\) e os valores associados de \(z_i\), o modelo completo torna-se:

\[ \begin{aligned} Y_1 & = \beta_0z_{10} + \beta_1z_{11} + \beta_2z_{12} + \cdots + \beta_rz_{1r} + \varepsilon_1 \\ Y_2 & = \beta_0z_{20} + \beta_1z_{21} + \beta_2z_{22} + \cdots + \beta_rz_{2r} + \varepsilon_2 \\ & \vdots \\ Y_n & = \beta_0z_{n0} + \beta_1z_{n1} + \beta_2z_{n2} + \cdots + \beta_rz_{nr} + \varepsilon_n \end{aligned} \]

ou

\[ Y_i = \sum_{k=0}^{r} \beta_k z_{ik} + \varepsilon_i\\ i = 1, 2, \ldots, n \]

Em que \(z_{10}=z_{20}=\cdots=z_{n0}=1\).

Os termos de erro são assumidos ter as seguintes propriedades:

1. \(\mathbb{E}(\varepsilon_i) = 0\);
2. \(\mathbb{V}(\varepsilon_i) = \sigma^2\) (homocedasticidade); e
3. \(\mathbb{C}(\varepsilon_i, \varepsilon_j) = 0\), para \(i \neq j\).

Em notação de matriz, a equação (7-1) torna-se:

\[ \begin{bmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{bmatrix} = \begin{bmatrix} 1 & z_{11} & z_{12} & \cdots & z_{1r} \\ 1 & z_{21} & z_{22} & \cdots & z_{2r} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & z_{n1} & z_{n2} & \cdots & z_{nr} \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_r \end{bmatrix} + \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{bmatrix} \tag{7-1} \]

ou

\[ \underset{n \times 1}{\mathbf{Y}} = \underset{n \times (r+1)}{\mathbf{z}}\underset{(r+1) \times 1}{\boldsymbol{\beta}} + \underset{n \times 1}{\boldsymbol{\varepsilon}} \]

E as especificações em (7-2) tornam-se:

1. \(\mathbb{E}(\boldsymbol{\varepsilon}) = \mathbf{0}\);
2. \(\mathbb{C}(\boldsymbol{\varepsilon}) = \mathbb{E}(\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^{\prime}) = \sigma^2\mathbf{I}\).

	Estocástico	Determinístico
Observável	\(\mathbf{Y}\)	\(\mathbf{z}\)
Estimável	\(\boldsymbol{\varepsilon}\)	\(\boldsymbol{\beta}\)

Observe que o valor 1 na primeira coluna da matriz de planejamento (design matrix) \(\mathbf{z}\) é o multiplicador do termo constante \(\beta_0\). É costumeiro introduzir a variável artificial igual a 1.

Cada coluna de \(\mathbf{z}\) consiste nos \(n\) valores da variável preditora correspondente, enquanto a \(i\)-ésima linha de \(\mathbf{z}\) contém os valores para todas as variáveis preditoras na \(i\)-ésima observação.

Modelo Clássico de Regressão Linear Univariada Múltipla (UMLR)

\[ \begin{align} \mathbf{Y} &= \mathbf{z}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\\\\ \mathbb{E}(\boldsymbol{\varepsilon}) &= \mathbf{0}\\ \mathbb{C}(\boldsymbol{\varepsilon}) &= \sigma^2\mathbf{I} \end{align} \tag{7-3} \]

em que \(\boldsymbol{\beta}\) e \(\sigma^2\) são parâmetros desconhecidos, e a matriz de planejamento \(\mathbf{z}\) tem como linha \(i\):

\[ [z_{i0} \; z_{i1} \; \cdots\; z_{ir}] \]

Embora as suposições do termo de erro em (7-2) sejam muito modestas, mais tarde precisaremos adicionar a suposição de multinormalidade para intervalos de confiança e testar hipóteses.

Agora, forneceremos alguns exemplos do modelo de regressão linear.

Exemplo 7.1: Ajustando um modelo de regressão linear para uma reta

Determine o modelo de regressão linear para ajustar uma reta

\[ \text{Resposta média} = \mathbb{E}(Y) = \beta_0 + \beta_1 z_1 \]

para os dados

\(z_1\)	0	1	2	3	4
\(y\)	1	4	3	8	9

Antes das respostas serem observadas, os erros são estocásticos:

\[ \mathbf{Y} = \mathbf{z}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

em que

\[ \boldsymbol{\varepsilon}^{\prime} = [\varepsilon_1 \; \varepsilon_2\; \cdots \; \varepsilon_5] \]

A equação para da \(i\)-ésimo observação é

\[ Y_i = \beta_0 + \beta_1 z_{i1} + \varepsilon_i \]

em que

\[ \mathbf{Y} = \begin{bmatrix} Y_1 \\ Y_2 \\ Y_3 \\ Y_4 \\ Y_5 \end{bmatrix} \]

é o vetor de respostas,

\[ \mathbf{z} = \begin{bmatrix} 1 & z_{11} \\ 1 & z_{21} \\ 1 & z_{31} \\ 1 & z_{41} \\ 1 & z_{51} \end{bmatrix} \]

é a matriz de design,

\[ \boldsymbol{\varepsilon} = \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \\ \varepsilon_4 \\ \varepsilon_5 \end{bmatrix} \]

é o vetor de erros aleatórios, e

\[ \boldsymbol{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \end{bmatrix} \]

é o vetor de parâmetros desconhecidos.

Os dados para esse modelo de regressão linear simples estão contidos no vetor de resposta observada \(\mathbf{y}\) e na matriz de design \(\mathbf{z}\), em que

\[ \mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \end{bmatrix} = \begin{bmatrix} 1 \\ 4 \\ 3 \\ 8 \\ 9 \end{bmatrix} \]

e

\[ \mathbf{z} = \begin{bmatrix} 1 & z_{11} \\ 1 & z_{21} \\ 1 & z_{31} \\ 1 & z_{41} \\ 1 & z_{51} \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 1 & 1 \\ 1 & 2 \\ 1 & 3 \\ 1 & 4 \end{bmatrix} \]

Note a seguir que podemos lidar com uma expressão quadrática para a resposta média introduzindo o termo \(\beta_2z_2\) com \(z_2 = z_{1}^2\). Os dados para este modelo de regressão linear múltipla estão contidos no vetor de resposta observável \(\mathbf{Y}\) e na matriz de design \(\mathbf{Z}\), em que

\[ \begin{aligned} \mathbf{Y} & = \begin{bmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{bmatrix} \\ \mathbf{z} & = \begin{bmatrix} 1 & z_{11} & z_{11}^2 \\ 1 & z_{21} & z_{21}^2 \\ \vdots & \vdots & \vdots \\ 1 & z_{n1} & z_{n1}^2 \end{bmatrix} \end{aligned} \]

O modelo de regressão linear múltipla para a \(i\)-ésimo observação neste último caso pode ser:

Modelo com duas variáveis preditoras:

\[ Y_i = \beta_0 + \beta_1 z_{i1} + \beta_2 z_{i2} + \varepsilon_i \]

ou

Modelo com uma variável preditora com efeitos linear e quadrático (polinômio de ordem 2):

\[ Y_i = \beta_0 + \beta_1 z_{i1} + \beta_2 z_{i1}^2 + \varepsilon_i \]

\[\Diamond\]

Exemplo 7.2: A matriz de delineamento para ANOVA univariada como um modelo de regressão

Determine a matriz de planejamento se o modelo de regressão linear for aplicado à situação de ANOVA unifatorial univariada independente no Exemplo 6.7.

• \(p=1, q=1, g=3\; (n_1=3, n_2=2, n_3=3)\)

Considere as seguintes amostras independentes:

População 1: 9, 6, 9

População 2: 0, 2

População 3: 3, 1, 2

Criamos as chamadas variáveis fictícias (dummy) para lidar com as três médias populacionais:

\[ \mu_1 = \mu + \tau_1\\ \mu_2 = \mu + \tau_2\\ \mu_3 = \mu + \tau_3 \]

Definimos

\[ z_{j1} = \begin{cases} 1 & \text{se a observação é da população 1} \\ 0 & \text{caso contrário} \end{cases} \]

\[ z_{j2} = \begin{cases} 1 & \text{se a observação é da população 2} \\ 0 & \text{caso contrário} \end{cases} \]

\[ z_{j3} = \begin{cases} 1 & \text{se a observação é da população 3} \\ 0 & \text{caso contrário} \end{cases} \]

e \(\beta_0 = \mu\), \(\beta_1 = \tau_1\), \(\beta_2 = \tau_2\), \(\beta_3 = \tau_3\). Então

\[ Y_j = \beta_0 + \beta_1z_{j1} + \beta_2z_{j2} + \beta_3z_{j3} + \epsilon_j\\ j = 1,2,\ldots,8 \]

em que organizamos as observações das três populações em sequência. Assim, obtemos o vetor de resposta observado e a matriz de planejamento:

\[ \begin{aligned} \underset{8 \times 1}{\mathbf{y}} &= \begin{bmatrix} 9 \\ 6 \\ 9 \\ 0 \\ 2 \\ 3 \\ 1 \\ -2 \end{bmatrix}\\ \underset{8 \times 4}{\mathbf{z}} &= \begin{bmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \end{bmatrix} \end{aligned} \]

• \(p=1, q=1, r=3\)

• JW6_Example7.2.R

invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
n1 <- 3
n2 <- 2
n3 <- 3
n <- n1 + n2 + n3
Grupo <- factor(c(rep("1", n1), rep("2", n2), rep("3", n3)))
y <- c(9, 6, 9, 0, 2, 3, 1, 2)
x <- data.frame(Grupo, y)
print(x)

  Grupo y
1     1 9
2     1 6
3     1 9
4     2 0
5     2 2
6     3 3
7     3 1
8     3 2

# One-way ANOVA (offset from reference group): 
# https://en.wikipedia.org/wiki/Design_matrix
X <- model.matrix(~Grupo, 
                  data=x) # ANOVA
print(X)

  (Intercept) Grupo2 Grupo3
1           1      0      0
2           1      0      0
3           1      0      0
4           1      1      0
5           1      1      0
6           1      0      1
7           1      0      1
8           1      0      1
attr(,"assign")
[1] 0 1 1
attr(,"contrasts")
attr(,"contrasts")$Grupo
[1] "contr.treatment"

# Multiple regression: 
# https://en.wikipedia.org/wiki/Design_matrix
Um <- model.matrix(~ 1, data = x)
print(Um)

  (Intercept)
1           1
2           1
3           1
4           1
5           1
6           1
7           1
8           1
attr(,"assign")
[1] 0

z.Grupo <- model.matrix(~ -1 + Grupo, data = x)
print(z.Grupo)

  Grupo1 Grupo2 Grupo3
1      1      0      0
2      1      0      0
3      1      0      0
4      0      1      0
5      0      1      0
6      0      0      1
7      0      0      1
8      0      0      1
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$Grupo
[1] "contr.treatment"

z <- cbind(Um, z.Grupo) # Regressão
print(z)

  (Intercept) Grupo1 Grupo2 Grupo3
1           1      1      0      0
2           1      1      0      0
3           1      1      0      0
4           1      0      1      0
5           1      0      1      0
6           1      0      0      1
7           1      0      0      1
8           1      0      0      1

out <- t(z) %*% z
print(out)

            (Intercept) Grupo1 Grupo2 Grupo3
(Intercept)           8      3      2      3
Grupo1                3      3      0      0
Grupo2                2      0      2      0
Grupo3                3      0      0      3

\[\Diamond\]

A construção de variáveis fictícias, como no Exemplo 7.2, permite que toda a análise de variância seja tratada dentro do arcabouço de regressão linear múltipla.

Análise de Regressão Multivariada por Mínimos Quadrados (MOLS)

Um dos objetivos da análise de regressão é desenvolver uma equação que permita ao pesquisador prever a resposta para valores dados das variáveis preditoras. Portanto, é necessário “ajustar” o modelo em (7-3) ao \(y_j\) observado correspondente aos valores conhecidos \(1, z_{j1}, \ldots, z_{jr}\). Ou seja, devemos determinar os valores para os coeficientes de regressão \(\boldsymbol{\beta}\) e a variância do erro \(\sigma^2\) consistentes com os dados disponíveis.

Seja \(\mathbf{b}\) os valores de tentativa para \(\boldsymbol{\beta}\). Considere a diferença \(y_j - b_0 - b_1z_{j1} - \cdots - b_rz_{jr}\) entre a resposta observada \(y_j\) e o valor \(b_0 + b_1z_{j1} + \cdots + b_rz_{jr}\) que seria esperado se \(\mathbf{b}\) fosse o vetor de parâmetro “verdadeiro”. Tipicamente, as diferenças \(y_j - b_0 - b_1z_{j1} - \cdots - b_rz_{jr}\) não serão nulas, porque a resposta flutua (de maneira caracterizada pelas suposições do termo de erro) em torno do seu valor esperado.

O método dos mínimos quadrados seleciona \(\mathbf{b}\) de modo a minimizar a soma dos quadrados dessas diferenças:

\[ \begin{align} \text{S}(\mathbf{b})&=\sum_{j=1}^{n} \left(y_j - \sum_{i=0}^{r}{b_i z_{ji}}\right)^2\\ \text{S}(\mathbf{b})&=(\mathbf{y} - \mathbf{z}\mathbf{b})^{\prime}(\mathbf{y} - \mathbf{z}\mathbf{b}) \end{align} \tag{7-4} \]

Os coeficientes \(\mathbf{b}\) escolhidos pelo critério dos mínimos quadrados são chamados de estimativas de mínimos quadrados dos parâmetros de regressão \(\boldsymbol{\beta}\). Daqui em diante, serão denotados por \(\hat{\boldsymbol{\beta}}\) para enfatizar seu papel como estimativas de \(\boldsymbol{\beta}\).

Os coeficientes \(\hat{\boldsymbol{\beta}}\) são consistentes com os dados no sentido de que produzem respostas médias estimadas (ajustadas), \(\hat{\beta}_0 + \hat{\beta}_1 z_{j1} + \cdots + \hat{\beta}_r z_{jr}\), cuja soma dos quadrados das diferenças em relação às observações \(y_j\) é a menor possível. As diferenças

\[ \hat{\varepsilon}_j = y_j - \sum_{i=0}^{r}{\hat{\beta}_i z_{ji}}\\ j = 1, 2, \ldots,n \tag{7-5} \]

são chamadas de resíduos. O vetor de resíduos \(\hat{\boldsymbol{\varepsilon}} = \mathbf{y} - \mathbf{z} \hat{\boldsymbol{\beta}}\) contém informações sobre o parâmetro desconhecido \(\sigma^2\). (Veja o Resultado 7.2.)

Resultado 7.1. Seja \(\mathbf{z}\) uma matriz com posto completo \(r + 1\). A estimativa de mínimos quadrados de \(\boldsymbol{\beta}\) em (7-3) é dada por

\[ \hat{\boldsymbol{\beta}} = (\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}\mathbf{Y} \]

Seja \(\hat{\mathbf{Y}} = \mathbf{z}\hat{\boldsymbol{\beta}} = \mathbf{H}\mathbf{Y}\) os valores estimados de \(\mathbf{Y}\), em que \(\mathbf{H} = \mathbf{z}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}\) é chamada de matriz “chapéu” (Hat). Então, os resíduos

\[ \begin{align} \hat{\boldsymbol{\varepsilon}} &= \mathbf{Y} - \hat{\mathbf{Y}} \\ &= \left[\mathbf{I} - \mathbf{z}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}\right]\mathbf{Y} \\ \hat{\boldsymbol{\varepsilon}} &= (\mathbf{I} - \mathbf{H})\mathbf{Y} \end{align} \]

satisfazem \(\mathbf{z}^{\prime}\hat{\boldsymbol{\varepsilon}} = \mathbf{0}\) e \(\mathbf{Y}^{\prime}\hat{\boldsymbol{\varepsilon}} = \mathbf{0}\). Além disso, a soma dos quadrados dos resíduos é

\[ \begin{align} \sum_{j=1}^{n}{\hat{\varepsilon}_j^2} &= \hat{\boldsymbol{\varepsilon}}^{\prime}\hat{\boldsymbol{\varepsilon}} \\ &=\mathbf{Y}^{\prime}[\mathbf{I} - \mathbf{z}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}]\mathbf{Y} \\ \sum_{j=1}^{n}{\hat{\varepsilon}_j^2} &= \mathbf{Y}^{\prime}\mathbf{Y} - \mathbf{Y}^{\prime}\mathbf{z}\hat{\boldsymbol{\beta}} \end{align} \]

Se \(\mathbf{z}\) não tem posto completo, \((\mathbf{z}^{\prime}\mathbf{z})^{-1}\) é substituído por \((\mathbf{z}^{\prime}\mathbf{z})^{-}\), uma inversa generalizada de \(\mathbf{z}\mathbf{z}\). (Veja o Exercício 7.6.)

\[\Diamond\]

Resultado 7.1 mostra como as estimativas de mínimos quadrados \(\boldsymbol{\beta}\) e os resíduos \(\boldsymbol{\varepsilon}\) podem ser obtidos a partir da matriz de planejamento \(\mathbf{Z}\) e das respostas \(\mathbf{Y}\) por operações simples de matrizes.

Exemplo 7.3: Calculando as estimativas de mínimos quadrados, os resíduos e a soma dos quadrados dos resíduos

Calcule as estimativas de mínimos quadrados \(\boldsymbol{\beta}\), os resíduos \(\boldsymbol{\varepsilon}\) e a soma dos quadrados dos resíduos para um modelo de reta

\[ Y_j=\beta_0 + \beta_1 z_{j1}+\varepsilon_j \]

ajustado para os dados

\[ \begin{array}{cc} z_1 & y \\ 0 & 1 \\ 1 & 4 \\ 2 & 3 \\ 3 & 8 \\ 4 & 9 \\ \end{array} \]

Sejam \(\mathbf{z}\) a matriz de planejamento e \(\mathbf{y}\) o vetor de respostas:

\[ \mathbf{z} = \begin{bmatrix} 1 & 0 \\ 1 & 1 \\ 1 & 2 \\ 1 & 3 \\ 1 & 4 \end{bmatrix}, \quad \mathbf{y} = \begin{bmatrix} 1 \\ 4 \\ 3 \\ 8 \\ 9 \end{bmatrix} \]

A estimativa de mínimos quadrados para \(\boldsymbol{\beta}\) é dada por:

\[ \hat{\boldsymbol{\beta}} = (\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}\mathbf{y} \]

Os resíduos são calculados como:

\[ \hat{\boldsymbol{\varepsilon}} = \mathbf{y} - \mathbf{z}\hat{\boldsymbol{\beta}} \]

E a soma dos quadrados dos resíduos é:

\[ \hat{\boldsymbol{\varepsilon}}^{\prime}\hat{\boldsymbol{\varepsilon}} \]

• \(p=1, q=1, r=1\)

• JW6_Example7.3.R

# Dados
z <- matrix(c(1, 0, 
              1, 1,
              1, 2, 
              1, 3, 
              1, 4), 
            ncol=2, byrow=TRUE)
z

     [,1] [,2]
[1,]    1    0
[2,]    1    1
[3,]    1    2
[4,]    1    3
[5,]    1    4

y <- c(1, 4, 3, 8, 9)
y

[1] 1 4 3 8 9

# Estimativas de mínimos quadrados para beta
beta_hat <- solve(t(z) %*% z) %*% t(z) %*% y
beta_hat

     [,1]
[1,]    1
[2,]    2

# Resíduos
y_hat <- z %*% beta_hat
y_hat

     [,1]
[1,]    1
[2,]    3
[3,]    5
[4,]    7
[5,]    9

residuals <- y - y_hat
round(residuals,2)

     [,1]
[1,]    0
[2,]    1
[3,]   -2
[4,]    1
[5,]    0

round(sum(residuals), 6)

[1] 0

# Soma dos quadrados dos resíduos
sqr <- sum(residuals^2)
sqr

[1] 6

plot(z[, 2], y, pch=16, col='gray', 
     main='Gráfico de Regressão', xlab='z', ylab='y')
curve(beta_hat[1]+beta_hat[2]*x, min(z), max(z), col='black', add=TRUE)
legend('topleft', legend=c('Dados', 'Reta de regressão'), 
       col=c('gray', 'black'), pch=c(16, NA), lty=c(NA, 1), bty="n")

Equação de regressão linear simples univariada por OLS:

\[ \hat{y} = 1 + 2z_1\\ z_1\in[0,4] \]

Este segmento de reta representa a média de \(y\) para cada valor de \(z_1\in[0,4]\). Um dos valores pelos quais a reta passa é o centróide \((\bar{y}, \bar{z}_1)\), sendo que \(\bar{y}=\frac{1}{5}(1+4+3+8+9)=5\) e \(\bar{z}_1=\frac{1}{5}(0+1+2+3+4)=2\):

\[ \begin{align} \hat{y}(\bar{z}_1)&=\bar{y}\\ \hat{y}(2)&=1+2\times 2=5 \end{align} \]

Visualização de OLS

Infinitas retas passam pelo centróide. Construindo quadrados cujos lados sejam a distância entre os valores observados \(y\) e valores os estimados \(\hat{y}\) por uma reta, a que melhor se ajusta é aquela cuja somatória das áreas dos quadrados for mínima. Esta reta ajustada pelo método de mínimos quadrados ordinários é a que foi representada no gráfico (toda reta ajustada por este método passa obrigatoriamente pelo centróide).

Execute demo_minimos_quadrados.R para uma demonstração:

source("demo_minimos_quadrados.R")

Método dos mínimos quadrados em WolframAlpha, R e SCILAB

LeastSquares[{{1,0}, {1,1}, {1,2}, {1,3}, {1,4}}, {1,4,3,8,9}]

• JW6_Example7.3_R.R

Y <- c(1,4,3,8,9)
X <- matrix(c(1, 0, 
              1, 1,
              1, 2, 
              1, 3, 
              1, 4), 
            ncol=2, byrow=TRUE)
Y

[1] 1 4 3 8 9

X

     [,1] [,2]
[1,]    1    0
[2,]    1    1
[3,]    1    2
[4,]    1    3
[5,]    1    4

b <- solve(t(X)%*%X)%*%t(X)%*%Y
b

     [,1]
[1,]    1
[2,]    2

solve(crossprod(X), crossprod(X,Y))

     [,1]
[1,]    1
[2,]    2

y <- X%*%b
e <- Y - y
P <- X%*%solve(t(X)%*%X)%*%t(X)
matrixcalc::is.symmetric.matrix(P)

[1] TRUE

matrixcalc::is.idempotent.matrix(P)

[1] TRUE

Q <- t(Y)%*%e
Q

     [,1]
[1,]    6

Q <- t(e)%*%e
Q

     [,1]
[1,]    6

e.soma <- round(t(e)%*%rep(1,length(e)), 6)
e.soma

     [,1]
[1,]    0

ortg <-  round(t(y)%*%e, 6)
ortg

     [,1]
[1,]    0

n <- length(Y)
R2 <- (t(y)%*%Y/n-mean(Y)^2)/(t(Y)%*%Y/n-mean(Y)^2)
R2

          [,1]
[1,] 0.8695652

plot(X[,2], Y,
     main=paste0("OLS: y = Xb\nR^2 = ", round(R2,2)),
     ylab="VD",
     xlab="X",
     ylim=1.1*c(min(y), max(y)))
curve(b[1]+b[2]*x, min(X[,2]), max(X[,2]), add=TRUE)

Y=[1;4;3;8;9]
X=[1 0;1 1;1 2;1 3;1 4]
b=inv(X'*X)*X'*Y
lsq(X,Y)
y=X*b
e=Y-y
P=X*inv(X'*X)*X'
Q=Y'*e
Q=e'*e
esum=e'*[1 1 1 1 1]'
ortg=y'*e
n=length(Y)
R2=(y'*Y/n-mean(Y)^2)/(Y'*Y/n-mean(Y)^2)
quit

Scilab 2023.1.0 (May 23 2023, 09:23:00)
 Y  = 

   1.
   4.
   3.
   8.
   9.
 X  = 

   1.   0.
   1.   1.
   1.   2.
   1.   3.
   1.   4.
 b  = 

   1.0000000
   2.0000000
 ans  =

   1.0000000
   2.0000000
 y  = 

   1.0000000
   3.
   5.0000000
   7.0000000
   9.0000000
 e  = 

   3.331D-16
   1.
  -2.0000000
   1.0000000
  -1.776D-15
 P  = 

   0.6   0.4         0.2  -2.776D-17  -0.2      
   0.4   0.3         0.2   0.1        -5.551D-17
   0.2   0.2         0.2   0.2         0.2      
   0.    0.1         0.2   0.3         0.4      
  -0.2   2.776D-17   0.2   0.4         0.6      
 Q  = 

   6.0000000
 Q  = 

   6.0000000
 esum  = 

  -3.220D-15
 ortg  = 

  -2.798D-14
 n  = 

   5.
 R2  = 

   0.8695652

\[\Diamond\]

Demonstração do método de mínimos quadrados

• \(p=1, q=1, r=2\)

• Proofs involving ordinary least squares: Wikipedia

A álgebra matricial é indispensável em áreas como análise estatística multivariada, planejamento de experimentos e análise de variância e covariância.

Para dar uma ideia de tal aplicação estatística, consideramos um modelo de regressão. Suponha que nos seja dado um diagrama de dispersão consistindo de pontos \((z_i,y_i)\), \(i = 1,2, \ldots, n\), em que \(x_i\) e \(y_i\) são observações intervalares. Suponha ainda que os pontos estejam próximos de uma certa curva e que queremos ajustar uma função quadrática. a equação é:

\[ y=\beta_0+\beta_1z+\beta_2z^2 \]

em que as constantes \(\beta_0\), \(\beta_1\) e \(\beta_2\) são desconhecidas. As coordenadas \(z_i\) e \(y_i\) dos pontos não satisfazem exatamente uma equação quadrática. Portanto, escrevemos:

\[ Y_i=\beta_0+\beta_1z_i+\beta_2z_i^2+\varepsilon_i \]

em que \(\epsilon_i\) é o termo de erro.

O sistema com \(n\) equações lineares nos parâmetros pode ser escrito na forma matricial:

O sistema com \(n\) equações lineares nos parâmetros pode ser escrito na forma matricial:

\[ \mathbf{Y} = \mathbf{z}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

Em que:

\[ \mathbf{Y} = \begin{bmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{bmatrix} \]

\[ \mathbf{z} = \begin{bmatrix} 1 & z_1 & z_1^2 \\ 1 & z_2 & z_2^2 \\ \vdots & \vdots & \vdots \\ 1 & z_n & z_n^2 \end{bmatrix} \]

\[ \boldsymbol{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{bmatrix} \]

\[ \boldsymbol{\varepsilon} = \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{bmatrix} \]

A equação matricial \(\mathbf{Y} = \mathbf{z}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\) é chamada de modelo linear, pois é linear em \(\boldsymbol{\beta}\).

Para encontrar coeficientes adequados, aplicamos o método dos mínimos quadrados (OLS), que minimiza a soma dos quadrados dos erros:

\[ \sum_{i=1}^{n} \varepsilon_i^2 = \boldsymbol{\varepsilon}^\prime \boldsymbol{\varepsilon} \]

Isso significa que os coeficientes desconhecidos \(\beta_0\), \(\beta_1\) e \(\beta_2\) são determinados de tal forma que:

\[ \begin{align} Q &= \boldsymbol{\varepsilon}^\prime \boldsymbol{\varepsilon}\\ Q &= (\mathbf{Y} - \mathbf{z}\boldsymbol{\beta})^\prime (\mathbf{Y} - \mathbf{z}\boldsymbol{\beta}) \end{align} \]

assume o menor valor possível. A solução pode ser encontrada por métodos de álgebra matricial ou podemos aplicar cálculo diferencial. Conforme Teichroew (1964, p. 580-3) e Tay (2018):

\[ \begin{align} Q &= \boldsymbol{\varepsilon}^\prime \boldsymbol{\varepsilon} \\ &= (\mathbf{Y} - \mathbf{z}\boldsymbol{\beta})^\prime (\mathbf{Y} - \mathbf{z}\boldsymbol{\beta}) \\ Q &= \boldsymbol{\beta}^\prime \mathbf{z}^\prime \mathbf{z}\boldsymbol{\beta} - 2\mathbf{Y}^\prime \mathbf{z}\boldsymbol{\beta} + \mathbf{Y}^\prime \mathbf{Y} \end{align} \]

A condição de primeira ordem (necessária) para encontrar o mínimo de \(Q\) é:

\[ \begin{aligned} \dfrac{\partial Q}{\partial \boldsymbol{\beta}} &= 2\mathbf{z}^\prime \mathbf{z}\boldsymbol{\beta} - 2\mathbf{z}^\prime \mathbf{Y} \\ \mathbf{0} &= 2\mathbf{z}^\prime \mathbf{z}\hat{\boldsymbol{\beta}} - 2\mathbf{z}^\prime \mathbf{Y} \\ \hat{\boldsymbol{\beta}} &= (\mathbf{z}^\prime \mathbf{z})^{-1} \mathbf{z}^\prime \mathbf{Y} \end{aligned} \]

A condição de segunda ordem (suficiente) para um mínimo é:

\[ \left.\dfrac{\partial^2 Q}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^\prime}\right|_{\boldsymbol{\beta}=\hat{\boldsymbol{\beta}}} = 2\mathbf{z}^\prime \mathbf{z} \]

A matriz de informação \(\mathbf{z}^\prime \mathbf{z}\) é uma matriz simétrica \(3 \times 3\) e definida positiva (todos os autovalores são positivos), ou seja, \(\mathbf{z}^\prime \mathbf{z} > 0\). Se a segunda derivada é positiva, então o ponto crítico é de mínimo. Portanto, \(\hat{\boldsymbol{\beta}}\) é a estimativa de mínimos quadrados:

\[ \begin{align} Q_{\min} &= \mathbf{Y}^\prime (\mathbf{Y} - \mathbf{z}\hat{\boldsymbol{\beta}}) \\ &= \mathbf{Y}^\prime \hat{\boldsymbol{\varepsilon}} \\ Q_{\min} &= \hat{\boldsymbol{\varepsilon}}^\prime \hat{\boldsymbol{\varepsilon}} \end{align} \]

Em que \(\hat{\boldsymbol{\varepsilon}} = \mathbf{Y} - \mathbf{z}\hat{\boldsymbol{\beta}}\) é o vetor de resíduos. As relações adicionais incluem:

\[ \begin{aligned} \mathbf{\hat{Y}} &= \mathbf{z}\hat{\boldsymbol{\beta}} \\ \mathbf{Y} - \hat{\boldsymbol{\varepsilon}} &= \mathbf{z}\hat{\boldsymbol{\beta}} \\ \mathbf{Y} &= \mathbf{z}\hat{\boldsymbol{\beta}} + \hat{\boldsymbol{\varepsilon}} \end{aligned} \]

Além disso, a soma dos quadrados dos resíduos é nula:

\[ \hat{\boldsymbol{\varepsilon}}^\prime \mathbf{1} = 0 \]

A matriz de planejamento \(\mathbf{z}\) e o vetor das estimativas de resíduos são ortogonais (\(\mathbf{z}\) e \(\hat{\boldsymbol{\varepsilon}}\) são matematicamente independentes):

\[ \begin{aligned} \mathbf{Y} &= \mathbf{z}\hat{\boldsymbol{\beta}} \\ \mathbf{Y} - \mathbf{z}\hat{\boldsymbol{\beta}} &= 0 \\ \mathbf{z}^\prime \mathbf{Y} - \mathbf{z}^\prime \mathbf{z}\hat{\boldsymbol{\beta}} &= 0 \\ \mathbf{z}^\prime (\mathbf{Y} - \mathbf{z}\hat{\boldsymbol{\beta}}) &= 0 \\ \mathbf{z}^\prime \hat{\boldsymbol{\varepsilon}} &= 0 \end{aligned} \]

O estimador OLS de \(\beta\) é não-viesado se \(\mathbb{E}(\boldsymbol{\varepsilon}) = 0\):

\[ \begin{aligned} \hat{\boldsymbol{\beta}} &= (\mathbf{z}^\prime \mathbf{z})^{-1} \mathbf{z}^\prime \mathbf{Y} \\ &= (\mathbf{z}^\prime \mathbf{z})^{-1} \mathbf{z}^\prime (\mathbf{z}\boldsymbol{\beta} + \boldsymbol{\varepsilon}) \\ &= (\mathbf{z}^\prime \mathbf{z})^{-1} \mathbf{z}^\prime \mathbf{z}\boldsymbol{\beta} + (\mathbf{z}^\prime \mathbf{z})^{-1} \mathbf{z}^\prime \boldsymbol{\varepsilon} \\ \hat{\boldsymbol{\beta}} &= \boldsymbol{\beta} + (\mathbf{z}^\prime \mathbf{z})^{-1} \mathbf{z}^\prime \boldsymbol{\varepsilon} \\ \mathbb{E}(\hat{\boldsymbol{\beta}}) &= \boldsymbol{\beta} + (\mathbf{z}^\prime \mathbf{z})^{-1} \mathbf{z}^\prime \mathbb{E}(\boldsymbol{\varepsilon}) \\ \mathbb{E}(\hat{\boldsymbol{\beta}}) &= \boldsymbol{\beta} \end{aligned} \]

Os valores preditos de \(\mathbf{Y}\) são:

\[ \begin{align} \mathbf{\hat{Y}} &= \mathbf{z}\hat{\boldsymbol{\beta}} \\ &= \mathbf{z}(\mathbf{z}^\prime \mathbf{z})^{-1} \mathbf{z}^\prime \mathbf{Y} \\ \mathbf{\hat{Y}} &= \mathbf{H}\mathbf{Y} \end{align} \]

Em que \(\mathbf{H} = \mathbf{z}(\mathbf{z}^\prime \mathbf{z})^{-1} \mathbf{z}^\prime\) é chamada de matriz de projeção (Projection matrix), matriz de influência (influence matrix) ou matriz chapéu (Hat matrix): https://en.wikipedia.org/wiki/Projection_matrix.

\(\mathbf{H}\) é simétrica (\(\mathbf{H}^\prime = \mathbf{H}\)) e idempotente (\(\mathbf{H}^2 = \mathbf{H}\)).

O vetor da variável dependente predita \(\mathbf{\hat{Y}}\) e o vetor das estimativas de resíduos são ortogonais (\(\mathbf{\hat{Y}}\) e \(\hat{\boldsymbol{\varepsilon}}\) são matematicamente independentes):

\[ \mathbf{\hat{Y}}^\prime \hat{\boldsymbol{\varepsilon}} = 0 \]

Decomposição da Soma dos Quadrados

De acordo com o Resultado 7.1, \(\hat{\mathbf{y}}^{\prime}\hat{\boldsymbol{\varepsilon}} = 0\), então a soma total dos quadrados da resposta \(\mathbf{y}^{\prime}\mathbf{y}=\sum_{j=1}^{n}{y^2_j}\) satisfaz:

\[ \begin{align} \mathbf{y}^{\prime}\mathbf{y} &= (\hat{\mathbf{y}} + \mathbf{y} - \hat{\mathbf{y}})^{\prime}(\hat{\mathbf{y}} + \mathbf{y} - \hat{\mathbf{y}})\\ &= (\hat{\mathbf{y}} + \hat{\boldsymbol{\varepsilon}})^{\prime}(\hat{\mathbf{y}} + \hat{\boldsymbol{\varepsilon}})\\ \mathbf{y}^{\prime}\mathbf{y}&= \hat{\mathbf{y}}^{\prime}\hat{\mathbf{y}} + \hat{\boldsymbol{\varepsilon}}^{\prime}\hat{\boldsymbol{\varepsilon}} \end{align} \tag{7-7} \]

Uma vez que a primeira coluna de \(\mathbf{z}\) é \(\mathbf{1}\), a condição \(\mathbf{z}^{\prime}\hat{\boldsymbol{\varepsilon}} = \mathbf{0}\) inclui o requisito

\[ \begin{align} 0 &= \mathbf{1}^{\prime}\hat{\boldsymbol{\varepsilon}} \\ &= \sum_{j=1}^n \hat{\varepsilon}_j \\ 0 &= \sum_{j=1}^n (y_j - \hat{y}_j) \end{align} \]

ou

\[ \bar{y} = \bar{\hat{y}} \]

Subtraindo \(n\bar{y}^2 = n\left(\bar{\hat{y}}\right)^2\) de ambos os lados da decomposição em (7-7), obtemos a decomposição básica da soma dos quadrados em relação à média:

\[ \mathbf{y}^{\prime}\mathbf{y} - n\bar{y}^2 = \hat{\mathbf{y}}^{\prime}\hat{\mathbf{y}} - n\left(\bar{\hat{y}}\right)^2 + \hat{\boldsymbol{\varepsilon}}^{\prime}\hat{\boldsymbol{\varepsilon}} \]

ou

\[ \sum_{j=1}^n (y_j - \bar{y})^2 = \sum_{j=1}^n (\hat{y}_j - \bar{y})^2 + \sum_{j=1}^n \hat{\varepsilon}_j^2 \tag{7-8} \]

\[\text{SS em relação à média} = \text{SS da regressão} + \text{SS dos resíduos}\]

A decomposição da soma dos quadrados precedente sugere que a qualidade do ajuste do modelo pode ser medida pelo coeficiente de determinação \(R^2\):

\[ R^2 = \dfrac{\sum_{j=1}^n (\hat{y}_j - \bar{y})^2}{\sum_{j=1}^n (y_j - \bar{y})^2} \tag{7-9} \]

A quantidade \(R^2\) fornece a proporção da variação total de \(y\) “explicada” por, ou atribuível a, as variáveis preditoras \(z_1, z_2, \ldots, z_r\). Aqui \(R^2\) (ou o coeficiente de correlação múltipla \(R = +\sqrt{R^2}\)) é igual a 1 se a equação ajustada passar por todos os pontos observados, de modo que \(\hat{\varepsilon}_j = 0\) para todos \(j\). No outro extremo, \(R^2\) é 0 se \(\hat{\beta}_0 = \bar{y}\) e \(\hat{\beta}_1 = \hat{\beta}_2 = \cdots = \hat{\beta}_r = 0\). Nesse caso, as variáveis preditoras \(z_1, z_2, \ldots, z_r\) não têm influência na resposta.

O \(R^2\) é uma medida estatística de quão próximos os dados estão da função de regressão ajustada. Esta medida varia entre 0 (ausência de ajuste) e 1 (ajuste perfeito). \(R^2\) é conhecido como o coeficiente de determinação ou o coeficiente de determinação múltipla para a regressão múltipla.

\[ R^2 = \dfrac{\dfrac{\mathbf{\hat{Y}}^\prime \mathbf{Y}}{n} - \bar{Y}^2}{\dfrac{\mathbf{Y}^\prime \mathbf{Y}}{n} - \bar{Y}^2} \]

Geometria dos Mínimos Quadrados

Uma interpretação geométrica da técnica de mínimos quadrados destaca a natureza do conceito. De acordo com o modelo clássico de regressão linear,

Vetor resposta média é:

\[ \begin{align} \mathbb{E}(\mathbf{Y})&=\mathbf{z}\boldsymbol{\beta}\\ \mathbb{E}(\mathbf{Y})&=\beta_0\begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix} + \beta_1 \begin{bmatrix} z_{11} \\ z_{21} \\ \vdots \\ z_{n1} \end{bmatrix} + \cdots + \beta_r \begin{bmatrix} z_{1r} \\ z_{2r} \\ \vdots \\ z_{nr} \end{bmatrix} \end{align} \]

Assim, \(\mathbb{E}(\mathbf{Y})\) é uma combinação linear das colunas de \(\mathbf{z}\). À medida que \(\boldsymbol{\beta}\) varia, \(\mathbf{z}\boldsymbol{\beta}\) abrange o plano do modelo de todas as combinações lineares. Geralmente, o vetor de observação \(\mathbf{Y}\) não estará no plano do modelo, devido ao erro aleatório \(\boldsymbol{\varepsilon}\); ou seja, \(\mathbf{Y}\) não é (exatamente) uma combinação linear das colunas de \(\mathbf{z}\). Lembre-se de que

\[ \mathbf{Y} = \begin{array}{c} \mathbf{z}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \end{array} \]

Figura 7.1 Mínimos quadrados com uma projeção para n = 3 e r = 1.

A estimativa OLS pode ser vista como uma projeção no espaço linear estendido pelos regressores. X₁ e X₂ referem-se às colunas da matriz de X.) Wikipedia: https://en.wikipedia.org/wiki/Ordinary_least_squares

• JW6_Cap7_RLS3D.R

invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
vec <- rbind(diag(3), c(1,1,1))
rownames(vec) <- c("1", "2", "3", "y")
rgl::open3d()

wgl 
  1 

matlib::vectors3d(vec, color=c(rep("black",3), "red"), lwd=2)
# draw the XZ plane, whose equation is Y=0
rgl::planes3d(0, 0, 1, 0, col="gray", alpha=0.2)
matlib::vectors3d(c(1,1,0), col="green", lwd=2)
# show projections of the unit vector J
rgl::segments3d(rbind(c(1,1,1), c(1, 1, 0)))
rgl::segments3d(rbind(c(0,0,0), c(1, 1, 0)))
rgl::segments3d(rbind(c(1,0,0), c(1, 1, 0)))
rgl::segments3d(rbind(c(0,1,0), c(1, 1, 0)))
# show some orthogonal vectors
p1 <- c(0,0,0)
p2 <- c(1,1,0)
p3 <- c(1,1,1)
p4 <- c(1,0,0)
p5 <- c(0,1,0)
matlib::corner(p1, p2, p3, col="red")
matlib::corner(p1, p4, p2, col="red")
matlib::corner(p1, p5, p2, col="blue")
rgl::rglwidget()

Uma vez que as observações estejam disponíveis, a solução de mínimos quadrados é derivada do vetor de desvio

\[ \mathbf{y} - \mathbf{z}\mathbf{b} = (\text{vetor de observação}) - (\text{vetor no plano do modelo}) \]

O comprimento ao quadrado \((\mathbf{y} - \mathbf{z}\mathbf{b})^{\prime}(\mathbf{y} - \mathbf{z}\mathbf{b})\) é a soma dos quadrados \(S(\mathbf{b})\). Como ilustrado na Figura 7.1, \(S(\mathbf{b})\) é o menor possível quando \(\mathbf{b}\) é selecionado de modo que \(\mathbf{z}\mathbf{b}\) seja o ponto no plano do modelo mais próximo de \(\mathbf{y}\). Esse ponto ocorre na ponta da projeção perpendicular de \(\mathbf{y}\) no plano. Ou seja, para a escolha \(\mathbf{b} = \hat{\boldsymbol{\beta}}\), \(\hat{\mathbf{y}} = \mathbf{z}\hat{\boldsymbol{\beta}}\) é a projeção de \(\mathbf{y}\) no plano consistindo de todas as combinações lineares das colunas de \(\mathbf{z}\). O vetor residual \(\hat{\boldsymbol{\varepsilon}} = \mathbf{y} - \hat{\mathbf{y}}\) é perpendicular a esse plano. Essa geometria se mantém mesmo quando \(\mathbf{z}\) não é de posto completo.

Quando \(\mathbf{z}\) tem posto completo, a operação de projeção é expressa analiticamente como multiplicação pela matriz \(\mathbf{z}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}\). Para ver isso, usamos a decomposição espectral (2-16) para escrever

\[ \mathbf{z}^{\prime}\mathbf{z} = \sum_{i=1}^{r+1} \lambda_i \mathbf{e}_i \mathbf{e}_i^{\prime} \]

em que \(\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_{r+1} > 0\) são os autovalores de \(\mathbf{z}^{\prime}\mathbf{z}\) e \(\mathbf{e}_1, \mathbf{e}_2, \ldots, \mathbf{e}_{r+1}\) são os autovetores correspondentes. Se \(\mathbf{z}\) é de posto completo,

\[ (\mathbf{z}^{\prime}\mathbf{z})^{-1} = \sum_{i=1}^{r+1} \frac{1}{\lambda_i} \mathbf{e}_i \mathbf{e}_i^{\prime} \]

Considere \(\mathbf{q}_i = \lambda_i^{-1/2} \mathbf{z} \mathbf{e}_i\), que é uma combinação linear das colunas de \(\mathbf{z}\). Então \(\mathbf{q}_i^{\prime}\mathbf{q}_j = 0\). Isso significa que os \(r + 1\) vetores \(\mathbf{q}_i\) são mutuamente perpendiculares e têm comprimento unitário. Suas combinações lineares abrangem o espaço de todas as combinações lineares das colunas de \(\mathbf{z}\). Além disso,

\[ \mathbf{z}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime} = \sum_{i=1}^{r+1} \mathbf{q}_i \mathbf{q}_i^{\prime} \]

De acordo com o Resultado 2A.2 e a Definição 2A.12, a projeção de \(\mathbf{y}\) em uma combinação linear de \(\{\mathbf{q}_1, \mathbf{q}_2, \ldots, \mathbf{q}_{r+1}\}\) é

\[ \sum_{i=1}^{r+1}{(\mathbf{q}_{i}^{\prime}\mathbf{y})\mathbf{q}_{i}}=\mathbf{z}\hat{\boldsymbol{\beta}} \]

Assim, a multiplicação por \(\mathbf{z}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}\) projeta um vetor sobre o espaço abrangido pelas colunas de \(\mathbf{z}\).

Da mesma forma, \([\mathbf{I} - \mathbf{z}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}]\) é a matriz para a projeção de \(\mathbf{y}\) no plano perpendicular ao plano abrangido pelas colunas de \(\mathbf{z}\).

Se \(\mathbf{z}\) não for de posto completo, podemos usar a inversa generalizada \((\mathbf{z}^{\prime}\mathbf{z})^{-}=\sum_{i=1}^{r+1}{\lambda_i^{-1}\mathbf{e}_i\mathbf{e}_i^{\prime}}\) - um tipo especial de inversa que pode ser usada quando a matriz não é invertível no sentido usual. No caso de \(\mathbf{z}^{\prime}\mathbf{z}\) não ser de posto completo, os autovalores \(\lambda_1\ge \lambda_2\ge \cdots \ge \lambda_{r_1+1}>0=\lambda_{r_1+2} = \cdots = \lambda_{r+1} = 0\), como descrito no Exercício 7.6. Então, \(\mathbf{z}(\mathbf{z}^{\prime}\mathbf{z})^{-}\mathbf{z}^{\prime}"=\sum_{i=1}^{r+1}{\mathbf{q}_i\mathbf{q}_i^{\prime}}\) gera a projeção única de \(\mathbf{y}\) no espaço abrangido pelas colunas linearmente independentes de \(\mathbf{z}\). Isso é verdadeiro para qualquer escolha da inversa generalizada. (Veja [23].)

Propriedades Amostrais dos Estimadores de Mínimos Quadrados Clássicos

O estimador de mínimos quadrados \(\hat{\boldsymbol{\beta}}\) e os resíduos \(\hat{\boldsymbol{\varepsilon}}\) possuem as propriedades amostrais detalhadas no próximo resultado.

Resultado 7.2. Sob o modelo de regressão linear geral em (7-3), o estimador de mínimos quadrados \(\hat{\boldsymbol{\beta}} = (\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}\mathbf{Y}\) possui

\[ \mathbb{E}(\hat{\boldsymbol{\beta}}) = \boldsymbol{\beta} \]

e

\[ \mathbb{C}(\hat{\boldsymbol{\beta}}) = \sigma^2(\mathbf{z}^{\prime}\mathbf{z})^{-1} \]

Os resíduos \(\hat{\boldsymbol{\varepsilon}}\) possuem as propriedades

\[ \begin{align} \mathbb{E}(\hat{\boldsymbol{\varepsilon}}) &= \mathbf{0}\\\\ \mathbb{C}(\hat{\boldsymbol{\varepsilon}}) &= \sigma^2[\mathbf{I} - \mathbf{z}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}] \\ \mathbb{C}(\hat{\boldsymbol{\varepsilon}}) &= \sigma^2[\mathbf{I} - \mathbf{H}] \end{align} \]

A soma de quadrados dos resíduos é:

\[ \text{SS}_{\text{res}}=\hat{\boldsymbol{\varepsilon}}^{\prime}\hat{\boldsymbol{\varepsilon}} \]

Além disso,

\[ \mathbb{E}(\hat{\boldsymbol{\varepsilon}}^{\prime}\hat{\boldsymbol{\varepsilon}}) = (n - r - 1)\sigma^2 \]

portanto, definindo

\[ \begin{align} s^2 &= \dfrac{\hat{\boldsymbol{\varepsilon}}^{\prime}\hat{\boldsymbol{\varepsilon}}}{n-\text{rank}(\mathbf{z})}\\ &= \dfrac{\hat{\boldsymbol{\varepsilon}}^{\prime}\hat{\boldsymbol{\varepsilon}}}{n-(r+1)}\\ &= \dfrac{\text{SS}_{\text{res}}}{n - r - 1}\\ &=\dfrac{\mathbf{Y}^{\prime}\left(\mathbf{I} - \mathbf{z}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}\right)\mathbf{Y}}{n - r - 1} \\ s^2&= \dfrac{\mathbf{Y}^{\prime}(\mathbf{I} - \mathbf{H})\mathbf{Y}}{\text{tr}(\mathbf{I} - \mathbf{H})} \end{align} \]

temos

\[ \mathbb{E}(s^2) = \sigma^2 \]

Além disso, \(\mathbb{C}(\hat{\boldsymbol{\beta}},\hat{\boldsymbol{\varepsilon}})=\mathbf{0}\).

Demonstração: J&W6e, Solutions Manual, 2007, p. 164

\[\Diamond\]

O estimador de mínimos quadrados \(\hat{\boldsymbol{\beta}}\) possui uma propriedade de variância mínima que foi estabelecida por Gauss. O seguinte resultado diz respeito aos estimadores “melhores” de funções lineares paramétricas da forma \(\mathbf{c}^{\prime} \boldsymbol{\beta} = \sum_{i=1}^r c_i \beta_i\) para qualquer vetor \(\mathbf{c}\).

Resultado 7.3 (Teorema dos mínimos quadrados de Gauss). Seja \(\mathbf{Y} = \mathbf{z}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\), onde \(\mathbb{E}(\boldsymbol{\varepsilon}) = 0\), \(\mathbb{C}(\boldsymbol{\varepsilon}) = \sigma^2\mathbf{I}\), e \(\mathbf{z}\) tem posto completo \(r + 1\). Para qualquer vetor \(\mathbf{c}\), o estimador

\[ \mathbf{c}^{\prime}\hat{\boldsymbol{\beta}} = \sum_{i=1}^r c_i\hat{\beta}_i \]

de \(\mathbf{c}^{\prime} \boldsymbol{\beta}\) possui a menor variância possível entre todos os estimadores lineares da forma

\[ \mathbf{a}^{\prime}\mathbf{Y} = \sum_{i=1}^n a_i Y_i \]

que são não viesados para \(\mathbf{c}^{\prime}\boldsymbol{\beta}\).

\[\Diamond\]

Este resultado poderoso afirma que a substituição de \(\boldsymbol{\beta}\) por \(\hat{\boldsymbol{\beta}}\) leva ao melhor estimador de \(\mathbf{c}^{\prime}\boldsymbol{\beta}\) para qualquer \(\mathbf{c}\) de interesse. Em terminologia estatística, o estimador \(\mathbf{c}^{\prime}\hat{\boldsymbol{\beta}}\) é chamado de melhor estimador linear não-viesado de variância mínima (BLUE) de \(\mathbf{c}^{\prime}\boldsymbol{\beta}\).

Inferência sobre o Modelo de Regressão

Descrevemos procedimentos inferenciais baseados no modelo de regressão linear clássico em (7-3) com a suposição adicional (tentativa) de que os erros \(\boldsymbol{\varepsilon}\) têm uma distribuição normal. Métodos para verificar a adequação geral do modelo são considerados na Seção 7.6.

Inferências Relativas aos Parâmetros de Regressão

Antes de podermos avaliar a importância de variáveis específicas na função de regressão \[\mathbb{E}(\mathbf{Y}) = \sum_{i=0}^r \beta_i z_i \tag{7-10}\] devemos determinar as distribuições amostrais de \(\hat{\boldsymbol{\beta}}\) e a soma dos quadrados dos resíduos, \(\hat{\boldsymbol{\varepsilon}}^{\prime}\hat{\boldsymbol{\varepsilon}}\). Para fazer isso, vamos supor que os erros \(\boldsymbol{\varepsilon}\) tenham uma distribuição multinormal.

Resultado 7.4. Seja \(\mathbf{Y} = \mathbf{z}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\), em que \(\mathbf{z}\) tem posto completo \(r + 1\) e \(\boldsymbol{\varepsilon}\) é distribuído como \(\mathcal{N}_n(\mathbf{0}, \sigma^2\mathbf{I})\). Então, o estimador de máxima verossimilhança de \(\boldsymbol{\beta}\) é o mesmo que o estimador de mínimos quadrados \(\hat{\boldsymbol{\beta}}\). Além disso,

\[ \hat{\boldsymbol{\beta}} = (\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}\mathbf{Y} \sim \mathcal{N}_{r+1}(\boldsymbol{\beta}, \sigma^2(\mathbf{z}^{\prime}\mathbf{z})^{-1}) \]

\(\hat{\boldsymbol{\beta}}\) é distribuído independentemente dos resíduos \(\hat{\boldsymbol{\varepsilon}} = \mathbf{Y} - \mathbf{z}\hat{\boldsymbol{\beta}}\). Adicionalmente,

\[ \begin{align} \hat{\sigma}^2 &= \dfrac{\hat{\boldsymbol{\varepsilon}}^{\prime}\hat{\boldsymbol{\varepsilon}}}{n}\\ n\hat{\sigma}^2 &= \hat{\boldsymbol{\varepsilon}}^{\prime}\hat{\boldsymbol{\varepsilon}} \sim\chi^2_{n - r - 1} \end{align} \]

em que \(\hat{\sigma}^2\) é o estimador de máxima verossimilhança de \(\sigma^2\).

Demonstração: J&W6e, Solutions Manual, 2007, p. 165

\[\Diamond\]

Um elipsoide de confiança para \(\boldsymbol{\beta}\) é facilmente construído. Ele é expresso em termos da matriz de covariância estimada \(s^2 (\mathbf{z}^{\prime}\mathbf{z})^{-1}\), em que \(s^2 = \hat{\boldsymbol{\varepsilon}}^{\prime}\hat{\boldsymbol{\varepsilon}} / (n - r - 1)\).

Resultado 7.5. Seja \(\mathbf{Y} = \mathbf{z}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\), em que \(\mathbf{z}\) tem posto completo \(r + 1\) e \(\boldsymbol{\varepsilon} \sim \mathcal{N}_n(\mathbf{0}, \sigma^2\mathbf{I})\). Então,

1. Uma região de confiança de \(100(1 - \alpha)\%\) para \(\boldsymbol{\beta}\) é dada por

\[ (\boldsymbol{\beta}-\hat{\boldsymbol{\beta}})^{\prime}\mathbf{z}^{\prime}\mathbf{z}(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}}) \leq s^2(r + 1)F_{r + 1, \,n - r - 1}(1-\alpha)=s^2\dfrac{n-r-1}{n-1}T^2_{r+1,\,n-1}(1-\alpha) \]

2. Intervalos de confiança simultâneos de \(100(1 - \alpha)\%\) para os \(\beta_i\) são dados por

\[ \text{IC}^{95\%}(\beta_i)=\left[\hat{\beta}_i \pm \sqrt{\widehat{\mathbb{V}}(\hat{\beta}_i)}\sqrt{(r + 1)F_{r + 1, \,n - r - 1}(1-\alpha)}\right]\\ i=1,2,\ldots,r \]

em que \(\widehat{\mathbb{V}}(\hat{\beta}_i)\) é o elemento diagonal de \(s^2(\mathbf{z}^{\prime}\mathbf{z})^{-1}\) correspondente a \(\hat{\beta}_i\).

\[\Diamond\]

O elipsoide de confiança é centrado na estimativa de máxima verossimilhança \(\hat{\boldsymbol{\beta}}\), e sua orientação e tamanho são determinados pelos autovalores e autovetores de \(\mathbf{z}^{\prime}\mathbf{z}\). Se um autovalor é quase zero, o elipsoide de confiança será muito longo na direção do autovetor correspondente.

Os praticantes muitas vezes ignoram a propriedade de confiança “simultânea” das estimativas de intervalo no Resultado 7.5. Em vez disso, eles substituem \((r + 1)F_{r + 1, n - r - 1}(1-\alpha)\) pelo valor t um de cada vez \(t_{n-r-1}(1-\alpha/2)\) e usam os intervalos

\[ \text{IC}^{95\%}(\beta_i)=\left[\hat{\beta}_i \pm t_{n-r-1}\left(1-\dfrac{\alpha}{2}\right)\sqrt{\widehat{\mathbb{V}}(\hat{\beta}_i)}\right]\\ i=1,2,\ldots,r \tag{7-11} \]

quando procuram por variáveis preditoras importantes.

Exemplo 7.4: Ajustando um modelo de regressão a dados imobiliários

Os dados de avaliação na Tabela 7.1 foram coletados de 20 residências em um bairro de Milwaukee, Wisconsin. Ajuste o modelo de regressão

\[ Y_j = \beta_0 + \beta_1 z_{j1} + \beta_2 z_{j2} + \varepsilon_j \]

em que \(z_{j1}\) é o tamanho total da habitação (em centenas de pés quadrados), \(z_{j2}\) é o valor avaliado (em milhares de dólares), e \(Y\) é o preço de venda (em milhares de dólares), a esses dados usando o método de mínimos quadrados. Um cálculo computacional produz

\[ (\mathbf{z}^{\prime}\mathbf{z})^{-1} = \begin{bmatrix} 5.1523 & \\ 0.2544 & 0.0512 \\ -0.1463 & -0.0172 & 0.0067 \end{bmatrix} \]

e

\[ \hat{\boldsymbol{\beta}} = (\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}\mathbf{y} = \begin{bmatrix} 30.967 \\ 2.634 \\ 0.045 \end{bmatrix} \]

Assim, a equação ajustada é

\[ \hat{y} = 30.967 + 2.634\;z_1 + 0.045\;z_2 \]

com \(s = 3.473\). Os números entre parênteses são os desvios-padrão estimados dos coeficientes de mínimos quadrados. Além disso, \(R^2 = 0.834\), indicando que os dados apresentam uma forte relação de regressão. Se os resíduos \(\hat{\boldsymbol{\varepsilon}}\) passarem nas verificações de diagnóstico descritas na Seção 7.6, a equação ajustada poderá ser usada para prever o preço de venda de outra casa no bairro a partir de seu tamanho e valor avaliado. Observamos que um intervalo de confiança de 95% para \(\beta_2\) é dado por

\[ \begin{align} \text{IC}^{95\%}(\beta_2)&=\left[\hat{\beta}_2 \pm t_{17}\left(1-\dfrac{0.05}{2}\right)\sqrt{\widehat{\mathbb{V}}(\hat{\beta}_2)}\right] \\ &= \left[0.045 \pm 2.110\times 0.285\right]\\ \text{IC}^{95\%}(\beta_2)&=[−0.556,0.647] \end{align} \]

O intervalo de confiança inclui \(\beta_2 = 0\).

Se o modelo é preditivo, a variável \(z_2\) pode ser eliminada do modelo de regressão e a análise repetida com a única variável preditora \(z_1\). Dado o tamanho da habitação, o valor avaliado parece agregar pouco à previsão do preço de venda. Se o modelo é explicativo, mesmo a variável \(z_2\) não sendo significante, ela permanece no modelo.

• JW6_Example7.4.R

invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))

Dados <- read.csv("JW6Data/T7-1.dat", sep = "", header = FALSE)
names(Dados) <- c("z1", "z2", "y")
print.data.frame(Dados)

      z1   z2     y
1  15.31 57.3  74.8
2  15.20 63.8  74.0
3  16.25 65.4  72.9
4  14.33 57.0  70.0
5  14.57 63.8  74.9
6  17.33 63.2  76.0
7  14.48 60.2  72.0
8  14.91 57.7  73.5
9  15.25 56.4  74.5
10 13.89 55.6  73.5
11 15.18 62.6  71.5
12 14.44 63.4  71.0
13 14.87 60.2  78.9
14 18.63 67.2  86.5
15 15.20 57.1  68.0
16 25.76 89.6 102.0
17 19.05 68.6  84.0
18 15.37 60.1  69.0
19 18.06 66.3  88.0
20 16.35 65.8  76.0

print(psych::describe(Dados))

   vars  n  mean   sd median trimmed  mad   min    max range skew kurtosis   se
z1    1 20 16.22 2.68  15.22   15.71 1.13 13.89  25.76 11.87 2.25     5.20 0.60
z2    2 20 63.06 7.39  62.90   61.94 4.67 55.60  89.60 34.00 2.16     5.52 1.65
y     3 20 76.55 8.07  74.25   75.25 3.71 68.00 102.00 34.00 1.67     2.45 1.80

boxplot(Dados)

GGally::ggpairs(Dados)

z <- model.matrix(~z1+z2,
                  data=Dados)
print(z)

   (Intercept)    z1   z2
1            1 15.31 57.3
2            1 15.20 63.8
3            1 16.25 65.4
4            1 14.33 57.0
5            1 14.57 63.8
6            1 17.33 63.2
7            1 14.48 60.2
8            1 14.91 57.7
9            1 15.25 56.4
10           1 13.89 55.6
11           1 15.18 62.6
12           1 14.44 63.4
13           1 14.87 60.2
14           1 18.63 67.2
15           1 15.20 57.1
16           1 25.76 89.6
17           1 19.05 68.6
18           1 15.37 60.1
19           1 18.06 66.3
20           1 16.35 65.8
attr(,"assign")
[1] 0 1 2

try(infomatrix.inv <- solve(t(z)%*%z))
round(svd(z)$d, 2)

[1] 293.22   4.93   0.44

# condition index < 30
CI <- sqrt(svd(z)$d[1]/svd(z)$d[length(svd(z)$d)])
CI

[1] 25.81951

y <- as.matrix(Dados["y"])
beta_hat <- infomatrix.inv%*%t(z)%*%y
beta_hat

                      y
(Intercept) 30.96656634
z1           2.63439962
z2           0.04518386

y_hat <- z%*%beta_hat
print(round(y_hat, 2))

        y
1   73.89
2   73.89
3   76.73
4   71.29
5   72.23
6   79.48
7   71.83
8   72.85
9   73.69
10  70.07
11  73.79
12  71.87
13  72.86
14  83.08
15  73.59
16 102.88
17  84.25
18  74.17
19  81.54
20  77.01

res <- y-y_hat
n <- dim(Dados)[1]
r <- length(grep("z",colnames(Dados)))
s <- as.numeric(sqrt(crossprod(res)/(n-r-1)))
s

[1] 3.472539

R2 <- 1 - crossprod(res)/((n-1)*var(y))
R2

         y
y 0.834397

alfa <- 0.05
LL_beta2 <- beta_hat[r+1] - 
  qt(1-alfa/(2*r), n-r-1)*sqrt(diag((s^2)*infomatrix.inv)[r+1])
UL_beta2 <- beta_hat[r+1] + 
  qt(1-alfa/(2*r), n-r-1)*sqrt(diag((s^2)*infomatrix.inv)[r+1])
cat("IC98.75%(beta_2) = [", 
    round(LL_beta2,3), ",", 
    round(UL_beta2,3), "]", 
    sep="")

IC98.75%(beta_2) = [-0.656,0.746]

# GLM
fit <- lm(y~z1+z2,
          data=Dados)
print(fit)

Call:
lm(formula = y ~ z1 + z2, data = Dados)

Coefficients:
(Intercept)           z1           z2  
   30.96657      2.63440      0.04518  

print(summary(fit))

Call:
lm(formula = y ~ z1 + z2, data = Dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.5894 -1.5411 -0.0718  1.3507  6.4605 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 30.96657    7.88221   3.929  0.00108 **
z1           2.63440    0.78560   3.353  0.00377 **
z2           0.04518    0.28518   0.158  0.87598   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.473 on 17 degrees of freedom
Multiple R-squared:  0.8344,    Adjusted R-squared:  0.8149 
F-statistic: 42.83 on 2 and 17 DF,  p-value: 2.302e-07

print(confint(fit), digits=2)

            2.5 % 97.5 %
(Intercept) 14.34  47.60
z1           0.98   4.29
z2          -0.56   0.65

print(round(predict(fit),2))

     1      2      3      4      5      6      7      8      9     10     11 
 73.89  73.89  76.73  71.29  72.23  79.48  71.83  72.85  73.69  70.07  73.79 
    12     13     14     15     16     17     18     19     20 
 71.87  72.86  83.08  73.59 102.88  84.25  74.17  81.54  77.01 

print(round(t(y_hat), 2))

      1     2     3     4     5     6     7     8     9    10    11    12    13
y 73.89 73.89 76.73 71.29 72.23 79.48 71.83 72.85 73.69 70.07 73.79 71.87 72.86
     14    15     16    17    18    19    20
y 83.08 73.59 102.88 84.25 74.17 81.54 77.01

# plot(fit)

alfa <- 0.05
# https://en.wikipedia.org/wiki/Partial_correlation
# https://en.wikiversity.org/wiki/Semi-partial_correlation
# IC95% de beta padronizado e sr2 (r semi-parcial^2) entre -1 e 1.
print(apaTables::apa.reg.table(fit,
                               prop.var.conf.level = 1-alfa))

Regression results using y as the criterion
 

   Predictor       b       b_95%_CI beta   beta_95%_CI sr2  sr2_95%_CI     r
 (Intercept) 30.97** [14.34, 47.60]                                         
          z1  2.63**   [0.98, 4.29] 0.88  [0.32, 1.43] .11 [-.03, .25] .91**
          z2    0.05  [-0.56, 0.65] 0.04 [-0.51, 0.59] .00 [-.01, .01] .85**
                                                                            
                                                                            
                                                                            
             Fit
                
                
                
     R2 = .834**
 95% CI[.60,.89]
                

Note. A significant b-weight indicates the beta-weight and semi-partial correlation are also significant.
b represents unstandardized regression weights. beta indicates the standardized regression weights. 
sr2 represents the semi-partial correlation squared. r represents the zero-order correlation.
Square brackets are used to enclose the lower and upper limits of a confidence interval.
* indicates p < .05. ** indicates p < .01.
 

fitst <- lm(scale(y)~scale(z1)+scale(z2),
          data=Dados)
print(fitst, digits=4)

Call:
lm(formula = scale(y) ~ scale(z1) + scale(z2), data = Dados)

Coefficients:
(Intercept)    scale(z1)    scale(z2)  
  6.319e-17    8.750e-01    4.134e-02  

# partial correlation
prz1 <- (cor(Dados$y,Dados$z1)-cor(Dados$y,Dados$z2)*cor(Dados$z2,Dados$z1))/
         sqrt((1-cor(Dados$y,Dados$z2)^2)*(1-cor(Dados$z2,Dados$z1)^2))
round(prz1,4)

[1] 0.631

pr.z1 <- as.numeric(ppcor::pcor.test(Dados$y, Dados$z1, Dados$z2)$estimate)
round(pr.z1,4)

[1] 0.631

prz2 <- (cor(Dados$y,Dados$z2)-cor(Dados$y,Dados$z1)*cor(Dados$z1,Dados$z2))/
  sqrt((1-cor(Dados$y,Dados$z1)^2)*(1-cor(Dados$z1,Dados$z2)^2))
round(prz2,4)

[1] 0.0384

pr.z2 <- as.numeric(ppcor::pcor.test(Dados$y, Dados$z2, Dados$z1)$estimate)
round(pr.z2,4)

[1] 0.0384

# semipartial = part correlation
srz1 <- (cor(Dados$y,Dados$z1)-cor(Dados$y,Dados$z2)*cor(Dados$z2,Dados$z1))/
         sqrt((1-cor(Dados$z2,Dados$z1)^2))
round(srz1,4)

[1] 0.331

sr.z1 <- as.numeric(ppcor::spcor.test(Dados$y, Dados$z1, Dados$z2)$estimate)
round(sr.z1,4)

[1] 0.331

sr2.z1 <- sr.z1^2
round(sr2.z1,4)

[1] 0.1095

srz2 <- (cor(Dados$y,Dados$z2)-cor(Dados$y,Dados$z1)*cor(Dados$z1,Dados$z2))/
         sqrt((1-cor(Dados$z2,Dados$z1)^2))
round(srz2,4)

[1] 0.0156

sr.z2 <- as.numeric(ppcor::spcor.test(Dados$y, Dados$z2, Dados$z1)$estimate)
round(sr.z2,4)

[1] 0.0156

sr2.z2 <- sr.z2^2
round(sr2.z2,6)

[1] 0.000245

# Kutner et al. (2005, p. 276)
# standardized beta
b1.std <- (cor(Dados$y,Dados$z1)-cor(Dados$y,Dados$z2)*cor(Dados$z2,Dados$z1))/
           (1-cor(Dados$z1,Dados$z2)^2)
round(b1.std,4)

[1] 0.875

b1std <-  srz1/sqrt((1-cor(Dados$z2,Dados$z1)^2))
round(b1std,4)

[1] 0.875

b2.std <- (cor(Dados$y,Dados$z2)-cor(Dados$y,Dados$z1)*cor(Dados$z1,Dados$z2))/
          (1-cor(Dados$z1,Dados$z2)^2)
round(b2.std,4)

[1] 0.0413

b2std <-  srz2/sqrt((1-cor(Dados$z2,Dados$z1)^2))
round(b2std,4)

[1] 0.0413

# beta padronizado não é uma correlação e os limites superior/inferior de 
# beta padronizado são +-1/sqrt((1-cor(Dados$z2,Dados$z1)^2))

# Plano de regressão
rgl::plot3d(fit, size = 5)
rgl::planes3d(0, 0, 1, 0, col="gray", alpha=0.2)
rgl::planes3d(0, 1, 0, 0, col="gray", alpha=0.2)
rgl::planes3d(1, 0, 0, 0, col="gray", alpha=0.2)
centroid <- colMeans(Dados)
rgl::points3d(centroid, size = 5, col = "purple3")
rgl::text3d(x = centroid,col = "purple3", texts = "centroide")
rgl::rglwidget()

\[\Diamond\]

Testes de Razão de Verossimilhança para os Parâmetros de Regressão

Parte da análise de regressão está preocupada em avaliar os efeitos de variáveis preditoras particulares na variável resposta. Uma hipótese nula de interesse afirma que certas variáveis \(z\) não influenciam a resposta \(Y\). Esses preditores serão rotulados \(z_{q+1}, z_{q+2}, \ldots, z_r\). Esta afirmação se traduz na hipótese estatística

\[ H_0:\beta_{q+1} = \beta_{q+2} = \cdots = \beta_r = 0 \]

OU

\[ H_0:\boldsymbol{\beta}_{(2)} = \mathbf{0} \tag{7-12} \]

em que \(\boldsymbol{\beta}_{(2)} = [\beta_{q+1}\; \beta_{q+2}\; \cdots\; \beta_r]^{\prime}\).

Definindo

\[ \mathbf{z} = \begin{bmatrix} \mathbf{z}_1 & \mathbf{z}_2 \end{bmatrix} \]

em que \(\mathbf{z}_1\) tem dimensões \(n \times (q+1)\) e \(\mathbf{z}_2\) tem dimensões \(n \times (r-q)\) e

\[ \boldsymbol{\beta}= \begin{bmatrix} \boldsymbol{\beta}_{(1)}\\ \boldsymbol{\beta}_{(2)} \end{bmatrix} \]

em que \(\boldsymbol{\beta}_{(1)}\) tem dimensões \((q+1) \times 1\) e \(\boldsymbol{\beta}_{(2)}\) tem dimensões \((r-q) \times 1\)

Podemos expressar o modelo linear geral como

\[ \begin{align} \mathbf{Y} &= \mathbf{z}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \\ &= \begin{bmatrix} \mathbf{z}_1 & \mathbf{z}_2 \end{bmatrix} \begin{bmatrix} \boldsymbol{\beta}_{(1)} \\ \boldsymbol{\beta}_{(2)} \end{bmatrix} + \boldsymbol{\varepsilon} \\ \mathbf{Y}&= \mathbf{z}_1\boldsymbol{\beta}_{(1)} + \mathbf{z}_2\boldsymbol{\beta}_{(2)} + \boldsymbol{\varepsilon} \end{align} \]

Sob a hipótese nula \(H_0: \boldsymbol{\beta}_{(2)} = 0\), temos \(\mathbf{Y} = \mathbf{z}_1\boldsymbol{\beta}_{(1)} + \boldsymbol{\varepsilon}\). O teste de razão de verossimilhança de \(H_0\) é baseado na

\[ \begin{align} \text{Soma Extra de Quadrados} &= \text{SS}_\text{res}(\mathbf{z}_1) - \text{SS}_\text{res}(\mathbf{z}) \\ \text{Soma Extra de Quadrados} &= \left(\mathbf{y} - \mathbf{z}_1\hat{\boldsymbol{\beta}}_{(1)}\right)^{\prime}\left(\mathbf{y} - \mathbf{z}_1\hat{\boldsymbol{\beta}}_{(1)}\right) - \left(\mathbf{y} - \mathbf{z}\hat{\boldsymbol{\beta}}\right)^{\prime}\left(\mathbf{y} - \mathbf{z}\hat{\boldsymbol{\beta}}\right) \end{align} \tag{7-13} \]

em que \(\hat{\boldsymbol{\beta}}_{(1)} = \left(\mathbf{z}_1^{\prime}\mathbf{z}_1\right)^{-1}\mathbf{z}_1^{\prime}\mathbf{Y}\).

Resultado 7.6. Seja \(\mathbf{z}\) de posto completo \(r + 1\) e \(\boldsymbol{\varepsilon}\) distribuído como \(\mathcal{N}_n(\mathbf{0},\sigma^2I)\). O teste de razão de verossimilhança de \(H_0: \boldsymbol{\beta}_{(2)} = \mathbf{0}\) é equivalente a um teste de \(H_0\) baseado na soma extra de quadrados em (7-13). Em particular, o teste de razão de verossimilhança rejeita \(H_0\) se

\[ \begin{align} F&=\dfrac{\dfrac{\mathbf{Y}^{\prime}(\mathbf{H}-\mathbf{H}_{(1)})\mathbf{Y}}{\text{tr}(\mathbf{H}-\mathbf{H}_{(1)})}}{\dfrac{\mathbf{Y}^{\prime}(\mathbf{I}-\mathbf{H})\mathbf{Y}}{\text{tr}(\mathbf{I}-\mathbf{H})}}> F_{\text{tr}(\mathbf{I}-\mathbf{H}),\,\text{tr}(\mathbf{H}-\mathbf{H}_{(1)})}(1-\alpha)\\ F&=\dfrac{\dfrac{\text{SS}_{\text{res}}(\mathbf{z}_1)-\text{SS}_{\text{res}}(\mathbf{z})}{\text{rank}(\mathbf{z})-\text{rank}(\mathbf{z}_1)}}{\dfrac{\text{SS}_{\text{res}}(\mathbf{z})}{n-\text{rank}(\mathbf{z})}}> F_{\text{rank}(\mathbf{z})-\text{rank}(\mathbf{z}_1),\,n-\text{rank}(\mathbf{z})}(1-\alpha) \end{align} \]

Sendo que \(\mathbf{I}-\mathbf{H}_{(1)} - (\mathbf{I}-\mathbf{H})= \mathbf{H}-\mathbf{H}_{(1)}\); \(\mathbf{z}\) e \(\mathbf{z}_1\) têm postos completos, \(\text{rank}(\mathbf{z})-\text{rank}(\mathbf{z}_1)=r+1-(q+1)=r-q\) e \(n-\text{rank}(\mathbf{z})=n-(r+1) = n-r-1\):

\[ F=\dfrac{\dfrac{\text{SS}_{\text{res}}(\mathbf{z}_1)-\text{SS}_{\text{res}}(\mathbf{z})}{r-q}}{\dfrac{\text{SS}_{\text{res}}(\mathbf{z})}{n-r-1}} > F_{r-q,\,n-r-1}(1-\alpha) \]

\[\Diamond\]

Comentário. O teste de razão de verossimilhança é implementado da seguinte forma. Para testar se todos os coeficientes em um subconjunto são zero, ajuste o modelo com e sem os termos correspondentes a esses coeficientes. A melhoria na soma dos quadrados dos resíduos (a soma extra dos quadrados) é comparada à soma dos quadrados dos resíduos para o modelo completo através da razão F. O mesmo procedimento se aplica mesmo em situações de análise de variância, em que \(\mathbf{z}\) não é de posto completo. Em situações em que \(\mathbf{z}\) não possui posto completo, \(\text{rank}(\mathbf{z})\) substitui \(r + 1\) e \(\text{rank}(\mathbf{z}_1)\) substitui \(q + 1\) no Resultado 7.6.

Mais geralmente, é possível formular hipóteses nulas concernentes a \(r - q\) combinações lineares de \(\boldsymbol{\beta}\) da forma \(H_0: \mathbf{C}\boldsymbol{\beta} = \mathbf{A}_0\). Seja a matriz \(\mathbf{C}\) de dimensões \((r - q) \times (r + 1)\) de posto completo, seja \(\mathbf{A}_0 = \mathbf{0}\), e considere

\[ H_0: \mathbf{C}\boldsymbol{\beta} = \mathbf{0} \]

Esta hipótese nula se reduz à escolha anterior quando \(\mathbf{C} = [\mathbf{0} \; \mathbf{I}]\), sendo que \(\mathbf{I}\) tem dimensões \((r-q)\times (r-q)\).

Sob o modelo completo,

\[ \mathbf{C}\hat{\boldsymbol{\beta}}\sim \mathcal{N}_{r-q}\left(\mathbf{C}\boldsymbol{\beta}, \sigma^2\mathbf{C}\mathbf{z}^{\prime}\mathbf{z}\mathbf{C}^{\prime}\right) \]

Rejeitamos \(H_0: \mathbf{C}\boldsymbol{\beta} = \mathbf{0}\) no nível \(\alpha\) se \(\mathbf{0}\) não estiver no elipsoide de confiança de 100(1 - \(\alpha\)% para \(\mathbf{C}\boldsymbol{\beta}\). Equivalentemente, rejeitamos \(H_0: \mathbf{C}\boldsymbol{\beta} = \mathbf{0}\) se

\[ F=\dfrac{\dfrac{(\mathbf{C}\hat{\boldsymbol{\beta}})^{\prime}\left(\mathbf{C}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{C}^{\prime}\right)\mathbf{C}\hat{\boldsymbol{\beta}}}{r-q}}{\dfrac{\text{SS}_{\text{res}}}{n-r-1}} > F_{r-q, \,n-r-1}(1-\alpha) \]

O teste em (7-14) é o teste de razão de verossimilhança, e o numerador na razão F é a soma extra dos quadrados dos resíduos incorridos ao ajustar o modelo, sujeito à restrição de que \(\mathbf{C}\boldsymbol{\beta} = \mathbf{0}\) (veja [23]).

O próximo exemplo ilustra como desenhos experimentais desbalanceados são facilmente tratados pela teoria geral que acabamos de descrever.

Exemplo 7.5: Testando a importância de preditores adicionais usando a abordagem de soma extra de quadrados

Clientes masculinos e femininos avaliaram o serviço em três estabelecimentos (locais) de uma grande rede de restaurantes. As avaliações de serviço foram convertidas em um índice. A Tabela 7.2 contém os dados para \(n = 18\) clientes. Cada ponto de dados na tabela é categorizado de acordo com o local (1, 2 ou 3) e sexo (masculino = 0 e feminino = 1). Essa categorização tem o formato de uma tabela bidimensional com números desiguais de observações por célula. Por exemplo, a combinação de local 1 e masculino tem 5 respostas, enquanto a combinação de local 2 e feminino tem 2 respostas. Introduzindo três variáveis dummy para representar o local e duas variáveis dummy para representar o gênero, podemos desenvolver um modelo de regressão relacionando o índice de serviço \(F\) ao local, gênero e sua “interação” usando a matriz de design.

O vetor de coeficientes pode ser representado como

\[ \boldsymbol{\beta}^{\prime} = [\beta_0, \beta_1, \beta_2, \beta_3, \tau_1, \tau_2, \gamma_{11}, \gamma_{12}, \gamma_{21}, \gamma_{22}, \gamma_{31}, \gamma_{32}] \]

em que os \(\beta_i (i > 0)\) representam os efeitos dos locais na determinação do serviço, os \(\tau_i\) representam os efeitos do sexo no índice de serviço, e os \(\gamma_{ij}\) representam os efeitos de interação entre local e sexo.

A matriz de design \(\mathbf{z}\) não é de posto completo. (Por exemplo, a coluna 1 é igual à soma das colunas 2-4 ou colunas 5-6.) Na verdade, \(\text{rank}(\mathbf{z}) = 6<12\).

Para o modelo completo, os resultados de um programa de computador dão

\[ \text{SS}_{\text{res}}(\mathbf{z}) = 2977.4 \]

e

\[ n - \text{rank}(\mathbf{z}) = 18 - 6 = 12 \]

O modelo sem os termos de interação tem a matriz de design \(\mathbf{z}_1\) consistindo nas primeiras seis colunas de \(\mathbf{z}\).

Descobrimos que

\[ \text{SS}_{\text{res}}(\mathbf{z}_1) = 3419.1 \]

com

\[ n - \text{rank}(\mathbf{z}_1) = 18 - 4 = 14 \]

Para testar

\[ H_0: \gamma_{11} = \gamma_{12} = \gamma_{21} = \gamma_{22} = \gamma_{31} = \gamma_{32} = 0 \]

(sem interação local-sexo), calculamos

\[ \begin{align} F&=\dfrac{\dfrac{\text{SS}_{\text{res}}(\mathbf{z}_1)-\text{SS}_{\text{res}}(\mathbf{z})}{\text{rank}(\mathbf{z})-\text{rank}(\mathbf{z}_1)}}{\dfrac{\text{SS}_{\text{res}}(\mathbf{z})}{n-\text{rank}(\mathbf{z})}}\\ &=\dfrac{\dfrac{3419.1 - 2977.4}{6-4}}{\dfrac{2977.4}{18-6}}\\ F&=0.89 \end{align} \]

\[ F_{2,12}(0.95)=3.89\quad p = 0.436 \]

Os valores p de soma extra de quadrados e do modelo completo com lm são iguais (\(p = 0.436\)). Note que os dados são desbalanceados.

• JW6_Example7.5.R

invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))

Dados <- read.csv("JW6Data/T7-2.dat", sep = "", header = FALSE)
names(Dados) <- c("F1", "F2", "y")
Dados$F1 <- factor(Dados$F1)
Dados$F2 <- factor(Dados$F2)
print.data.frame(Dados)

   F1 F2    y
1   1  0 15.2
2   1  0 21.2
3   1  0 27.3
4   1  0 21.2
5   1  0 21.2
6   1  1 36.4
7   1  1 92.4
8   2  0 27.3
9   2  0 15.2
10  2  0  9.1
11  2  0 18.2
12  2  0 50.0
13  2  1 44.0
14  2  1 63.6
15  3  0 15.2
16  3  0 30.3
17  3  1 36.4
18  3  1 40.9

# Dados desbalanceados
xtabs(~F1+F2, data=Dados)

   F2
F1  0 1
  1 5 2
  2 5 2
  3 2 2

print(psych::describeBy(y~F1+F2,
                        mat=1,
                        digits=2,
                        data=Dados))

   item group1 group2 vars n  mean    sd median trimmed   mad  min  max range
y1    1      1      0    1 5 21.22  4.28  21.20   21.22  0.00 15.2 27.3  12.1
y2    2      2      0    1 5 23.96 15.97  18.20   23.96 13.49  9.1 50.0  40.9
y3    3      3      0    1 2 22.75 10.68  22.75   22.75 11.19 15.2 30.3  15.1
y4    4      1      1    1 2 64.40 39.60  64.40   64.40 41.51 36.4 92.4  56.0
y5    5      2      1    1 2 53.80 13.86  53.80   53.80 14.53 44.0 63.6  19.6
y6    6      3      1    1 2 38.65  3.18  38.65   38.65  3.34 36.4 40.9   4.5
   skew kurtosis    se
y1 0.02    -1.40  1.91
y2 0.67    -1.41  7.14
y3 0.00    -2.75  7.55
y4 0.00    -2.75 28.00
y5 0.00    -2.75  9.80
y6 0.00    -2.75  2.25

boxplot(y~F1*F2,
        data=Dados)

GGally::ggpairs(Dados)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Solução 1: Soma extra de quadrados por lm e anova
fit <- lm(y~F1+F2+F2:F1,
          data=Dados)
print(summary(fit))

Call:
lm(formula = y ~ F1 + F2 + F2:F1, data = Dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-28.000  -7.168  -0.020   5.395  28.000 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   21.220      7.044   3.012  0.01082 * 
F12            2.740      9.962   0.275  0.78796   
F13            1.530     13.179   0.116  0.90950   
F21           43.180     13.179   3.276  0.00662 **
F12:F21      -13.340     18.638  -0.716  0.48784   
F13:F21      -27.280     20.538  -1.328  0.20879   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 15.75 on 12 degrees of freedom
Multiple R-squared:  0.5857,    Adjusted R-squared:  0.4131 
F-statistic: 3.393 on 5 and 12 DF,  p-value: 0.03849

print(car::Anova(fit))

Anova Table (Type II tests)

Response: y
          Sum Sq Df F value   Pr(>F)   
F1         247.0  2  0.4978 0.619890   
F2        3746.7  1 15.1005 0.002164 **
F1:F2      441.8  2  0.8902 0.436021   
Residuals 2977.4 12                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fit.princ <- lm(y~F1+F2,
                data=Dados)
print(summary(fit.princ))

Call:
lm(formula = y ~ F1 + F2, data = Dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-19.419  -9.185  -3.452   3.451  36.581 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   24.652      6.329   3.895  0.00162 **
F12           -1.071      8.353  -0.128  0.89976   
F13           -9.536      9.942  -0.959  0.35379   
F21           31.167      7.957   3.917  0.00155 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 15.63 on 14 degrees of freedom
Multiple R-squared:  0.5242,    Adjusted R-squared:  0.4223 
F-statistic: 5.142 on 3 and 14 DF,  p-value: 0.01323

print(car::Anova(fit.princ))

Anova Table (Type II tests)

Response: y
          Sum Sq Df F value  Pr(>F)   
F1         247.0  2  0.5057 0.61368   
F2        3746.7  1 15.3411 0.00155 **
Residuals 3419.1 14                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(anova(fit.princ, fit))

Analysis of Variance Table

Model 1: y ~ F1 + F2
Model 2: y ~ F1 + F2 + F2:F1
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     14 3419.1                           
2     12 2977.4  2    441.76 0.8902  0.436

# os valores p de soma extra de quadrados e do modelo completo são iguais (p = 0.436)

# Solução 2: Soma extra de quadrados por model.matrix
z0 <- model.matrix(~1,
                   data=Dados)
print(z0)

   (Intercept)
1            1
2            1
3            1
4            1
5            1
6            1
7            1
8            1
9            1
10           1
11           1
12           1
13           1
14           1
15           1
16           1
17           1
18           1
attr(,"assign")
[1] 0

zF1 <- model.matrix(~-1+F1,
                    data=Dados)
print(zF1)

   F11 F12 F13
1    1   0   0
2    1   0   0
3    1   0   0
4    1   0   0
5    1   0   0
6    1   0   0
7    1   0   0
8    0   1   0
9    0   1   0
10   0   1   0
11   0   1   0
12   0   1   0
13   0   1   0
14   0   1   0
15   0   0   1
16   0   0   1
17   0   0   1
18   0   0   1
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$F1
[1] "contr.treatment"

zF2 <- model.matrix(~-1+F2,
                    data=Dados)
print(zF2)

   F20 F21
1    1   0
2    1   0
3    1   0
4    1   0
5    1   0
6    0   1
7    0   1
8    1   0
9    1   0
10   1   0
11   1   0
12   1   0
13   0   1
14   0   1
15   1   0
16   1   0
17   0   1
18   0   1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$F2
[1] "contr.treatment"

zF12 <- model.matrix(~-1+F1:F2,
                     data=Dados)
print(zF12)

   F11:F20 F12:F20 F13:F20 F11:F21 F12:F21 F13:F21
1        1       0       0       0       0       0
2        1       0       0       0       0       0
3        1       0       0       0       0       0
4        1       0       0       0       0       0
5        1       0       0       0       0       0
6        0       0       0       1       0       0
7        0       0       0       1       0       0
8        0       1       0       0       0       0
9        0       1       0       0       0       0
10       0       1       0       0       0       0
11       0       1       0       0       0       0
12       0       1       0       0       0       0
13       0       0       0       0       1       0
14       0       0       0       0       1       0
15       0       0       1       0       0       0
16       0       0       1       0       0       0
17       0       0       0       0       0       1
18       0       0       0       0       0       1
attr(,"assign")
[1] 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$F1
[1] "contr.treatment"

attr(,"contrasts")$F2
[1] "contr.treatment"

z <- cbind(z0,zF1,zF2,zF12)
z

   (Intercept) F11 F12 F13 F20 F21 F11:F20 F12:F20 F13:F20 F11:F21 F12:F21
1            1   1   0   0   1   0       1       0       0       0       0
2            1   1   0   0   1   0       1       0       0       0       0
3            1   1   0   0   1   0       1       0       0       0       0
4            1   1   0   0   1   0       1       0       0       0       0
5            1   1   0   0   1   0       1       0       0       0       0
6            1   1   0   0   0   1       0       0       0       1       0
7            1   1   0   0   0   1       0       0       0       1       0
8            1   0   1   0   1   0       0       1       0       0       0
9            1   0   1   0   1   0       0       1       0       0       0
10           1   0   1   0   1   0       0       1       0       0       0
11           1   0   1   0   1   0       0       1       0       0       0
12           1   0   1   0   1   0       0       1       0       0       0
13           1   0   1   0   0   1       0       0       0       0       1
14           1   0   1   0   0   1       0       0       0       0       1
15           1   0   0   1   1   0       0       0       1       0       0
16           1   0   0   1   1   0       0       0       1       0       0
17           1   0   0   1   0   1       0       0       0       0       0
18           1   0   0   1   0   1       0       0       0       0       0
   F13:F21
1        0
2        0
3        0
4        0
5        0
6        0
7        0
8        0
9        0
10       0
11       0
12       0
13       0
14       0
15       0
16       0
17       1
18       1

try(solve(t(z)%*%z))

Error in solve.default(t(z) %*% z) : 
  sistema é computacionalmente singular: condição recíproca número = 8.73995e-19

round(svd(z)$d, 2)

 [1] 6.23 3.37 3.26 2.52 1.63 1.48 0.00 0.00 0.00 0.00 0.00 0.00

round(infomatrix.inv_z <- matlib::Ginv(t(z)%*%z), 2)

                                               
 [1,]  0.5 -0.5 -0.5 0 -0.5 0  0.5  0.5 0 0 0 0
 [2,] -0.5  1.0  0.5 0  0.5 0 -1.0 -0.5 0 0 0 0
 [3,] -0.5  0.5  1.0 0  0.5 0 -0.5 -1.0 0 0 0 0
 [4,]  0.0  0.0  0.0 0  0.0 0  0.0  0.0 0 0 0 0
 [5,] -0.5  0.5  0.5 0  1.0 0 -1.0 -1.0 0 0 0 0
 [6,]  0.0  0.0  0.0 0  0.0 0  0.0  0.0 0 0 0 0
 [7,]  0.5 -1.0 -0.5 0 -1.0 0  1.7  1.0 0 0 0 0
 [8,]  0.5 -0.5 -1.0 0 -1.0 0  1.0  1.7 0 0 0 0
 [9,]  0.0  0.0  0.0 0  0.0 0  0.0  0.0 0 0 0 0
[10,]  0.0  0.0  0.0 0  0.0 0  0.0  0.0 0 0 0 0
[11,]  0.0  0.0  0.0 0  0.0 0  0.0  0.0 0 0 0 0
[12,]  0.0  0.0  0.0 0  0.0 0  0.0  0.0 0 0 0 0

# condition index < 30
CI <- sqrt(svd(z)$d[1]/svd(z)$d[length(svd(z)$d)])
CI

[1] 1.208375e+16

y <- as.matrix(Dados["y"])
n <- nrow(Dados)
I <- diag(1,n)
SSres.z <- t(y)%*%(I-z%*%infomatrix.inv_z%*%t(z))%*%y
as.numeric(SSres.z)

[1] 2977.39

z1 <- z[,1:(ncol(z0)+ncol(zF1)+ncol(zF2))]
z1

   (Intercept) F11 F12 F13 F20 F21
1            1   1   0   0   1   0
2            1   1   0   0   1   0
3            1   1   0   0   1   0
4            1   1   0   0   1   0
5            1   1   0   0   1   0
6            1   1   0   0   0   1
7            1   1   0   0   0   1
8            1   0   1   0   1   0
9            1   0   1   0   1   0
10           1   0   1   0   1   0
11           1   0   1   0   1   0
12           1   0   1   0   1   0
13           1   0   1   0   0   1
14           1   0   1   0   0   1
15           1   0   0   1   1   0
16           1   0   0   1   1   0
17           1   0   0   1   0   1
18           1   0   0   1   0   1

try(solve(t(z1)%*%z1))

Error in solve.default(t(z1) %*% z1) : 
  sistema é computacionalmente singular: condição recíproca número = 3.65506e-18

round(svd(z1)$d, 2)

[1] 5.89 2.82 2.65 2.07 0.00 0.00

round(infomatrix.inv_z1 <- matlib::Ginv(t(z1)%*%z1), 2)

                                
[1,]  0.31 -0.22 -0.22 0 -0.13 0
[2,] -0.22  0.40  0.26 0 -0.06 0
[3,] -0.22  0.26  0.40 0 -0.06 0
[4,]  0.00  0.00  0.00 0  0.00 0
[5,] -0.13 -0.06 -0.06 0  0.26 0
[6,]  0.00  0.00  0.00 0  0.00 0

SSres.z1 <- t(y)%*%(I-z1%*%infomatrix.inv_z1%*%t(z1))%*%y
as.numeric(SSres.z1)

[1] 3419.149

# https://en.wikipedia.org/wiki/QR_decomposition
z.rank <- qr(z)$rank
z.rank

[1] 6

z1.rank <- qr(z1)$rank
z1.rank

[1] 4

dfnum <- z.rank-z1.rank
dfden <- n-z.rank
F <- as.numeric(((SSres.z1-SSres.z)/dfnum)/(SSres.z/dfden))
as.numeric(F)

[1] 0.8902283

alfa <- 0.05
Fcrit <- qf(1-alfa, dfnum, dfden)
F > Fcrit

[1] FALSE

pv <- 1-pf(F, dfnum, dfden)
cat("F(", dfnum,",",dfden,") = ", F, 
    " p = ", pv," ESS = ", SSres.z1-SSres.z,"\n", sep="")

F(2,12) = 0.8902283 p = 0.4360194 ESS = 441.7595

print(anova(fit.princ, fit))

Analysis of Variance Table

Model 1: y ~ F1 + F2
Model 2: y ~ F1 + F2 + F2:F1
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     14 3419.1                           
2     12 2977.4  2    441.76 0.8902  0.436

A razão \(F\) pode ser comparada com um ponto percentual apropriado de uma distribuição \(F\) com 2 e 12 graus de liberdade. Esta razão \(F\) não é significativa para qualquer nível de significância razoável \(\alpha\). Consequentemente, concluímos que o índice de serviço não depende de qualquer interação localização-sexo, e esses termos podem ser descartados do modelo.

Usando a abordagem de soma extra de quadrados, podemos verificar que não há diferença entre as localizações (sem efeito de localização), mas que o sexo é significante; ou seja, homens e mulheres não dão as mesmas classificações para o serviço.

Em situações de análise de variância em que as contagens de células são desiguais, a variação na resposta atribuível a diferentes variáveis preditoras e suas interações geralmente não pode ser separada em quantidades independentes. Para avaliar as influências relativas dos preditores na resposta neste caso, é necessário ajustar o modelo com e sem os termos em questão e calcular as estatísticas de teste \(F\) apropriadas.

\[\Diamond\]

Inferência a partir da Função de Regressão Estimada

Regressão é interpolação da VD no domínio observado das VI.

Uma vez que um investigador esteja satisfeito com o modelo de regressão ajustado, ele pode ser usado para resolver dois problemas de predição. Seja \(\mathbf{z}_0 = [1\;z_{01}\; \cdots\; z_{0r}]\) valores selecionados para as variáveis preditoras que pertencem ao seu domínio observado. Então, \(\mathbf{z}_0\) e \(\hat{\boldsymbol{\beta}}\) podem ser usados para (1) estimar a função de regressão \(\beta_0 + \beta_1 z_{01} + \cdots + \beta_r z_{0r}\) em \(\mathbf{z}_0\) e (2) para estimar o valor da resposta \(Y\) em \(\mathbf{z}_0\).

Estimando a Função de Regressão em \(\mathbf{z}_0\)

Seja \(Y\) o valor da resposta quando as variáveis preditoras têm valores \(\mathbf{z}_0 = [1\;z_{01}\; \cdots\; z_{0r}]\). De acordo com o modelo em (7-3), o valor esperado de \(Y\) é

\[ \mathbb{E}\left(Y | \mathbf{z}_0\right) = \sum_{i=0}^{r}{\beta_i z_{0i}} = \mathbf{z}_0^{\prime} \boldsymbol{\beta} \tag{7-15} \]

Sua estimativa de mínimos quadrados é \(\mathbf{z}_0^{\prime} \hat{\boldsymbol{\beta}}\).

Resultado 7.7. Para o modelo de regressão linear em (7-3), \(\mathbf{z}_0 \hat{\boldsymbol{\beta}}\) é o estimador linear não-viesado de \(\mathbb{E}(Y | \mathbf{z}_0)\) com variância mínima, \(\mathbb{V}(\mathbf{z}_0^{\prime} \hat{\boldsymbol{\beta}}) = \mathbf{z}_0 (\mathbf{z}^{\prime}\mathbf{z})^{-1} \mathbf{z}_0^{\prime} \sigma^2\). Se \(\boldsymbol{\varepsilon}\) têm distribuição multinormal, então um intervalo de confiança de \(100(1 - \alpha)\%\) para \(\mathbb{E}(Y | \mathbf{z}_0) = \mathbf{z}_0^{\prime} \boldsymbol{\beta}\) é fornecido por

\[ \text{IC}^{95\%}(\mathbb{E}(Y|\mathbf{z}_0))=\left[\mathbf{z}_0^{\prime} \hat{\boldsymbol{\beta}} \pm t_{n-r-1}(1-\alpha/2) s \sqrt{\mathbf{z}_0^{\prime} (\mathbf{z}^{\prime}\mathbf{z})^{-1} \mathbf{z}_0}\right] \]

\[\Diamond\]

Previsão de uma Nova Observação em \(\mathbf{z}_0\)

A previsão de uma nova observação, como \(Y\), em \(\mathbf{z}_0 = [1\; z_{01}\; \cdots\; z_{0r}]\) é mais incerta do que estimar o valor esperado de \(Y\). De acordo com o modelo de regressão de (7-3), temos

\[ Y = \mathbf{z}_0^{\prime}\boldsymbol{\beta} + \varepsilon_0 \]

ou

\[ \text{(nova resposta } Y\text{)} = \text{(valor esperado de } Y \text{ em } \mathbf{z}_0\text{)} + \text{(novo erro)} \]

em que \(\varepsilon_0\) é distribuído como \(\mathcal{N}(0, \sigma^2)\) e é independente de \(\varepsilon\) e, portanto, de \(\hat{\boldsymbol{\beta}}\) e \(s^2\). Os erros \(\boldsymbol{\varepsilon}\) influenciam os estimadores \(\hat{\boldsymbol{\beta}}\) e \(s^2\) através das respostas \(\mathbf{Y}\), mas \(\varepsilon_0\) não.

Resultado 7.8. Dado o modelo de regressão linear de (7-3), uma nova observação \(Y\) tem o preditor não viesado

\[ \mathbf{z}_0^{\prime} \hat{\boldsymbol{\beta}} = \sum_{i=0}^{r}{\hat{\beta}_i z_{0i}} \]

A variância do erro de previsão \(Y - \mathbf{z}_0^{\prime} \hat{\boldsymbol{\beta}}\) é

\[ \mathbb{V}\left(Y - \mathbf{z}_0^{\prime} \hat{\boldsymbol{\beta}}\right) = \sigma^2\left(1 + \mathbf{z}_0^{\prime}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}_0\right) \]

Quando os erros \(\boldsymbol{\varepsilon }\) têm uma distribuição multinormal, um intervalo de predição de \(100(1 - \alpha)\%\) para \(Y\) é dado por

\[ \text{IP}^{95\%}(Y|\mathbf{z}_0)=\left[\mathbf{z}_0^{\prime} \hat{\boldsymbol{\beta}} \pm t_{n-r-1}(1-\alpha/2)s\sqrt{1 + \mathbf{z}_0^{\prime}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}_0}\right] \]

\[\Diamond\]

O intervalo de predição para \(y_0\) é mais amplo do que o intervalo de confiança para estimar o valor da função de regressão \(\mathbb{E}(Y| \mathbf{z}_0) = \mathbf{z}_0^{\prime}\boldsymbol{\beta}\). A incerteza adicional na previsão de \(Y\), que é representada pelo termo extra \(s^2\) na expressão

\[ s^2\left(1 + \mathbf{z}_0^{\prime}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}_0\right) \]

vem da presença do termo de erro desconhecido \(\varepsilon_0\).

Exemplo 7.6: Intervalos de confiança para uma resposta média e uma resposta futura

Empresas considerando a compra de um computador devem primeiro avaliar suas futuras necessidades para determinar o equipamento adequado. Um cientista da computação coletou dados de sete locais semelhantes de empresas para que uma equação de previsão dos requisitos de hardware de computador para gestão de estoque pudesse ser desenvolvida. Os dados são fornecidos na Tabela 7.3 para

\[ \begin{align} z_1 &= \text{pedidos do cliente (milhar)}\\ z_2 &= \text{contagem de adição-remoção de itens (milhar)}\\ y &= \text{tempo da CPU (unidade central de processamento) (hora)} \end{align} \]

Construa um intervalo de confiança de 95% para o tempo médio da CPU, \(\mathbb{E}(Y|\mathbf{z}_0) = \beta_0 + \beta_1z_{01} + \beta_2z_{02}\) em \(\mathbf{z}_0^{\prime} = [1,130, 7.5]\). Além disso, encontre um intervalo de predição de 95% para o requisito de CPU de uma nova instalação correspondente ao mesmo \(\mathbf{z}_0\).

A função de regressão estimada é

\[ \hat{y} = 8.42 + 1.08\;z_1 + 0.42\;z_2 \]

e

\[ (\mathbf{z}^{\prime}\mathbf{z})^{-1} = \begin{bmatrix} 8.17969 & & \\ -.06411 & 0.00052 & \\ 0.08831 & -.00107 & 0.01440 \end{bmatrix} \]

com \(s = 1.204\) (residual.scale de predict). Consequentemente,

\[ \mathbf{z}_0^{\prime} \hat{\boldsymbol{\beta}} = 8.42 + 1.08\times 130 + 0.42\times7.5 = 151.84 \]

e

\[ s\sqrt{\mathbf{z}_0^{\prime}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}_0} = 1.204\times0.58928 = 0.73 \]

Temos \(t_4(0.025) = 2.776\), então o intervalo de confiança de 95% para o tempo médio da CPU em \(\mathbf{z}_0\) é

\[ \begin{align} \text{IC}^{95\%}(\mathbb{E}(Y|\mathbf{z}_0))&=\left[\mathbf{z}_0^{\prime} \hat{\boldsymbol{\beta}} \pm t_{n-r-1}(0.975) s\sqrt{\mathbf{z}_0^{\prime} (\mathbf{z}^{\prime}\mathbf{z})^{-1} \mathbf{z}_0}\right]\\ &= [151.84 \pm 2.776 \times 0.73]\\ \text{IC}^{95\%}(\mathbb{E}(Y|\mathbf{z}_0))&=[149.81, 153.87] \end{align} \]

Uma vez que

\[ s\sqrt{1 + \mathbf{z}_0^{\prime}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}_0} = 1.204\times1.16071 = 1.40 \]

um intervalo de predição de 95% para o tempo da CPU em uma nova instalação com condições \(\mathbf{z}_0\) é

\[ \begin{align} \text{IP}^{95\%}(Y|\mathbf{z}_0)&=\left[\mathbf{z}_0^{\prime} \hat{\boldsymbol{\beta}} \pm t_{n-r-1}(0.975)s\sqrt{1 + \mathbf{z}_0^{\prime}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}_0}\right] \\ &= [151.84 \pm 2.776\times1.40]\\ \text{IP}^{95\%}(Y|\mathbf{z}_0)&=[147.93,155.75] \end{align} \]

• JW6_Example7.6.R

invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))

Dados <- read.csv("JW6Data/T7-3.dat", sep = "", header = FALSE)
names(Dados) <- c("z1", "z2", "y")
print.data.frame(Dados)

     z1    z2     y
1 123.5 2.108 141.5
2 146.1 9.213 168.9
3 133.9 1.905 154.8
4 128.5 0.815 146.5
5 151.5 1.061 172.8
6 136.2 8.603 160.1
7  92.0 1.125 108.5

print(psych::describe(Dados))

   vars n   mean    sd median trimmed   mad    min    max range  skew kurtosis
z1    1 7 130.24 19.42 133.90  130.24 15.42  92.00 151.50  59.5 -0.83    -0.58
z2    2 7   3.55  3.70   1.91    3.55  1.25   0.81   9.21   8.4  0.72    -1.60
y     3 7 150.44 21.63 154.80  150.44 19.72 108.50 172.80  64.3 -0.79    -0.73
     se
z1 7.34
z2 1.40
y  8.18

boxplot(Dados)

GGally::ggpairs(Dados)

fit <- lm(y~z1+z2,
          data=Dados)
print(summary(fit))

Call:
lm(formula = y ~ z1 + z2, data = Dados)

Residuals:
      1       2       3       4       5       6       7 
-1.0632 -1.0315  1.1007 -0.9151  0.4650  1.1066  0.3375 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.42369    3.44328   2.446   0.0707 .  
z1           1.07898    0.02749  39.249 2.52e-06 ***
z2           0.41989    0.14447   2.906   0.0438 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.204 on 4 degrees of freedom
Multiple R-squared:  0.9979,    Adjusted R-squared:  0.9969 
F-statistic: 966.5 on 2 and 4 DF,  p-value: 4.265e-06

print(car::Anova(fit))

Anova Table (Type II tests)

Response: y
           Sum Sq Df   F value    Pr(>F)    
z1        2232.86  1 1540.4714 2.517e-06 ***
z2          12.24  1    8.4467   0.04384 *  
Residuals    5.80  4                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

alfa <- 0.05
# IC95% de beta (correlação parcial) e sr2 entre -1 e 1.
print(apaTables::apa.reg.table(fit,
                               prop.var.conf.level = 1-alfa))

Regression results using y as the criterion
 

   Predictor      b       b_95%_CI beta  beta_95%_CI sr2  sr2_95%_CI     r
 (Intercept)   8.42 [-1.14, 17.98]                                        
          z1 1.08**   [1.00, 1.16] 0.97 [0.90, 1.04] .80 [.26, 1.33] 1..**
          z2  0.42*   [0.02, 0.82] 0.07 [0.00, 0.14] .00 [-.00, .01]   .45
                                                                          
                                                                          
                                                                          
              Fit
                 
                 
                 
      R2 = .998**
 95% CI[.97,1.00]
                 

Note. A significant b-weight indicates the beta-weight and semi-partial correlation are also significant.
b represents unstandardized regression weights. beta indicates the standardized regression weights. 
sr2 represents the semi-partial correlation squared. r represents the zero-order correlation.
Square brackets are used to enclose the lower and upper limits of a confidence interval.
* indicates p < .05. ** indicates p < .01.
 

z <- model.matrix(~z1+z2,
                  data=Dados)
print(z)

  (Intercept)    z1    z2
1           1 123.5 2.108
2           1 146.1 9.213
3           1 133.9 1.905
4           1 128.5 0.815
5           1 151.5 1.061
6           1 136.2 8.603
7           1  92.0 1.125
attr(,"assign")
[1] 0 1 2

try(infomatrix.inv <- solve(t(z)%*%z))
round(svd(z)$d, 2)

[1] 348.01   8.62   0.35

# condition index < 30
CI <- sqrt(svd(z)$d[1]/svd(z)$d[length(svd(z)$d)])
CI

[1] 31.54983

z_0 <- c(1, 130, 7.5)
y <- as.matrix(Dados["y"])
beta_hat <- infomatrix.inv%*%t(z)%*%y 
beta_hat

                    y
(Intercept) 8.4236890
z1          1.0789825
z2          0.4198885

n <- nrow(Dados)
z.rank <- qr(z)$rank
df <- n - z.rank
df

[1] 4

I <- diag(1,n)
SSres <- t(y)%*%(I-z%*%infomatrix.inv%*%t(z))%*%y
s <- sqrt(SSres/df)
as.numeric(s)

[1] 1.203937

y.mean_z0 <- as.numeric(t(z_0)%*%beta_hat)
cat("E^(Y|z0) = ", round(y.mean_z0,2), sep="")

E^(Y|z0) = 151.84

new <- data.frame(z1=130, z2=7.5)
out <- predict(fit, new, se.fit = TRUE)
as.numeric(out$fit)

[1] 151.8406

print(out)

$fit
       1 
151.8406 

$se.fit
[1] 0.7322769

$df
[1] 4

$residual.scale
[1] 1.203937

LL.conf <- t(z_0)%*%beta_hat - qt(1-alfa/2, df)*s*
           sqrt(t(z_0)%*%infomatrix.inv%*%z_0)
UL.conf <- t(z_0)%*%beta_hat + qt(1-alfa/2, df)*s*
           sqrt(t(z_0)%*%infomatrix.inv%*%z_0)
cat("IC95%(E(Y|z0)) = [", round(LL.conf,2), ",", round(UL.conf,2),
    "]", sep="")

IC95%(E(Y|z0)) = [149.81,153.87]

predict(fit, new, interval="confidence", level=1-alfa)

       fit      lwr      upr
1 151.8406 149.8075 153.8737

LL.pred <- t(z_0)%*%beta_hat - qt(1-alfa/2, df)*s*
           sqrt(1+t(z_0)%*%infomatrix.inv%*%z_0)
UL.pred <- t(z_0)%*%beta_hat + qt(1-alfa/2, df)*s*
           sqrt(1+t(z_0)%*%infomatrix.inv%*%z_0)
cat("IP95%(E(Y|z0)) = [", round(LL.pred,2), ",", round(UL.pred,2),
    "]", sep="")

IP95%(E(Y|z0)) = [147.93,155.75]

predict(fit, new, interval="prediction", level=1-alfa)

       fit      lwr     upr
1 151.8406 147.9282 155.753

Verificação do Modelo e Outros Aspectos da Regressão

O Modelo é Adequado?

Assumindo que o modelo é “correto”, usamos a função de regressão estimada para fazer inferências. Claro, é imperativo examinar a adequação do modelo antes que a função estimada se torne parte permanente do aparelho de tomada de decisão.

Todas as informações da amostra sobre falta de ajuste estão contidas nos resíduos

\[ \hat{\epsilon_j} = Y_j - \sum_{i=1}^{r}{\hat{\beta_i} z_{ji}} \\ j=1,2,\ldots,n \]

ou

\[ \begin{align} \hat{\boldsymbol{\varepsilon}} &= \mathbf{Y} - \mathbf{z}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}\mathbf{Y} \\ \hat{\boldsymbol{\varepsilon}} &= (\mathbf{I} - \mathbf{H})\mathbf{Y} \end{align} \tag{7-16} \]

Se o modelo é válido, cada resíduo é uma estimativa do erro \(\varepsilon_j\), que se presume ser uma variável aleatória normal com média zero e variância \(\sigma^2\). Embora os resíduos \(\hat{\boldsymbol{\varepsilon}}\) tenham valor esperado 0, sua matriz de covariância \(\sigma^2[\mathbf{I} - \mathbf{z}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}] = \sigma^2[\mathbf{I} - \mathbf{H}]\) não é diagonal. Os resíduos têm variâncias desiguais e correlações diferentes de zero. Felizmente, as correlações geralmente são pequenas e as variâncias são quase iguais.

Como os resíduos \(\hat{\boldsymbol{\varepsilon}}\) têm matriz de covariância \(\sigma^2[\mathbf{I} - \mathbf{H}]\), as variâncias dos \(\varepsilon_j\) podem variar muito se os elementos diagonais de \(\mathbf{H}\), as alavancas \(h_{jj}\), forem substancialmente diferentes. Consequentemente, muitos estatísticos preferem diagnósticos gráficos baseados em resíduos studentizados. Usando o quadrado médio residual \(s^2\) como uma estimativa de \(\sigma^2\), temos

\[ \widehat{\mathbb{V}}(\hat{\varepsilon}_j) = s^2(1 - h_{jj}) \\ j = 1, 2, \ldots, n \tag{7-17} \]

e os resíduos studentizados são

\[ \hat{\varepsilon}_j^* = \dfrac{\hat{\varepsilon}_j}{\sqrt{\widehat{\mathbb{V}}(\hat{\varepsilon}_j)}}=\dfrac{\hat{\varepsilon}_j}{s\sqrt{1 - h_{jj}}} \tag{7-18} \]

Esperamos que os resíduos studentizados pareçam, aproximadamente, como amostras independentes de uma distribuição \(\mathcal{N}(0, 1)\). Alguns programas estatísticos vão um passo além e studentizam \(\hat{\varepsilon}_j\) usando a variância estimada excluindo uma observação \(s^2(j)\), que é o quadrado médio residual quando a observação \(j\)-ésima é removida da análise.

Alavancagem e Influência

Embora uma análise de resíduos seja útil para avaliar o ajuste de um modelo, desvios do modelo de regressão muitas vezes são ocultados pelo processo de ajuste. Por exemplo, pode haver “outliers” tanto na resposta quanto nas variáveis explicativas que podem ter um efeito considerável na análise, mas que não são facilmente detectados a partir de uma análise de gráficos de resíduos. Na verdade, esses outliers podem determinar o ajuste.

A alavancagem \(h_{jj}\), o elemento diagonal \((j, j)\) de \(\mathbf{H} = \mathbf{z}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}\), pode ser interpretada de duas maneiras relacionadas. Primeiro, a alavancagem está associada ao \(i\)-ésimo ponto dos dados e mede, no espaço das variáveis explicativas, o quão distante a observação \(j\)-ésima está das outras \(n - 1\) observações. Para a regressão linear simples com uma variável explicativa \(z\),

\[ \begin{align} h_{jj}&=\dfrac{1}{n}\left(1+\left(\dfrac{z_j-\bar{z}}{s}\right)^2\right) \\ \text{sendo que:}\\ s^2&=\dfrac{1}{n}\sum_{j=1}^{n}{(z_j-\bar{z})^2}\end{align} \]

A alavancagem média é \(\bar{h}=(r + 1)/n\). (Veja o Exercício 7.8.)

Segundo, a alavancagem é uma medida do “puxão” que um único caso exerce no ajuste. O vetor de valores previstos é

\[ \begin{align} \hat{\mathbf{Y}} &= \mathbf{z}\hat{\boldsymbol{\beta}}\\ &=\mathbf{z}(\mathbf{z}^{\prime}\mathbf{z})^{-1}\mathbf{z}^{\prime}\mathbf{Y}\\ \hat{\mathbf{Y}} &= \mathbf{H}\mathbf{Y} \end{align} \]

em que a \(j\)-ésima linha expressa o valor ajustado \(\hat{y}_j\) em termos das observações como

\[ \hat{y}_j = h_{jj}y_j+\sum_{k\ne j}^{}{h_{jk}y_k} \]

Desde que todos os outros valores de \(y\) sejam mantidos fixos,

\[ \text{(mudança em } \hat{y}_j) = h_{jj}\times\text{(mudança em } y_j) \]

Se a alavancagem é grande em relação aos outros \(h_{jk}\), então \(y_j\) será um grande contribuidor para o valor previsto \(\hat{y}_j\).

Observações que afetam significativamente as inferências extraídas dos dados são consideradas influentes. Métodos para avaliar a influência geralmente são baseados na mudança no vetor de estimativas de parâmetros, \(\hat{\boldsymbol{\beta}}\), quando observações são excluídas. Gráficos baseados em estatísticas de alavancagem e influência e seu uso na verificação diagnóstica de modelos de regressão são descritos em [3], [5] e [10]. Essas referências são recomendadas para qualquer pessoa envolvida em uma análise de modelos de regressão.

Se, após as verificações diagnósticas, nenhuma violação séria das suposições for detectada, podemos fazer inferências sobre \(\boldsymbol{\beta}\) e os futuros valores de \(Y\) com alguma certeza de que não seremos enganados.

• JW6_Cap7_leverage.R

# Criar dados simulados
set.seed(123)
n <- 10
x <- rnorm(n)
y <- 2 * x + rnorm(n)
Dados <- data.frame(x,y)
print(psych::describe(Dados))

  vars  n mean   sd median trimmed  mad   min  max range skew kurtosis   se
x    1 10 0.07 0.95  -0.08    0.04 0.76 -1.27 1.72  2.98 0.52    -1.08 0.30
y    2 10 0.36 2.65   0.00    0.36 1.51 -4.50 5.22  9.71 0.14    -0.50 0.84

boxplot(Dados)

GGally::ggpairs(Dados)

DescTools::PlotBag(Dados$x, Dados$y)

# Ajustar um modelo de regressão linear
modelo <- lm(y ~ x, data=Dados)
print(summary(modelo))

Call:
lm(formula = y ~ x, data = Dados)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.33303 -0.64421 -0.02448  0.49596  1.41472 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.1617     0.2852   0.567    0.586    
x             2.6287     0.3141   8.368 3.15e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8988 on 8 degrees of freedom
Multiple R-squared:  0.8975,    Adjusted R-squared:  0.8847 
F-statistic: 70.03 on 1 and 8 DF,  p-value: 3.153e-05

alfa <- 0.05
# IC95% de beta (correlação parcial) e sr2 entre -1 e 1.
print(apaTables::apa.reg.table(modelo,
                               prop.var.conf.level = 1-alfa))

Regression results using y as the criterion
 

   Predictor      b      b_95%_CI beta  beta_95%_CI sr2 sr2_95%_CI     r
 (Intercept)   0.16 [-0.50, 0.82]                                       
           x 2.63**  [1.90, 3.35] 0.95 [0.69, 1.21] .90 [.60, .94] .95**
                                                                        
                                                                        
                                                                        
             Fit
                
                
     R2 = .897**
 95% CI[.60,.94]
                

Note. A significant b-weight indicates the beta-weight and semi-partial correlation are also significant.
b represents unstandardized regression weights. beta indicates the standardized regression weights. 
sr2 represents the semi-partial correlation squared. r represents the zero-order correlation.
Square brackets are used to enclose the lower and upper limits of a confidence interval.
* indicates p < .05. ** indicates p < .01.
 

# Calcular valores de alavancagem
leverage <- hatvalues(modelo)
print.data.frame(cbind(Dados, leverage), digits=2)

        x     y leverage
1  -0.560  0.10     0.15
2  -0.230 -0.10     0.11
3   1.559  3.52     0.37
4   0.071  0.25     0.10
5   0.129 -0.30     0.10
6   1.715  5.22     0.43
7   0.461  1.42     0.12
8  -1.265 -4.50     0.32
9  -0.687 -0.67     0.17
10 -0.446 -1.36     0.13

cat("Limiar de alavancagem = ", 2*mean(leverage), sep="")

Limiar de alavancagem = 0.4

# Plotar valores de alavancagem
# plot(modelo)
plot(leverage, 
     ylab = "Valores de Alavancagem", 
     xlab = "Observação", pch = 16)
abline(h = 2*mean(leverage), col = "red")

# Criar um dataframe para os dados
dados <- data.frame(x = x, y = y, leverage = leverage)

# Criar o gráfico
grafico <- ggplot2::ggplot(dados, ggplot2::aes(x = x, y = y)) +
  ggplot2::geom_point() +
  ggplot2::geom_smooth(method = "lm", se = FALSE, color = "blue") +
  ggplot2::geom_text(ggplot2::aes(label = sprintf("%.2f", leverage)), vjust = 1.3) +
  ggplot2::labs(title = "Dados, Reta de Regressão e Valores de Alavancagem",
       x = "x",
       y = "y") 
# Exibir o gráfico
print(grafico)

`geom_smooth()` using formula = 'y ~ x'

Resultado 7.10

Demonstração: J&W6e, Solutions Manual, 2007, p. 166-7

Questões

7.1, 7.2, 7.3, 7.6, 7.8, 7.9, 7.11, 7.12, 7.21, 7.22, 7.24, 7.25, 7.27

• Regression effective degrees of freedom: Wikipedia

Códigos em R

invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
dados <- read.csv("JW6Data/T5-1.DAT", sep = "", header = FALSE)
colnames(dados) <- c("Suor", "Na", "K")
print(psych::describe(dados))

     vars  n  mean    sd median trimmed   mad  min  max range  skew kurtosis
Suor    1 20  4.64  1.70   4.50    4.58  1.48  1.5  8.5   7.0  0.40    -0.39
Na      2 20 45.40 14.13  47.30   45.81 12.82 13.5 71.6  58.1 -0.32    -0.44
K       3 20  9.96  1.90   9.75    9.87  2.22  7.1 14.0   6.9  0.34    -1.00
       se
Suor 0.38
Na   3.16
K    0.43

print(cov(dados), digits = 3)

      Suor     Na     K
Suor  2.88  10.01 -1.81
Na   10.01 199.79 -5.64
K    -1.81  -5.64  3.63

car::Boxplot(dados)

rgl::plot3d(lm(K~Suor+Na, data = dados), size = 5)
rgl::planes3d(0, 0, 1, 0, col="gray", alpha=0.2)
rgl::planes3d(0, 1, 0, 0, col="gray", alpha=0.2)
rgl::planes3d(1, 0, 0, 0, col="gray", alpha=0.2)
centroid <- colMeans(dados)
rgl::points3d(centroid, size = 5, col = "purple3")
rgl::text3d(x = centroid,col = "purple3", texts = "centroide")
rgl::rglwidget()

sigma <- matrix(c(4,3,3,4), ncol = 2, nrow = 2)
mu <- c(5, 5)
n <- 1000
set.seed(123)
x <- mvtnorm::rmvnorm(n = n, mean = mu, sigma = sigma)
d <- data.frame(x)
p2 <- ggplot2::ggplot(d, 
                      ggplot2::aes(x = X1, y = X2)) +
  ggplot2::geom_point(alpha = .5) +
  ggplot2::geom_density_2d()
p2

# Regressão3D
x1 <- rep(seq(0,10,0.5),100)
 x2 <- rep(seq(0,10,0.5),each=100)
 par(mfrow=c(2,2))
 Ey1 <- 83 + 9*x1 + 6*x2
 scatterplot3d::scatterplot3d(x1,x2,Ey1,highlight.3d=TRUE,
                              xlim=c(0,10),
                ylim=c(0,10),zlim=c(0,240),
                xlab=expression(x[1]),
                ylab=expression(x[2]),
                zlab="E(y)",main =
                expression(paste("A 3-d plot for ", 
                                 E(Y*"|"*x,beta) == 83 + 9*x[1]
                                     + 6*x[2])),z.ticklabs="")
 Ey2 <- 83 + 9*x1 + 6*x2 + 3*x1*x2
 scatterplot3d::scatterplot3d(x1,x2,Ey2,highlight.3d=TRUE,
                              xlim=c(0,10),
                ylim=c(0,10),zlim=c(0,600),
                xlab=expression(x[1]),
                ylab=expression(x[2]),zlab="E(y)",main =
                expression(paste("A 3-d plot for ",
                                 E(Y*"|"*x,beta)== 83 + 9*x[1]
                                     + 6*x[2] + 3*x[1]*x[2])),
                z.ticklabs="")
 Ey3 <- 83 + 9*x1 + 6*x2 + 2*x1^4 + 3*x2^3 + 3*x1*x2
 scatterplot3d::scatterplot3d(x1,x2,Ey3,highlight.3d=TRUE,
                              xlim=c(0,10),
                ylim=c(0,10),
                zlim=c(0,25000),
                xlab=expression(x[1]),
                ylab=expression(x[2]),
                zlab="E(y)",main =
                expression(paste("A 3-d plot for ",
                                 E(Y*"|"*x,beta)== 83 + 9*x[1]
                                    + 6*x[2] + 2*x[1]^4 + 3*x[2]^3 + 3*x[1]*x[2])),
                z.ticklabs="")
 Ey4 <- 83 + 9*x1 + 6*x2 - 2*x1^4 - 3*x2^3 + 3*x1*x2
 scatterplot3d::scatterplot3d(x1,x2,Ey4,highlight.3d=TRUE,
                              xlim=c(0,10),
                ylim=c(0,10),zlim=c(-23000,100),
                xlab=expression(x[1]),
                ylab=expression(x[2]),
                zlab="E(y)",main =
                expression(paste("A 3-d plot for ",
                                 E(Y*"|"*x,beta)==83 + 9*x[1]
                                     + 6*x[2] - 2*x[1]^4 - 3*x[2]^3 + 3*x[1]*x[2])),
                z.ticklabs="")

par(mfrow=c(2,2))
x1=x2=seq(from=0,to=10,by=0.2)
ey1 <- function(a,b) 83 + 9*a + 6*b
Ey1 <- outer(x1,x2,ey1)
contour(x1,x2,Ey1,main = expression(paste("Cantour plot for ",
                                          E(Y*"|"*x,beta) ==83 + 9*x[1]+ 6*x[2])))
ey2 <- function(a,b) 83 + 9*a + 6*b + 3*a*b
Ey2 <- outer(x1,x2,ey2)
contour(x1,x2,Ey2,main = expression(paste("Cantour plot for ",
                                          E(Y*"|"*x,beta)==83 + 9*x[1]+ 6*x[2] + 3*x[1]*x[2])))
ey3 <- function(a,b) 83 + 9*a + 6*b + 2*a^4 + 3*b^3 + 3*a*b
Ey3 <- outer(x1,x2,ey3)
contour(x1,x2,Ey3,main = expression(paste("Cantour plot for ",
                                          E(Y*"|"*x,beta)==83 + 9*x[1] + 6*x[2] + 2*x[1]^4 + 3*x[2]^3 + 3*x[1]*x[2])))
ey4 <- function(a,b) 83 + 9*a + 6*b - 2*a^4 - 3*b^3 + 3*a*b
Ey4 <- outer(x1,x2,ey4)
contour(x1,x2,Ey4,main = expression(paste("Cantour plot for ",
                                          E(Y*"|"*x,beta)==83 + 9*x[1] + 6*x[2] - 2*x[1]^4 - 3*x[2]^3 + 3*x[1]*x[2])))

Bibliografia

• Livro-texto:
  • JOHNSON, RA & WICHERN, DW (2007) Applied multivariate statistical analysis. 6^th ed. NJ: Prentice Hall.
  • JOHNSON, RA & WICHERN, DW (2007) Applied multivariate statistical analysis: Solutions Manual. 6^th ed. NJ: Prentice Hall.
  • JOHNSON, RA & WICHERN, DW (2002) Applied multivariate statistical analysis. 5^th ed. NJ: Prentice Hall.
• ANDERSON, TW (2003) An introduction to multivariate statistical analysis. 3^rd ed. NY: Wiley.
• EVERITT, BS & HOTHORN, T (2011) An introduction to applied multivariate analysis with R. USA: Springer.
• FLURY, B (1997) A first course in multivariate statistics. NY: Springer-Verlag.
• HARDLE, W & SIMAR, L (2015) Applied multivariate statistical analysis. 4^th ed. NY: Springer-Verlag.
• JAMES, G; WITTEN, D; HASTIE, T; & TIBSHIRANI, R (2021) An Introduction to statistical learning with applications in R. 2^nd ed. USA: Springer.
• KUTNER, MH et al. (2005) Applied linear statistical model. 5th ed. NY: McGraw-Hill/Irwin.
• MARDIA, KV, KENT, JT & BIBBY, JM (2003) Multivariate analysis. UK: Academic.
• MORRISON, DF (1990) Multivariate statistical methods. 3^rd ed. NY: McGraw-Hill.
• REIS, E (2001) Estatística multivariada aplicada. 2^a ed. Lisboa: Sílabo.
• REIS, SF (1988) Morfometria e Estatística Multivariada em Biologia Evolutiva. Revista Brasileira de Zoologia 5(4): 571-80.
• RENCHER, AC & CHRISTENSEN, WF (2012) Methods of multivariate analysis. 3^rd ed. NJ: Wiley.
• SAVILLE, DJ & WOOD, GR (1996) Statistical methods: a geometric primer.USA: Lawrence Erlbaum.
• SOUZA, J (1997) Estatística econômica e social. RJ: Campus.
• SRIVASTAVA, MS & CARTER, EM (1983) An Introduction to applied multivariate Statistics. North Holland.
• SRIVASTAVA, MS (2006) Methods of multivariate statistics. NY: Wiley.
• WICKENS, TD (1995) The geometry of multivariate Statistics. USA: Lawrence Erlbaum.

Bibliografia adicional

• Batschelet, E (1978) Introdução à matemática para biocientistas. Tradução da 2ª ed. São Paulo: EDUSP e Rio de Janeiro: Interciência.
• Batschelet, E (1979) Introduction to mathematics for life scientists. 3rd ed. NY: Springer.
• Bickel, PJ & Doksum, KA (2001) Mathematical Statistics: Basic Ideas and Selected Topics, Volume I. USA: CRC.
• Bilodeau, M & Brenner, D (1999) Theory of Multivariate Statistics. USA: Springer. T^2 de Hotelling
• Denis, DJ (2021) Applied Univariate, Bivariate, and Multivariate Statistics Using Python: A Beginners Guide to Advanced Data Analysis. NJ: Wiley.
• Genz, A (1992) Numerical Computation of Multivariate Normal Probabilities. Journal of Computational and Graphical Statistics 1(2): 141–149. https://doi.org/10.2307/1390838
• Hardle, W & Hlavka, Z (2007) Multivariate Statistics - Exercises and Solutions. USA: Springer.
• KASSAMBARA, A (2017) Practical guide to cluster analysis in R: Unsupervised machine learning. STHDA: http://www.sthda.com.
• Kollo. T & von Rosen, D (2005) Advanced Multivariate Statistics with Matrices. USA: Springer.
• MAIR, P (2018) Modern psychometrics with R. USA: Springer.
• Manly, BFJ & Alberto, JAN (2017) Multivariate Statistical Methods: A Primer using R. 4^th ed. USA: CRC.
  
• Matloff, N (2020) Probability and Statistics for Data Science: Math + R Data_. USA: CRC.
• Mirman, D (2014) Growth Curve Analysis and Visualization Using R. USA: CRC.
• Pessoa, DGC (2008) Exemplos do livro de Johnson e Wichern.
• Rawlings, JO et al. (1998) Applied Regression Analysis: A Research Tool. USA: Springer. T^2 de Hotelling
• Schumacker, RE (2016) Using R With Multivariate Statistics. USA: SAGE.
• Siqueira, JO (2012) Fundamentos de métodos quantitativos: aplicados em Administração, Economia, Contabilidade e Atuária usando WolframAlpha e SCILAB. São Paulo: Saraiva. Soluções dos exercícios em https://www.researchgate.net/publication/326533772_Fundamentos_de_metodos_quantitativos_-_Solucoes.
• Tay, A (2018) OLS using Matrix Algebra. http://www.mysmu.edu/faculty/anthonytay/MFE/OLS_using_Matrix_Algebra.pdf
• Teichroew, D (1964) An introduction to management science: deterministic models. NJ: Wiley.
• Timm, NH (2002) Applied Multivariate Analysis. USA: Springer. T^2 de Hotelling
• Zelterman, D (2015) Applied Multivariate Statistics with R. USA: Springer.

Aplicativos

• OpenAI (2021) ChatGPT. https://openai.com/research/chatgpt
• Wolfram Research (n.d.) Wolfram|Alpha. https://www.wolframalpha.com/
• Scilab Enterprises (n.d.) Scilab - Open source software for numerical computation. https://www.scilab.org/

Análise estatística multivariada em R e Python na internet

• Chaves Neto, A. CE-704 ANÁLISE MULTIVARIADA APLICADA À PESQUISA. https://docs.ufpr.br/~soniaisoldi/ce090/CE076AM_2010.pdf

• ME 731 - Métodos em Análise Multivariada em R. https://www.ime.unicamp.br/~cnaber/Material_AM_2S_2020.htm

• Nogueira, FE (2007) Modelos de regressão multivariada. Dissertação (Mestrado). https://teses.usp.br/teses/disponiveis/45/45133/tde-25062007-163150/publico/dissertacao_4.pdf

• Análise Multivariada. https://www.professores.uff.br/samuelcampos/analise-multivariada/

• Severn, K (2023) Multivariate Statistics. https://rich-d-wilkinson.github.io/MATH3030/

• Powell, J (2019) Multivariate statistics in R. https://www.hiercourse.com/multivariate

• Plotting PCA/clustering results using ggplot2 and ggfortify: https://rpubs.com/sinhrks/plot_pca

• ggfortify : Extension to ggplot2 to handle some popular packages - R software and data visualization: http://www.sthda.com/english/wiki/ggfortify-extension-to-ggplot2-to-handle-some-popular-packages-r-software-and-data-visualization

• A Little Book of Python for Multivariate Analysis: https://python-for-multivariate-analysis.readthedocs.io/index.html

Matriz e vetor no YouTube

• Zach Star: The Applications of Matrices | What I wish my teachers told me way earlier

• Zach Star: O determinante | A essência da Álgebra Linear, capítulo 5

• Zach Star: The applications of eigenvectors and eigenvalues | That thing you heard in Endgame has other uses