Prova 07

CONTEÚDO

Capítulo 16 - Regressão
Diagrama de Dispersão
Sumário da Regressão
Análise Estatística do Sumário da Regressão
Matemática envolvida no teste de Regressão
Testes Necessários para a Análise de Regressão
Teste de Shapiro-Wilk
Teste de Jarque-Bera
Adequação da Forma Funcional
Homoscedasticidade
Segundo Exemplo
Exercícios do Freund, 2006
Exercício 16.3
Exercício 16.7
Exercício 16.11
Exercício 16.13
Exercício 16.39
Exercício 16.47
Capítulo 17 - Correlação
Conceitos Fundamentais sobre a Correlação
Exemplo de Aplicação
A significância da correlação
Correlação parcial
Intervalo de confiança de correlações
Matriz de correlação
Exemplo de aplicação 01
Capítulo 18 - Testes não-paramétricos
Teste de Sinais
Exercício 18.1
Exercício 18.2
Exercício 18.5
Exercício 18.6
Exercício 18.8
Exercício 18.11
Exercício 18.12
Exercício 18.15
Exercícios da Lista (Julho/2017)
Exercício 01
Exercício 02
Exercício 03

Capítulo 16 - Regressão

Primeiro Exemplo

Aplique os conceitos de Regressão nos seguintes pares de dados:

x = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15)
y = c(5, 1, 4, 2, 10, 3, 5, 12, 9, 12, 6, 15, 22, 19, 21)
cbind(x, y)

##        x  y
##  [1,]  1  5
##  [2,]  2  1
##  [3,]  3  4
##  [4,]  4  2
##  [5,]  5 10
##  [6,]  6  3
##  [7,]  7  5
##  [8,]  8 12
##  [9,]  9  9
## [10,] 10 12
## [11,] 11  6
## [12,] 12 15
## [13,] 13 22
## [14,] 14 19
## [15,] 15 21

Diagrama de Dispersão

Uma regressão linear é a verificação estatística de que entre essas duas variáveis exista uma relação linear, ou seja, $y = β_0 + β_1x + Ꞓ$.Conhecidos os valores de yi e xi, devem ser calculados os valores de $β_0$, $β_1$ e $Ꞓ_i$. Na equação comentada neste parágrafo a variável do lado esquerdo (y) é chamada de variável dependente e a variável do lado direito (x) é chamada de variável independente. Geometricamente, calcula-se uma reta tal que a distância de cada $y_i$ a cada $y$ da reta minimiza a soma dos $Ꞓ^2$.A reta da figura abaixo (em azul) é denominada de reta dos mínimos quadrados, pois obedece à equação abaixo.

\[ min(\sum_{1=1}^n(Ꞓ_i^2)) \]

m = lm(y~x)
plot(x, y, xlim = c(0, max(x)), main="Gráfico de Dispersão", ylim = c(floor(m$coefficients[1]), max(y)))
abline((m = lm(y~x)), col = "blue")
for (i in 1:length(x)) lines(c(x[i], x[i]), c(y[i], m$fitted.values[i]), col = "red", lty = 3, lwd = 1.5)

Sumário da Regressão

options(scipen = 999)
(s=summary(lm(y~x)))

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.7190 -2.2405 -0.0476  2.1167  5.6238 
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  -0.8952     2.0258  -0.442     0.666    
## x             1.3286     0.2228   5.963 0.0000472 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.728 on 13 degrees of freedom
## Multiple R-squared:  0.7323, Adjusted R-squared:  0.7117 
## F-statistic: 35.56 on 1 and 13 DF,  p-value: 0.00004724

Análise Estatística do Sumário da Regressão

Iniciando a leitura pela última linha do resultado apresentado acima, iremos verificar se existe regressão ou não.

options(scipen = 999)
s$fstatistic

##    value    numdf    dendf 
## 35.55508  1.00000 13.00000

1-pf(s$fstatistic[1], s$fstatistic[2], s$fstatistic[3])

##         value 
## 0.00004724112

Esse teste F tem como hipótese nula a ausência de regressão, ou seja, não há uma relação linear entre x e y. Como a $H_0$ foi rejeitada, pode-se dizer que a regressão funcionou.

O que mais podemos afirmar sobre o teste de regressão utilizado? Podemos compreender melhor isso através do comando $summary$ realizado logo acima. Veja o comando a seguir:

c(s$r.squared, s$adj.r.squared)

## [1] 0.7322628 0.7116676

Onde há referência ao coeficiente $R^2$ , que representa a porção de variância da variável y que é explicada pela variável x. O segundo valor, $R^2$ ajustado, é o valor de $R^2$ ajustado aos graus de liberdade do sistema, sendo o valor comumente utilizado.

Os graus de liberdade do sistema dependem do número de observações e do número de variáveis que participam da relação. No presente caso, há apenas uma variável, mas podem existir mais; no caso, a regressão se chama regressão múltipla. O $R^2$ ajustado é penalizado em relação ao $R^2$, pelo fato de que o valor de $R^2$ aumenta à medida em que são acrescentadas variáveis independentes na regressão, mesmo que essas não guardem associação com a variável dependente.

Matemática envolvida no teste de Regressão

A solução das equações para uma regressão baseiam-se na equação fundamental $XB=Y$ onde $X$ é o vetor de observações da variável independente, acrescido de uma coluna igual a 1 para representar a constante de regressão; $B$ é o vetor de estimativas dos coeficientes de regressão e $Y$ é o vetor de observações da variável dependente.Se X fosse uma matriz quadrada, a solução seria imediata, mas como não é, torna-se necessária alguma álgebra matricial, onde $XB=Y$, $X'XB=X'Y$ e $B=(X'X)^-1$.

# número de observações
n = length(x)
# número de estimadores (exceto a constante)
k = 1
# matriz x aumentada
x1 = cbind(Intercept = 1, x = x)
# pré-multiplica pela sua transposta
tx1x = t(x1)%*%x1
# inverte o resultado
txt1m1 = solve(tx1x)
# calcula os coeficientes b
(b = (txt1m1%*%t(x1)%*%y))

##                 [,1]
## Intercept -0.8952381
## x          1.3285714

Para os demais cálculos, calcula-se os somatórios dos quadrados e seus respectivos graus de liberdade:

\[ SST=\sum_{1=1}^n(y_i-y)^2 \] \[ SSTr=\sum_{1=1}^n(ў_i-ӯ)^2 \] \[ SSTE=\sum_{1=1}^n(y_i-ў_i)^2 \]

# calcula os valores estimados
yhat = x1%*%b
ymean = mean(y)
SST = sum((y-ymean)^2)
SSTr = sum((yhat-ymean)^2)
SSE = sum((y-yhat)^2)
SST == (SSTr+SSE)

## [1] TRUE

## [1] TRUE
dfT = n-1
dfTr = k
dfE = n-k-1

Com isso calcula-se a estatística F:

(F = (SSTr/dfTr)/(SSE/dfE))

## [1] 35.55508

1-pf(F, dfTr, dfE)

## [1] 0.00004724112

Calcula-se os valores de $R^2$:

(R2 = 1-SSE/SST)

## [1] 0.7322628

(R2.adj = 1-(SSE/dfE)/(SST/dfT))

## [1] 0.7116676

Podemos calcular também o valor do erro-padrão dos resíduos:

sigma2 = SSE/dfE
sqrt(sigma2)

## [1] 3.72832

dfE

## [1] 13

Calculamos o erro padrão dos coeficientes:

Vb = sigma2*txt1m1
(Se = sqrt(diag(Vb)))

## Intercept         x 
## 2.0258128 0.2228097

E em seguida calculamos as significâncias dos coeficientes:

(t = b/Se)

##                 [,1]
## Intercept -0.4419155
## x          5.9628077

(p = 2*(1-pt(abs(t), n-2)))

##                    [,1]
## Intercept 0.66581247322
## x         0.00004724112

Colocando tudo isso em uma matriz temos:

data.frame('b' = b, Se = Se, t = t, p = p)

##                    b        Se          t             p
## Intercept -0.8952381 2.0258128 -0.4419155 0.66581247322
## x          1.3285714 0.2228097  5.9628077 0.00004724112

Esse procedimento pode ser expandido para quantas variáveis independentes for necessário, mas com certeza é bem mais fácil usar a função $lm$ do R.

Testes Necessários para a Análise de Regressão

Normalidade dos resíduos

O teste de normalidade garante que os resíduos da regressão $(y_i − ў_i)$ sigam uma distribuição normal.Qualquer teste de normalidade é adequado. Costuma-se usar o teste de Shapiro-Wilk. No exemplo apresentado, observa-se a forma dos resíduos.

hist(m$residuals, freq = FALSE, xlim = c(-8, 8), ylim = c(0, 0.12), main="Histograma de Resíduos")
curve(dnorm(x, mean(m$residuals), sd(m$residuals)), xlim = c(-8, 8), add = TRUE, col = "blue")

Teste de Shapiro-Wilk

shapiro.test(m$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  m$residuals
## W = 0.97913, p-value = 0.9632

Como o valor p do teste é maior que 0.05, aceita-se a H0 de que os resíduos seguem uma distribuição normal.Quando o número de observações é pequeno, os testes de normalidade geralmente não possuem potência suficiente para rejeitar a hipótese nula, resultando em um aumento do Erro Tipo II.

Teste de Jarque-Bera

O Teste Jarque-Bera, um tipo de teste de normalidade. O teste de Jarque-Bera é normalmente usado para grandes conjuntos de dados, porque outros testes de normalidade não são confiáveis quando n é grande (por exemplo, Shapiro-Wilk não é confiável com n mais de 2.000). O teste de Jarque-Bera tem como hipótese nula a normalidade. Assim, se o p-valor for menor do que 5% (ou 10%), p<0,05 (p<0,10), então o autor rejeita a normalidade. Já se p>0,05, aceita-se a normalidade.

library(moments)
jarque.test(m$residuals)

## 
##  Jarque-Bera Normality Test
## 
## data:  m$residuals
## JB = 0.36963, p-value = 0.8313
## alternative hypothesis: greater

library(DescTools)
JarqueBeraTest(m$residuals, robust = FALSE)

## 
##  Jarque Bera Test
## 
## data:  m$residuals
## X-squared = 0.36963, df = 2, p-value = 0.8313

**Da mesma forma, aceita-se a $H_0$ de normalidade do resíduos.

Adequação da Forma Funcional

A linearidade da relação funcional entre as variáveis x e y é atestada pelo teste RESET de Ramsey, aplicado ao modelo de regressão. Esta função está na biblioteca lmtest.

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(zoo)
m = lm(y~x)
resettest(m)

## 
##  RESET test
## 
## data:  m
## RESET = 1.5879, df1 = 2, df2 = 11, p-value = 0.2478

O valor $p ≥ 0.05$ atesta que a forma funcional está adequada, ou seja, que a relação entre x e y é linear - ou seja, possui a forma $y = β0 + β1x + Ꞓ$.

Um contra-exemplo

Suponha a massa de dados gerada a seguir, consistindo de 300 observações de $x_a$ e $y_a$.

set.seed(12)
xa = 1:300
ya = 30*xa+(xa/2)^2+1800*rnorm(length(xa))
m1 = lm(ya~xa)
m2 = lm(ya~xa+I(xa^2))
yahat2 = cbind(1, xa, xa^2)%*%m2$coefficients

No gráfico a seguir observa-se o ajuste de uma reta $(y = b_0 + b_1x + Ꞓ)$ e de uma curva de segundo grau $(y = b_0 + b_1x + b_2x2 + Ꞓ)$.

plot(xa, ya, cex = 0.2)
abline(m1, col = "red")
lines(xa, yahat2, col = "blue", lwd = 2)

O desajuste da reta obtida pela equação do primeiro grau não é evidente pelo sumário da regressão que indica, inclusive, um bom ajuste.

summary(m1)

## 
## Call:
## lm(formula = ya ~ xa)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6588.5 -1592.1  -124.7  1599.9  6526.8 
## 
## Coefficients:
##              Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) -3837.699    273.564  -14.03 <0.0000000000000002 ***
## xa            105.377      1.575   66.89 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2363 on 298 degrees of freedom
## Multiple R-squared:  0.9375, Adjusted R-squared:  0.9373 
## F-statistic:  4474 on 1 and 298 DF,  p-value: < 0.00000000000000022

A aplicação do teste RESET, todavia, indica a inadequação do modelo linear.

resettest(m1)

## 
##  RESET test
## 
## data:  m1
## RESET = 131.02, df1 = 2, df2 = 296, p-value < 0.00000000000000022

Com um p < 0.05 o teste rejeita a H0 de que a forma funcional aplicada é adequada. Modificando o modelo para um modelo de segundo grau, observa-se a adequação:

resettest(m2)

## 
##  RESET test
## 
## data:  m2
## RESET = 1.7946, df1 = 2, df2 = 295, p-value = 0.168

summary(m2)

## 
## Call:
## lm(formula = ya ~ xa + I(xa^2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5399.4 -1115.2   -31.4  1188.8  3839.5 
## 
## Coefficients:
##               Estimate Std. Error t value             Pr(>|t|)    
## (Intercept) -197.91516  301.01416  -0.657                0.511    
## xa            33.06383    4.61806   7.160     0.00000000000638 ***
## I(xa^2)        0.24024    0.01486  16.169 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1726 on 297 degrees of freedom
## Multiple R-squared:  0.9668, Adjusted R-squared:  0.9666 
## F-statistic:  4322 on 2 and 297 DF,  p-value: < 0.00000000000000022

Homoscedasticidade

Homocedasticidade é o termo para designar variância constante dos erros experimentais $\varepsilon_i$ para observações distintas. Caso a suposição de homocedasticidade não seja válida, podemos listar alguns efeitos no ajuste do modelo, como: Os erros padrões dos estimadores, obtidos pelo Método dos Mínimos Quadrados, são incorretos e portanto a inferência estatística não é válida, e que não podemos mais dizer que os Estimadores de Mínimos Quadrados são os melhores estimadores de mínima variância para β, embora ainda possam ser não viciados.A condição de homoscedasticidade é obtida quando a variância da amostra é constante ao longo das observações. Segue um exemplo de dados homoscedásticos.

set.seed(14)
xa = 1:100
ya1 = xa+rnorm(length(xa), 0, 10)
plot(xa, ya1)

A seguir tem-se um exemplo de dados heteroscedásticos (ou não homoscedásticos), ou seja, cuja variância é diferente ao longo da amostra. No caso, observa-se uma variância maior no final da amostra que no início dela.

set.seed(14)
xa = 1:100
ya2 = xa+rnorm(length(xa), 0, 10+xa/3)
plot(xa, ya2)

Existem alguns testes para a heteroscedasticidade e, para nossos estudos, usaremos o teste de Breush-Pagan, também disponível na biblioteca lmtest.

bptest(m)

## 
##  studentized Breusch-Pagan test
## 
## data:  m
## BP = 0.15452, df = 1, p-value = 0.6943

O valor p > 0.05 atesta a homoscedasticidade do modelo de regressão. A título de contra-prova, segue a análise de homoscedasticidade aplicada ao modelo linear para os dados $x_a$ e $y_a2$

bptest(lm(ya2~xa))

## 
##  studentized Breusch-Pagan test
## 
## data:  lm(ya2 ~ xa)
## BP = 11.506, df = 1, p-value = 0.0006936

Rejeitando a hipótese nula de homoscedasticidade.

Segundo Exemplo

Um empregador trabalha com a hipótese de que o número de erros de montagem cresce linearmente com relação ao número de dias que uma equipe trabalha ininterruptamente em regime de turno. Ele afirma que essa relação é menor que 0.2 erros/dia. Considerando os dados coletados apresentados a seguir, teste essa hipótese:

x = 1:23
y = c(1,2,2,3,4,3,1,7,2,6,4,4,5,6,7,8,6,9,8,5,11,7,7)
cbind(x,y)

##        x  y
##  [1,]  1  1
##  [2,]  2  2
##  [3,]  3  2
##  [4,]  4  3
##  [5,]  5  4
##  [6,]  6  3
##  [7,]  7  1
##  [8,]  8  7
##  [9,]  9  2
## [10,] 10  6
## [11,] 11  4
## [12,] 12  4
## [13,] 13  5
## [14,] 14  6
## [15,] 15  7
## [16,] 16  8
## [17,] 17  6
## [18,] 18  9
## [19,] 19  8
## [20,] 20  5
## [21,] 21 11
## [22,] 22  7
## [23,] 23  7

plot(x,y)

Pressupondo uma relação linear, monta-se o modelo e ele é testado quanto à sua adequação:

(m=lm(y~x))

## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##      1.3123       0.3182

shapiro.test(m$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  m$residuals
## W = 0.97758, p-value = 0.8612

resettest(m)

## 
##  RESET test
## 
## data:  m
## RESET = 0.45586, df1 = 2, df2 = 19, p-value = 0.6407

bptest(m)

## 
##  studentized Breusch-Pagan test
## 
## data:  m
## BP = 1.6574, df = 1, p-value = 0.198

Nada havendo que desabone o modelo, passa-se à análise dos coeficientes:

(s=summary(m))

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6759 -0.9713 -0.2213  1.0059  3.1423 
## 
## Coefficients:
##             Estimate Std. Error t value   Pr(>|t|)    
## (Intercept)   1.3122     0.7048   1.862     0.0767 .  
## x             0.3182     0.0514   6.190 0.00000385 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.635 on 21 degrees of freedom
## Multiple R-squared:  0.646,  Adjusted R-squared:  0.6291 
## F-statistic: 38.31 on 1 and 21 DF,  p-value: 0.000003853

Para testar a hipótese do empregador $H_0 : β ≤ 0.2$, e elabora-se a hipótese alternativa $H_1 : β > 0.2$. Na regressão encontramos um valor de b1 maior que 0.2, então torna-se necessário verificar se ele é suficientemente maior para rejeitar a hipótese nula. Um teste t resolve isso:

\[ t= \frac{ḇ-b}{S_{eb}} \]

onde ḇ é a estimativa de b obtida a partir dos dados (0.3182), b é o valor proposto pelo empregador(0.2) e $Se_b$ é o erro-padrão de ḇ (0.0514).

(b1 = s$coefficients[2,])

##       Estimate     Std. Error        t value       Pr(>|t|) 
## 0.318181818182 0.051403302983 6.189909980807 0.000003852972

t = (b1[1]-0.2)/b1[2]
1-pt(t, length(x)-2)

##   Estimate 
## 0.01593082

Com isso, a 95% de confiança, rejeita-se a H0 proposta pelo empregador.

Exercícios do Freund, 2006

Exercício 16.3

Para verificar se um conservante de alimentos largamente utilizado contribui para a hiperatividade de crianças em idade pré-escolar, um nutricionista escolheu uma amostra aleatória de dez crianças de quatro anos reconhecidas como bastante hiperativas de várias escolinhas e observou seu comportamento 45 minutos depois de terem ingerido quantidades controladas de comida contendo conservante. Na tabela a seguir, é a quantidade de comida consumida contendo o conservante (em gramas) e é uma medição subjetiva de hiperatividade (numa escala de 1 a 20) baseada na agitação da criança e na interação com outras crianças:

x = c(36, 82, 45, 49, 21, 24, 58, 73, 85, 52)
y = c(6, 14, 5, 13, 5, 8, 14, 11, 18, 6)

m = lm(y~x)
shapiro.test(m$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  m$residuals
## W = 0.89929, p-value = 0.2152

resettest(m)

## 
##  RESET test
## 
## data:  m
## RESET = 0.23535, df1 = 2, df2 = 6, p-value = 0.7973

bptest(m)

## 
##  studentized Breusch-Pagan test
## 
## data:  m
## BP = 0.015088, df = 1, p-value = 0.9022

Feita a regressão, nada foi encontrado que coloque a validade do modelo em dúvida. Passa-se à análise.

(s = summary(m))

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9193 -2.0659 -0.3386  2.7156  3.5650 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  1.52570    2.56461   0.595  0.56835   
## x            0.16142    0.04528   3.565  0.00735 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.045 on 8 degrees of freedom
## Multiple R-squared:  0.6137, Adjusted R-squared:  0.5654 
## F-statistic: 12.71 on 1 and 8 DF,  p-value: 0.007346

Constata-se a existência de regressão, explicando aproximadamente 56.54% da variância da hiperatividade, por meio de uma relação significante e positiva (b1 = 0.1614) com a ingestão de alimentos com conservante.

Exercício 16.7

Os dados a seguir mostram os tempos médios semanais, em horas, que seis estudantes dedicaram aos seus trabalhos para case os índices de pontuação para suas disciplinas que fizeram naquele semestre.Admitindo que todas as suposições subjacentes à análise de regressão normal tenham sido satisfeitas, construa um intervalo de 99% de confiança para β, a quantidade pela qual um estudante da população amostrada poderia aumentar seu índice de pontuação estudando uma hora extra por semana.

x = c(15, 28, 13, 20, 4, 10)       # horas gastas em deveres de casa
y = c(20, 27, 13, 19, 09, 17)/10   # Índice de pontuação
(s = summary((m=lm(y~x))))

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6 
##  0.25000  0.05814 -0.31279 -0.19302 -0.09535  0.29302 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.72093    0.24641   2.926  0.04300 * 
## x            0.06860    0.01467   4.678  0.00946 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.272 on 4 degrees of freedom
## Multiple R-squared:  0.8455, Adjusted R-squared:  0.8068 
## F-statistic: 21.88 on 1 and 4 DF,  p-value: 0.009461

alfa=0.01
(limb1 = s$coefficients[2,1]+qt(1-alfa/2, length(x)-2)*s$coefficients[2,2]*c(-1, 1))

## [1] 0.001085074 0.136124228

Comentário: Qual a quantidade que o estudante poderia aumentar seu índice de pontuação? …

Exercício 16.11

Os dados referentes aos lucros líquidos (em milhares de unidades monetárias) de uma companhia durante os seis primeiros anos de operação são os seguintes:

ano = 1:6
ll = c(112, 149, 238, 354, 580, 867)
par(mfrow=c(1, 2))  # divide graph area in 2 columns
plot(ano, ll, type = "b", xlab="Anos em Operação", ylab="Lucro Líquido")
plot(ano, log10(ll), type = "b", xlab="Anos em Operação", ylab="Lucro Líquido")

options(scipen = 999)
(s = summary((m=lm(log10(ll)~ano))))

## 
## Call:
## lm(formula = log10(ll) ~ ano)
## 
## Residuals:
##         1         2         3         4         5         6 
##  0.030538 -0.027984 -0.007083 -0.017147  0.014787  0.006888 
## 
## Coefficients:
##             Estimate Std. Error t value    Pr(>|t|)    
## (Intercept) 1.836190   0.022449   81.79 0.000000134 ***
## ano         0.182490   0.005764   31.66 0.000005934 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02411 on 4 degrees of freedom
## Multiple R-squared:  0.996,  Adjusted R-squared:  0.995 
## F-statistic:  1002 on 1 and 4 DF,  p-value: 0.000005934

(a = 10^s$coeff[1,1])

## [1] 68.57875

(b = 10^s$coeff[2,1])

## [1] 1.522265

Os valores de a e b resultantes foram respectivamente 68.57875 e 1.522265 para a equação $y=a*b^x$.

Exercício 16.13

O tempo de secagem (em horas) de um verniz e a quantidade (em gramas) de um certo aditivo químico são os seguintes:

x = 1:8
y = c(72, 67, 47, 37, 47, 42, 52, 57)/10

Ajuste os dados a uma parábola.

Faça uma previsão dos tempos de secagem quando se adiciona 6.5 gramas de aditivo químico.

plot(x,y, xlab="Quantidade de Aditivo (gr)", ylab="Tempo de Secagem (hrs)")

(s = summary((m = lm(y~x+I(x^2)))))

## 
## Call:
## lm(formula = y ~ x + I(x^2))
## 
## Residuals:
##       1       2       3       4       5       6       7       8 
## -0.2292  0.6875 -0.2946 -0.6756  0.5446 -0.1339  0.2887 -0.1875 
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  9.24464    0.76453  12.092 0.0000683 ***
## x           -2.01488    0.38979  -5.169   0.00356 ** 
## I(x^2)       0.19940    0.04228   4.716   0.00526 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.548 on 5 degrees of freedom
## Multiple R-squared:  0.8531, Adjusted R-squared:  0.7943 
## F-statistic: 14.51 on 2 and 5 DF,  p-value: 0.008276

yhat = predict(m)
plot(x, y, ylim = c(min(c(y, yhat)), max(c(y, yhat))),xlab="Quantidade de Aditivo (gr)", ylab="Tempo de Secagem (hrs)")
points(x, yhat, pch = 8, col = "#F02020")
lines(x, yhat, pch = 8, col = "#F02020", lty = 3, lwd = 1.5)

newx = data.frame(x = c(6.5))
predict(m, newx)

##        1 
## 4.572768

A previsão do tempo de secagem quando se adiciona 6.5 gramas de aditivo químico resultou em 4.572768 hrs.

Exercício 16.39

Os dados a seguir referem-se a uma colônia de bactérias num meio de cultura.

x = c(1:5*2)
y = c(112, 148, 241, 363, 585)

Monte as duas equações para normais para o ajuste de uma curva exponencial
Transforme a equação obtida em 1 no tipo $y=a.b^x$
Use a equação obtida em 2 para estimar a contagem de bactérias cinco dias depois da inoculação.

par(mfrow=c(1, 2))  # divide graph area in 2 columns
plot(x,y, xlab="Dias desde a inoculação", ylab="Contagem de bactérias (milhares")
plot(x, log10(y), xlab="Dias desde a inoculação", ylab="Contagem de bactérias (milhares")

(s = summary((m = lm(log10(y)~x))))

## 
## Call:
## lm(formula = log10(y) ~ x)
## 
## Residuals:
##         1         2         3         4         5 
##  0.028610 -0.032898 -0.003695 -0.008357  0.016340 
## 
## Coefficients:
##             Estimate Std. Error t value   Pr(>|t|)    
## (Intercept) 1.838056   0.028731   63.97 0.00000842 ***
## x           0.091276   0.004331   21.07   0.000234 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02739 on 3 degrees of freedom
## Multiple R-squared:  0.9933, Adjusted R-squared:  0.9911 
## F-statistic: 444.1 on 1 and 3 DF,  p-value: 0.0002338

\[ log_{10}(y)=log_{10}(a)+log_{10}(b).x \] \[ log_{10}(y)=1.8381+0.0913.x \]

(a = 10^m$coeff[1])

## (Intercept) 
##    68.87406

(b = 10^m$coeff[2])

##        x 
## 1.233889

\[ y=a.b^x \] \[ y=68.8741 . (1.2339)^x \]

(yp = a*b^5)

## (Intercept) 
##    196.9862

#ou

newx = data.frame(x = c(5))
10^predict(m, newx)

##        1 
## 196.9862

Exercício 16.47

Os dados a seguir se referem à demanda y (em milhares de unidades) por um produto x (em unidades monetárias) em cinco áreas comerciais bastante semelhantes. Baseado nestes dados, qual será a demanda quando o preço for de 12 unidades monetárias? Use o R ou outro aplicativo para ajustar uma parábola a esses dados, além de estimar a demanda quando o preço do produto é de 12 unidades monetárias.

x = c(20, 16, 10, 11, 14)
y = c(22, 41, 120, 89, 56)

# verifica-se o formato da curva
plot(x, y, xlab="Preço", ylab="Demanda")
# aproximação de primeiro grau
(s = summary((m = lm(y~x))))

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##       1       2       3       4       5 
##   9.978  -7.972  15.602  -6.160 -11.448 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  196.775     25.206   7.807  0.00438 **
## x             -9.238      1.721  -5.369  0.01265 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.85 on 3 degrees of freedom
## Multiple R-squared:  0.9057, Adjusted R-squared:  0.8743 
## F-statistic: 28.82 on 1 and 3 DF,  p-value: 0.01265

mlinear = m
plot(x, y, xlab="Preço", ylab="Demanda"); abline(m, col = "red")

Fazendo uma aproximação com segundo grau tem-se:

(s = summary((m = lm(y~x+I(x^2)))))

## 
## Call:
## lm(formula = y ~ x + I(x^2))
## 
## Residuals:
##       1       2       3       4       5 
## -1.0111  3.0837  5.9404 -7.8869 -0.1261 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 384.3925    65.1637   5.899   0.0276 *
## x           -35.9975     9.1421  -3.938   0.0589 .
## I(x^2)        0.8964     0.3047   2.942   0.0987 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.35 on 2 degrees of freedom
## Multiple R-squared:  0.9823, Adjusted R-squared:  0.9646 
## F-statistic: 55.51 on 2 and 2 DF,  p-value: 0.0177

plot(x, y, xlab="Preço", ylab="Demanda")
x1 = 100:200/10
yhat = predict(m, data.frame(x = x1))
lines(x1, yhat, col = "blue")
abline(mlinear, col = "red")

Estimando a demanda quando o preço do produto é de 12 unidades monetárias.

predict(m, data.frame(x = 12))

##        1 
## 81.50713

Capítulo 17 - Correlação

Este capítulo aborda a correlação chamada de produto-momento de Pearson, que é a correlação paramétrica. Para que se possa usá-la, os dados devem estar em uma escala no mínimo intervalar e possuírem distribuição normal.Observa-se que as variáveis x e y devem possuir o mesmo número de observações.

Conceitos Fundamentais sobre a Correlação

Através do gráfico de dispersão pode-se indicar se a correlação linear é positiva, negativa ou a inexistência de correlação. Além disso, podemos identificar, por meio do coeficiente de correlação (r), o grau da correlação.Por exemplo, para r = +0,953, temos uma correlação positiva forte, ou seja, a variável dependente y cresce quase na mesma proporção que a variável independente x - Para avaliar se esse resultado é signiﬁcativo, pode-se realizar um Teste de Hipóteses para a o Coeﬁciente de Correlação.

Equação (17.1) \[ S_{xx}=\sum_{1=1}^n(x_i-ẍ)^2 \]

Equação (17.2)

\[S_{yy}=\sum_{1=1}^n(y_i-ӯ)^2\]

Equação (17.3)

\[S_{xy}=\sum_{1=1}^n(x_i-ẍ)(y_i-ӯ)\]

Equação (17.4)

\[cor(x,y)=\frac{S{xy}}{\sqrt{S_{xx}S_{yy}}}\]

Exemplo de Aplicação

No exemplo que consta no livro de Freund na página 435, são apresentadas notas que 12 estudantes obtidas em exames finais de Economia e Antropologia. Para se calcular o Coeficiente de correlação do momento do produto (r) temos que:

econom = c(51, 68, 72, 97, 55, 73, 95, 74, 20, 91, 74, 80)
antrop = c(74, 70, 88, 93, 67, 73, 99, 73, 33, 91, 80, 86)
cbind(econom, antrop)

##       econom antrop
##  [1,]     51     74
##  [2,]     68     70
##  [3,]     72     88
##  [4,]     97     93
##  [5,]     55     67
##  [6,]     73     73
##  [7,]     95     99
##  [8,]     74     73
##  [9,]     20     33
## [10,]     91     91
## [11,]     74     80
## [12,]     80     86

Calculando de forma manual no R temos:

Sxx = sum((econom-mean(econom))^2)
Syy = sum((antrop-mean(antrop))^2)
Sxy = sum((econom-mean(econom))*(antrop-mean(antrop)))
c(Sxy, Sxx, Syy)

## [1] 3790.500 5021.667 3272.250

r = Sxy/sqrt(Sxx*Syy)
r

## [1] 0.9350812

Esse procedimento pode ser implementado no R por meio da função cor:

cor(econom, antrop)

## [1] 0.9350812

A significância da correlação

Para que se considere a correlação obtida como diferente de zero, é necessário testar a sua significância, que é dada pela equação abaixo, onde t é uma variável que segue a distribuição de Student com $n − 2$ graus de liberdade.Veja o código a seguir:

$$ t=r{}

options(scipen = 999)
n = length(antrop)
t = r*sqrt((n-2)/(1-r^2))
t

## [1] 8.34285

p = 2*(1-pt(t, n-2))
p

## [1] 0.000008138117

O procedimento é implementado no R por meio da função $cor.test$:

cor.test(econom, antrop)

## 
##  Pearson's product-moment correlation
## 
## data:  econom and antrop
## t = 8.3429, df = 10, p-value = 0.000008138
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7794872 0.9819986
## sample estimates:
##       cor 
## 0.9350812

É importante observar que, a partir da equação 17.4, $cor(x, y) = cor(y, x)$.

Correlação parcial

Em teoria das probabilidades e estatística, a correlação parcial mede o grau de associação entre duas variáveis aleatórias, excluindo o efeito de uma terceira variável.Em certos casos, o coeficientee de correlação linear simples pode nos levar a equvocos na interpretação da associação entre duas variáveis, pois este não considera a influência das demais variáveis contidas no conjunto de dados. O coeficiente de correlação parcial é uma técnica baseada em operações matriciais que nos permite identificar a associação entre duas variáveis retirando-se os efeitos das demais variáveis presentes.

Considere os seguintes dados:

x1 = número de taças de chocolate quente vendidas

x2 = número de pessoas que frequentam a estação de veraneio

x3 = temperatura no momento da observação

r12 = -0.30
r13 = -0.70
r23 = 0.80
r123 = (r12-r13*r23)/(sqrt(1-r13^2)*sqrt(1-r23^2))
r123

## [1] 0.606788

Usando o R temos:

pcor = function(r12, r13, r23) {
r123 = (r12-r13*r23)/(sqrt(1-r13^2)*sqrt(1-r23^2))
r123
}
pcor(r12, r13, r23)

## [1] 0.606788

Ou, quando de posse dos dados brutos, pode-se usar uma das funções das bibliotecas do R, por exemplo, a função CorPart da biblioteca MASS ou a função partial.r da biblioteca psych:

m = matrix(c(1, r12, r13, r12, 1, r23, r13, r23, 1), nrow = 3)
m

##      [,1] [,2] [,3]
## [1,]  1.0 -0.3 -0.7
## [2,] -0.3  1.0  0.8
## [3,] -0.7  0.8  1.0

library(DescTools)
library(MASS)
x = mvrnorm(300, c(0, 0, 0), m, empirical = TRUE)
CorPart(m = m, 1:2, 3)

##          [,1]     [,2]
## [1,] 1.000000 0.606788
## [2,] 0.606788 1.000000

library(psych)

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:DescTools':
## 
##     AUC, ICC, SD

(a = psych::partial.r(data = m, x = 1:2, y = 3))

## partial correlations 
##      [,1] [,2]
## [1,] 1.00 0.61
## [2,] 0.61 1.00

a[1,2]

## [1] 0.606788

Outra forma é correlacionar os resíduos de duas regressões lineares:

lm1 = lm(x[,1]~x[,3])
lm2 = lm(x[,2]~x[,3])
cor(lm1$res, lm2$res)

## [1] 0.606788

# ou fazer uma regressão sobre os resíduos padronizados, onde o coeficiente angular é a correlação parcial

lm(scale(lm1$res)~I(scale(lm2$res)))$coeff[2]

## I(scale(lm2$res)) 
##          0.606788

Intervalo de confiança de correlações

Em estatística, um intervalo de confiança (CI) é um tipo de estimativa de intervalo, calculado a partir das estatísticas dos dados observados, que pode conter o valor real de um parâmetro populacional desconhecido. O intervalo tem um nível de confiança associado que, falando de maneira simples, quantifica o nível de confiança de que o parâmetro está no intervalo. Mais estritamente falando, o nível de confiança representa a frequência (ou seja, a proporção) de possíveis intervalos de confiança que contêm o verdadeiro valor do parâmetro da população desconhecida. Em outras palavras, se os intervalos de confiança forem construídos usando um determinado nível de confiança de um número infinito de estatísticas de amostra independentes, a proporção desses intervalos que contêm o valor real do parâmetro será igual ao nível de confiança.

Inicialmente iremos utilizar o método de Fisher:

\[ z=\frac{1}{2}log\frac{1+r}{1-r} \]

r = 0.5
(z = 1/2*log((1+r)/(1-r)))

## [1] 0.5493061

# ou 
atanh(r)

## [1] 0.5493061

# ou 
FisherZ(0.5)

## [1] 0.5493061

E sua inversa:

\[ r=\frac{e^{2z}-1}{e^{2z}+1} \]

(r = (exp(2*z)-1)/(exp(2*z)+1))

## [1] 0.5

# ou
tanh(z)

## [1] 0.5

# ou 
FisherZInv(z)

## [1] 0.5

Sabendo que o erro-padrão de Z é:

\[ \sigma_Z=\frac {1}{\sqrt{n-3}} \]

pode-se construir um intervalo de confiança paramétrico seguindo os passos:

Calcular o Z associado a r
Calcular o σZ com base em n
Calcular o quantil da distribuição normal para a significância desejada
Calcular o intervalo de confiança de Z
Calcular os valores de r associados ao intervalo de confiança de Z

Por exemplo, supondo que se queira construir um intervalo de confiança de 95% para uma correlação r = 0.5 com n = 30.

r = 0.5
(Z = atanh(r))

## [1] 0.5493061

n=30
(sigmaZ = 1/sqrt(n-3))

## [1] 0.1924501

alfa = 0.05
(z = qnorm(1-alfa/2))

## [1] 1.959964

(E = z*sigmaZ)

## [1] 0.3771952

(intZ = Z+E*c(-1,1))

## [1] 0.1721109 0.9265014

tanh(intZ)

## [1] 0.1704314 0.7289586

Empacotando isso em uma função no R temos:

corIntConf = function(r, n, alfa = 0.05) {
Z = atanh(r)
sigmaZ = 1/sqrt(n-3)
z = qnorm(1-alfa/2)
E = z*sigmaZ
intZ = Z+E*c(-1,1)
tanh(intZ)
}

corIntConf(0.5, 30)

## [1] 0.1704314 0.7289586

Ou usando diretamente a função r.con da biblioteca psych tem-se que:

library(psych)
r.con(r, n)

## [1] 0.1704314 0.7289586

Ou ainda, usando a função $cor.test$, se houver acesso às observações dos dados:

(m = cbind(c(1, r), c(r, 1)))

##      [,1] [,2]
## [1,]  1.0  0.5
## [2,]  0.5  1.0

library(MASS)
x = mvrnorm(n, c(0, 0), m, empirical = TRUE)
cor.test(x[,1], x[,2])

## 
##  Pearson's product-moment correlation
## 
## data:  x[, 1] and x[, 2]
## t = 3.0551, df = 28, p-value = 0.0049
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1704314 0.7289586
## sample estimates:
## cor 
## 0.5

Matriz de correlação

Uma matriz de correlação é uma tabela de coeficientes de correlação para um conjunto de variáveis usadas para determinar se existe uma relação entre as variáveis. O coeficiente indica tanto a força do relacionamento quanto a direção (correlações positivas versus negativas).Existem diferentes métodos para análise de correlação: teste de correlação paramétrica de Pearson, análise de correlação baseada em classificação de Spearman e Kendall.

Para determinar se uma matriz de correlação é significativamente diferente de uma matriz identidade (Uma matriz identidade é resultante do cálculo de correlação entre variáveis ortogonais, ou seja, cujas correlações par a par são iguais a zero), usa-se o teste de Bartlett:

\[ X^2=log|R|(N-1-\frac{2p+5}{6}) \]

onde |R| é o determinante da matriz de correlações, N é o número de observações e p é o número de variáveis envolvidas.

N = 80
sigma1 = matrix(c(1, 0.2, 0.4, 0.2, 1, 0.5, 0.4, 0.5, 1), nrow = 3)
sigma2 = matrix(c(1, 0.3, 0.45, 0.3, 1, 0.7, 0.45, 0.7, 1), nrow = 3)
sigma3 = sigma1-sigma2
diag(sigma3) = 1
sigma3

##       [,1] [,2]  [,3]
## [1,]  1.00 -0.1 -0.05
## [2,] -0.10  1.0 -0.20
## [3,] -0.05 -0.2  1.00

(p = nrow(sigma1))

## [1] 3

(det1 = det(sigma1))

## [1] 0.63

(det2 = det(sigma2))

## [1] 0.4065

(det3 = det(sigma3))

## [1] 0.9455

(chisq = -log(abs(det3))*(N-1-(2*p+5)/6))

## [1] 4.324527

1-pchisq(chisq, p*(p-1)/2)

## [1] 0.2284863

library(psych)
psych::cortest.bartlett(R = sigma3, n = N)

## $chisq
## [1] 4.324527
## 
## $p.value
## [1] 0.2284863
## 
## $df
## [1] 3

Uma variante do teste de Bartlett é calculada por meio da transformada r para Z de Fisher.

\[ Rd = R_1 − R_2 \]

\[ Z = atanh (R_d) \]

\[ X_2=(\sum_{i=1}^n\sum_{j=1}^n(z_{ij}^2))\frac{N-3}{2} \]

E obtendo a probabilidade associada ao valor de $χ_2$ obtido com $p(p − 1)/2$ graus de liberdade.

zsigma3 = sigma3
diag(zsigma3) = 0
zsigma3 = atanh(zsigma3)
zsigma3 = sum(zsigma3^2)
(chisq = zsigma3*(N-3)/2)

## [1] 4.132732

1-pchisq(chisq, p*(p-1)/2)

## [1] 0.2474836

cortest(R1 = sigma1, R2 = sigma2, n1 = N, n2 = N)

## Tests of correlation matrices 
## Call:cortest(R1 = sigma1, R2 = sigma2, n1 = N, n2 = N)
##  Chi Square value 4.3  with df =  3   with probability < 0.23 
## z of differences =  0.02

cortest(R1 = sigma3, n1 = N)

## Tests of correlation matrices 
## Call:cortest(R1 = sigma3, n1 = N)
##  Chi Square value 4.13  with df =  3   with probability < 0.25

psych::cortest.normal(R1 = sigma1, R2 = sigma2, n1 = N, n2 = N)

## Tests of correlation matrices 
## Call:psych::cortest.normal(R1 = sigma1, R2 = sigma2, n1 = N, n2 = N)
##  Chi Square value 4.3  with df =  3   with probability < 0.23

psych::cortest.jennrich(R1 = sigma1, R2 = sigma2, n1 = N, n2 = N)

## $chi2
## [1] 3.375026
## 
## $prob
## [1] 0.3373361

psych::cortest.mat(R1 = sigma1, R2 = sigma2, n1 = N, n2 = N)

## Tests of correlation matrices 
## Call:psych::cortest.mat(R1 = sigma1, R2 = sigma2, n1 = N, n2 = N)
##  Chi Square value 11.56  with df =  6   with probability < 0.072

Exemplo de aplicação 01

Em dois testes de campo, a correlação entre duas variáveis foi medida como r1 = 0.7 no primeiro teste, com n1 = 44 observações, e como r2 = 0.3 no segundo teste, com n2 = 66 observações. É possível afirmar, com uma confiança de 95%, que essas correlações são iguais?

A solução para esse problema é construir os intervalos de confiança das duas variáveis e verificar se eles estão sobrepostos. Se não estiverem, a hipótese nula de igualdade das correlações deve ser rejeitada na significância especificada.

r1 = 0.7
n1 = 40
r2 = 0.3
n2 = 66
corIntConf(r1, n1)

## [1] 0.4968270 0.8304289

corIntConf(r2,n2)

## [1] 0.06250576 0.50534011

Foi por pouco, mas aceitou-se a H0 de igualdade das correlações.Pode-se definir uma função para verificar se os intervalos de confiança estão sobrepostos:

corEqual = function(r1, r2, n1, n2, alfa = 0.05, plt = FALSE, xlim = NULL) {
equal = FALSE
ic1 = corIntConf(r1, n1, alfa)
ic2 = corIntConf(r2, n2, alfa)
if(plt) {
if(is.null(xlim)) xlim = c(min(c(ic1, ic2)), max(c(ic1, ic2)))
plot(0, 0, xlim = xlim, ylim = c(0, 2), type = "n")
lines(ic1, c(1, 1), col = "blue", lwd = 3)
lines(ic2, c(1.5, 1.5), col = "red", lwd = 3)
abline(v = ic1, col = "blue", lty = 2)
abline(v = ic2, col = "red", lty = 2)
}
if((ic1[1]>ic2[1])&(ic1[1]<ic2[2])) equal = TRUE
if((ic1[2]>ic2[1])&(ic1[2]<ic2[2])) equal = TRUE
return(invisible(list(equal = equal, ic1 = ic1, ic2 = ic2)))
}

(corEqual(r1, r2, n1, n2, alfa = 0.05, plt = TRUE, xlim = c(0, 1)))

## $equal
## [1] TRUE
## 
## $ic1
## [1] 0.4968270 0.8304289
## 
## $ic2
## [1] 0.06250576 0.50534011

Para examinar melhor a sobreposição dos intervalos:

(corEqual(r1, r2, n1, n2, alfa = 0.05, plt = TRUE, xlim = c(0.48, 0.52)))

## $equal
## [1] TRUE
## 
## $ic1
## [1] 0.4968270 0.8304289
## 
## $ic2
## [1] 0.06250576 0.50534011

Pronto, agora ficou mais fácil. Constata-se que as correlações propostas, para esses tamanhos de amostra, mostram-se incapazes de rejeitar a H0 de que elas sejam iguais.Uma alternativa é construir duas matrizes com as correlações indicadas e usar um teste de Bartlett usando:

html image example

(R1 = matrix(c(1, r1, r1, 1), nrow = 2))

##      [,1] [,2]
## [1,]  1.0  0.7
## [2,]  0.7  1.0

(R2 = matrix(c(1, r2, r2, 1), nrow = 2))

##      [,1] [,2]
## [1,]  1.0  0.3
## [2,]  0.3  1.0

cortest.bartlett(R1-R2, n1*n2/(n1+n2))

## $chisq
## [1] 3.906503
## 
## $p.value
## [1] 0.04809959
## 
## $df
## [1] 1

Note-se que o teste de Bartlett é um pouco mais sensível que a construção de intervalos de confiança, o que implica — neste caso — na rejeição da H0, também por pequena margem. De qualquer forma, os valores estão tão próximos do valor crítico que a prudência indica a necessidade de refazer os testes com amostras maiores.

Capítulo 18 - Testes não-paramétricos

Os testes não-paramétricos são uma opção para quando os dados coletados violarem o pressuposto de normalidade, que é o requisito para que se possa usar um teste paramétrico. Ao contrário da estatística paramétrica, esses testes sem distribuição podem ser usados com dados quantitativos e qualitativos. Na estatística, os testes não paramétricos são métodos de análise estatística que não requerem uma distribuição para atender às premissas necessárias a serem analisadas (especialmente se os dados não são normalmente distribuídos). Devido a essa razão, eles são às vezes chamados de testes sem distribuição.

Note que testes não paramétricos são usados como um método alternativo para testes paramétricos, não como seus substitutos. Em outras palavras, se os dados atenderem às premissas necessárias para a realização dos testes paramétricos, o teste paramétrico relevante deve ser aplicado.Além disso, em alguns casos, mesmo que os dados não atendam às suposições necessárias, mas o tamanho da amostra dos dados seja grande o suficiente, ainda podemos aplicar os testes paramétricos em vez dos testes não paramétricos.

Teste de Sinais

O teste de sinais é usado para a comparação de uma amostra com um valor predeterminado, ou de duas amostras pareadas.

Exercício 18.1 - Página 451

Para conferir a alegação de um professor de que o valor publicado de 9.050 para coeficiente de fricção de metais bem engraxados deve ser muito pequeno, numa turma de ciências fazem 18 determinações do coeficiente (conforme dados abaixo).Normalmente o teste $t$ de uma amostra seria a escolha lógica para testar a alegação, mas a assimetria dos dados suere o uso de uma alternativa não paramétrica. Portanto o professor sugere que a turma use o teste de sinais para uma amostra para atestar a hipótese nula ū>0.05 ao nivel 0.05 de significância.

a = c(54, 52, 44, 56, 50, 51, 55, 53, 47, 53, 52, 50, 51, 51, 54, 46, 53, 43)/1000
a

##  [1] 0.054 0.052 0.044 0.056 0.050 0.051 0.055 0.053 0.047 0.053 0.052
## [12] 0.050 0.051 0.051 0.054 0.046 0.053 0.043

thresh = 0.050
w = which(a==thresh)
w

## [1]  5 12

a = a[-w]
a

##  [1] 0.054 0.052 0.044 0.056 0.051 0.055 0.053 0.047 0.053 0.052 0.051
## [12] 0.051 0.054 0.046 0.053 0.043

n = length(a)
n

## [1] 16

a>thresh

##  [1]  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE
## [12]  TRUE  TRUE FALSE  TRUE FALSE

paste(c("-", "+")[(a>thresh)+1], collapse = " ")

## [1] "+ + - + + + + - + + + + + - + -"

s = sum(a>thresh)
s

## [1] 12

binom.test(s, n, p = 0.5, alternative = "greater")

## 
##  Exact binomial test
## 
## data:  s and n
## number of successes = 12, number of trials = 16, p-value = 0.03841
## alternative hypothesis: true probability of success is greater than 0.5
## 95 percent confidence interval:
##  0.5156036 1.0000000
## sample estimates:
## probability of success 
##                   0.75

Este teste pode ser realizado diretamente pela função $SignTest$ por meio da biblioteca $DescTools$.

library(DescTools)
SignTest(a, mu = thresh, alternative = "greater")

## 
##  One-sample Sign-Test
## 
## data:  a
## S = 12, number of differences = 16, p-value = 0.03841
## alternative hypothesis: true median is greater than 0.05
## 96.2 percent confidence interval:
##  0.051   Inf
## sample estimates:
## median of the differences 
##                     0.052

**Ou então através da biblioteca $BSDA$ com a função $SIGN.test$.

library(BSDA)

## Loading required package: lattice

## 
## Attaching package: 'BSDA'

## The following object is masked from 'package:datasets':
## 
##     Orange

library(lattice)
SIGN.test(a, md = thresh, alternative = "greater")

## 
##  One-sample Sign-Test
## 
## data:  a
## s = 12, p-value = 0.03841
## alternative hypothesis: true median is greater than 0.05
## 95 percent confidence interval:
##  0.051   Inf
## sample estimates:
## median of x 
##       0.052 
## 
## Achieved and Interpolated Confidence Intervals: 
## 
##                   Conf.Level L.E.pt U.E.pt
## Lower Achieved CI     0.8949  0.051    Inf
## Interpolated CI       0.9500  0.051    Inf
## Upper Achieved CI     0.9616  0.051    Inf

Em ambos os casos, a 95% de confiança, rejeita-se a hipótese nula de que o professor esteja correto quanto ao índice de atrito.

Comparando com um teste t, verifica-se que o teste de sinais apresentou maior potência que o teste paramétrico de mesma finalidade, ou seja:

t.test(a, mu = 0.050, alternative = "greater")

## 
##  One Sample t-test
## 
## data:  a
## t = 0.9641, df = 15, p-value = 0.1751
## alternative hypothesis: true mean is greater than 0.05
## 95 percent confidence interval:
##  0.04923282        Inf
## sample estimates:
## mean of x 
## 0.0509375

Exercício 18.2 - Página 452

No exemplo 12.9 (pg.318) foram fornecidos os dados relativos às perdas semanais médias devidas a acidentes em dez indústrias, antes e depois da adoção de um programa de segurança abrangente.Use o teste de sinais com pares de dados para refazer este exercício.

acid = c(45, 36, 73, 60, 46, 44, 124, 119, 33, 35, 57, 51, 83, 77, 34, 29, 26, 24, 17, 11)
acid = matrix(acid, nrow = 2)
acid

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]   45   73   46  124   33   57   83   34   26    17
## [2,]   36   60   44  119   35   51   77   29   24    11

SignTest(x = acid[1,], y = acid[2,], paired = TRUE, alternative = "greater")

## 
##  Dependent-samples Sign-Test
## 
## data:  acid[1, ] and acid[2, ]
## S = 9, number of differences = 10, p-value = 0.01074
## alternative hypothesis: true median difference is greater than 0
## 98.9 percent confidence interval:
##    2 Inf
## sample estimates:
## median of the differences 
##                       5.5

Como 0.01074 é menor do que 0.05, a hipótese nula deve ser rejeitada (hipótese nula geralmente afirma que não existe relação entre dois fenômenos medidos) - e como no exemplo 12.9 pode-se concluir que o programa de seguraça é eficaz.

Teste de Sinais com Posto ou Teste Wicoxon

O teste de Wilcoxon, como o teste de sinais, compara uma amostra com um parâmetro estabelecido ou testa duas amostras pareadas.

wilcox.test(a, mu = thresh, alternative = "greater", exact = FALSE)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  a
## V = 84, p-value = 0.211
## alternative hypothesis: true location is greater than 0.05

O teste Wilcoxon de amostras pareadas (também conhecido como teste de postos sinalizados de Wilcoxon) é uma alternativa não paramétrica ao teste t pareado usado para comparar dados pareados. É usado quando seus dados não são normalmente distribuídos.

Exercício 18.5 - Página 458

Feitas as medições de octanagem de uma certa marca de gasolina, testar a hipótese nula de que a octanagem seja maior ou igual a 98.5.

gas = c(975, 952, 973, 960, 968, 1003, 974, 953, 932, 991, 961, 976, 982, 985, 949)/10
wilcox.test(x = gas, mu = 98.5, alternative = "less", exact = FALSE)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  gas
## V = 10, p-value = 0.004187
## alternative hypothesis: true location is less than 98.5

Como V=10 e menor do que 16 a hipótese nula deve ser rejeitada - portanto, conclui-se que a octanagem é menor do que 98.5

Exercício 18.6 - Página 459

Use o teste de Wilcoxon para refazer o exemplo 18.2, referente às perdas semanais médias de horas de trabalho.

a = c(45, 73, 46, 124, 33, 57, 83, 34, 26, 17) 
b = c(36, 60, 44, 119, 35, 51, 77, 29, 24, 11) 

wilcox.test(a, b, paired = TRUE, alternative = "greater", exact = FALSE, correct = TRUE)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  a and b
## V = 53, p-value = 0.005185
## alternative hypothesis: true location shift is greater than 0

O resultado nos leva a entender que a hipótes nula deve ser rejeitada, e, portanto conclui-se que o programa de segurança é eficaz.

Grandes amostras

Usando o R, não há necessidade de adotar um procedimento alternativo para quando houver grandes amostras. Basta aplicar o teste diretamente.

set.seed(11) # para replicabilidade
x = rnorm(200, 10, 2)
y = rnorm(250, 10.4, 1)
wilcox.test(x, y)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  x and y
## W = 20988, p-value = 0.00343
## alternative hypothesis: true location shift is not equal to 0

Exercício 18.8 - Página 465

O teste U (e Mann-Whitney) é adotado para a comparação de duas amostras independentes. Na verdade, o nome completo é teste de Mann-Whitney-Wilcoxon e é uma extensão do teste de Wilcoxon a amostras independentes.Esse teste é implementado no R na mesma função que o teste de Wilcoxon.

Exercício 18.8Suponha que se queira comparar o tamanho do grão de areia obtido em duas localidades diferentes da superfície da lua, com base nos seguintes diâmetros (em milímetros):

Por meio do teste U ao nível 0.05 de significância para testar se as amostras provêm, ou não de populações com médias iguais.

l1 = c(37, 70, 75, 30, 45, 16, 62, 73, 33)/100
l2 = c(86, 55, 80, 42, 97, 84, 24, 51, 92, 69)/100

wilcox.test(l1, l2, alternative = "two.sided")

## 
##  Wilcoxon rank sum test
## 
## data:  l1 and l2
## W = 24, p-value = 0.09472
## alternative hypothesis: true location shift is not equal to 0

Não há como rejeitar a hipótese nula de que os diâmetros dos grãos de areia sejam iguais.

Exercício 18.11 - Página 468

Distribuem-se os estudantes aleatoriamente por grupos que estudam espanhol por três métodos diferentes: (1) ensino em sala de aula e laboratório de linguagem, (2) apenas ensino em sala de aula e (3) apenas auto-estudo em laboratório de linguagem. São apresentadas, a seguir, as notas do exame final de amostras de estudantes dos três grupos:

Use o teste H ao nível 0.05 de significância para testar a hipótese nula de que as populações originais são idênticas, contra a hipótese alternativa de que suas médias não são todas iguais.

m1 = c(94, 88, 91, 74, 86, 97)
m2 = c(85, 82, 79, 84, 61, 72, 80)
m3 = c(89, 67, 72, 76, 69)
kruskal.test(list(m1, m2, m3))

## 
##  Kruskal-Wallis rank sum test
## 
## data:  list(m1, m2, m3)
## Kruskal-Wallis chi-squared = 6.6731, df = 2, p-value = 0.03556

O teste rejeita a H0 de que as médias das populações sejam todas iguais. Para saber qual média (ou médias) é diferente das demais, usa-se pares de testes de Mann-Whitney (implementados pela função wilcox.test):

wilcox.test(m1, m2)

## 
##  Wilcoxon rank sum test
## 
## data:  m1 and m2
## W = 37, p-value = 0.02214
## alternative hypothesis: true location shift is not equal to 0

wilcox.test(m1, m3)

## 
##  Wilcoxon rank sum test
## 
## data:  m1 and m3
## W = 26, p-value = 0.05195
## alternative hypothesis: true location shift is not equal to 0

#wilcox.test(m2, m3)

Portanto o que posso afirmar é que a amostra A se destaca das outras

Exercício 18.12 - Página 472

Considerando uma alameda de árvores plantadas ao longo de uma estrada, verificou-se um padrão de árvores saudáveis (H) e doentes (D): HHHH, DDD, HHHHHHH, DD, HH, DDDD.O problema consiste em testar se esse padrão pode ser considerado aleatório ou não,com significância α = 0.05.

library(randtests) 
a = c(1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 2, 2, 2, 2) 
runs.test(a, threshold = 1.5)

## 
##  Runs Test
## 
## data:  a
## statistic = -2.5513, runs = 6, n1 = 9, n2 = 13, n = 22, p-value =
## 0.01073
## alternative hypothesis: nonrandomness

Exercício 18.15 - Página 477

O número de horas durante as quais 10 alunos estudaram para um exame, e as notas obtidas foram:

a = c(9, 5, 11, 13, 10, 5, 18, 15, 2, 8) 
b = c(56, 44, 79, 72, 70, 54, 94, 85, 33, 65) 
#ct = cor.test(a, b, method = "spearman") 
#ct$p.value
#ct$estimate

Exercícios da Lista (Julho/2017)

Exercício 01

Um teste psicológico compreende três habilidades: 1) compreensão de leitura; 2) habilidade aritmética; e 3) domínio de vocabulário. As correlações medidas entre as habilidades são $r_{12}$ = 0.05, $r_{13}$ = 0.73 e $r_{23}$ = 0.20. Qual a correlação entre a compreensão de leitura e a habilidade aritmética, controlando o domínio de vocabulário?

Sabendo que:

$$ r_{1,2,3}=

tem-se que:

r12 = 0.05 
r13 = 0.73 
r23 = 0.20 
r_123 = (r12-r13*r23)/(sqrt(1-r13^2)*sqrt(1-r23^2)) 
r_123

## [1] -0.1433609

library(MASS)
Sigma = matrix(c( 
  1, r12, r13, 
  r12, 1, r23, 
  r13, r23, 1 ), nrow = 3) 
x = mvrnorm(100, mu = c(0, 0, 0), Sigma = Sigma, empirical = TRUE) 

cor(x)

##      [,1] [,2] [,3]
## [1,] 1.00 0.05 0.73
## [2,] 0.05 1.00 0.20
## [3,] 0.73 0.20 1.00

cor(x)[1,2]

## [1] 0.05

cor(x)[1, 2]

## [1] 0.05

Exercício 02

Determine a correlação, a sua signiﬁcância e o seu intervalo de conﬁança de 95%para os conjuntos de dados:

x1 = c(3, 3, 3, 4, 5, 6, 7, 8, 9, 10, 7, 6, 8, 4) 
x2 = c(5, 5, 3, 5, 2, 7, 6, 5, 4, 11, 4, 2, 3, 2) 
ct = cor.test(x1, x2) 
ct$estimate

##      cor 
## 0.431634

ct$p.value

## [1] 0.1232926

ct$conf.int

## [1] -0.1283366  0.7829137
## attr(,"conf.level")
## [1] 0.95

Exercício 03

O que se pode aﬁrmar sobre a sua igualdade de médias? É um caso para uso de estatística paramétrica ou nãoparamétrica? Discuta o resultado.

a = c(9, 5, 11, 13, 10, 5, 18, 15, 2, 8) 
b = c(56, 44, 79, 72, 70, 54, 94, 85, 33, 65) 
set.seed(123) 
a = rexp(120, 1) 
b = rexp(330, 2) 
b = max(b)-b 
a = a+mean(b)-mean(a)

shapiro.test(a)

## 
##  Shapiro-Wilk normality test
## 
## data:  a
## W = 0.75793, p-value = 0.0000000000008738

shapiro.test(b)

## 
##  Shapiro-Wilk normality test
## 
## data:  b
## W = 0.86722, p-value = 0.0000000000000003202

median(b)

## [1] 1.83681

x = c(a, b) 
g = as.factor(c(rep(1, length(a)), rep(2, length(b)))) 
g

##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [71] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [106] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [141] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [176] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [211] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [246] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [281] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [316] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [351] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [386] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [421] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## Levels: 1 2

summary(aov(x~g))

##              Df Sum Sq Mean Sq F value Pr(>F)
## g             1    0.0  0.0000       0      1
## Residuals   448  184.4  0.4116

vt = var.test(a, b) 
t.test(a, b, var.equal = vt$p.value > 0.05)

## 
##  Welch Two Sample t-test
## 
## data:  a and b
## t = 0, df = 138.48, p-value = 1
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1841191  0.1841191
## sample estimates:
## mean of x mean of y 
##  1.682405  1.682405

wilcox.test(a, b)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  a and b
## W = 16857, p-value = 0.01587
## alternative hypothesis: true location shift is not equal to 0

kruskal.test(list(a, b))

## 
##  Kruskal-Wallis rank sum test
## 
## data:  list(a, b)
## Kruskal-Wallis chi-squared = 5.8196, df = 1, p-value = 0.01585