Façamos painel de dados. No painel, um valor se repete diversas vezes na coluna = dados longitudinais <> dados transversais (ou cross-section).

regressão logística é diferente da múltipla pois depende da variável que se quer explicar.

1. Regressão em painel

Quando se usa um painel de dados, é possível controlar fontes não observadas ou imensuráveis de heterogeneidade individual que varia para cada indivíduo mas não com o tempo.

Faremos regressão por dados em painel empilhado (ou pooled), quando o indivíduo é observado ao longo do tempo.

p<-read.csv2("Dados_Painel.csv")
summary(p)
##     empresa        ano            eva                 roa         
##  Ambev  : 4   Min.   :2010   Min.   :-58811845   Min.   : -3.700  
##  Brf    : 4   1st Qu.:2011   1st Qu.: -1968079   1st Qu.:  6.175  
##  Cbd    : 4   Median :2012   Median :  -159436   Median : 13.450  
##  Ccr    : 4   Mean   :2012   Mean   : -3942511   Mean   : 26.850  
##  Cielo  : 4   3rd Qu.:2012   3rd Qu.:   541395   3rd Qu.: 29.200  
##  Jbs    : 4   Max.   :2013   Max.   : 15251592   Max.   :178.200  
##  (Other):16                                                       
##       roe               ml              gaf             imobpl     
##  Min.   :-2.700   Min.   :-4.900   Min.   :-0.500   Min.   : 12.5  
##  1st Qu.: 2.850   1st Qu.: 1.975   1st Qu.: 1.150   1st Qu.: 28.3  
##  Median : 5.400   Median : 7.450   Median : 1.430   Median : 71.2  
##  Mean   : 8.805   Mean   :14.385   Mean   : 1.423   Mean   : 62.4  
##  3rd Qu.:13.075   3rd Qu.:26.750   3rd Qu.: 1.725   3rd Qu.: 79.0  
##  Max.   :48.500   Max.   :45.900   Max.   : 2.710   Max.   :152.8  
##                                                                    
##     ccapter      
##  Min.   :0.0000  
##  1st Qu.:0.0875  
##  Median :0.1150  
##  Mean   :0.1610  
##  3rd Qu.:0.1700  
##  Max.   :0.7500  
## 
dim(p)
## [1] 40  9
boxplot(p$eva ~ p$ano,col="orange",xlab="ano",ylab="eva",
        ylim=c(-6000000,9000000))

Objetivo: explicar valor econômico (economical value ou eva) a partir das outras variáveis. Note que o eva mediano é negativo!

Comecemos criando um objeto que receberá a equação.

eq=eva ~ roa + roe + ml + gaf + imobpl + ccapter
eq
## eva ~ roa + roe + ml + gaf + imobpl + ccapter

1.1 Regressão múltipla

Faremos a seguir uma regressão múltipla, como vimos na aula passada.

mult=lm(eq,data=p)
summary(mult)
## 
## Call:
## lm(formula = eq, data = p)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -29566622  -5228631   1526712   6010800  24830720 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 17023516    9005231   1.890 0.067513 .  
## roa          -185924     144154  -1.290 0.206098    
## roe           736644     772179   0.954 0.347029    
## ml            114613     308623   0.371 0.712738    
## gaf         -6366668    3510100  -1.814 0.078809 .  
## imobpl       -251698      66875  -3.764 0.000654 ***
## ccapter      4074153   14852336   0.274 0.785556    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11590000 on 33 degrees of freedom
## Multiple R-squared:  0.4549, Adjusted R-squared:  0.3558 
## F-statistic: 4.589 on 6 and 33 DF,  p-value: 0.001723
par(mar=c(2,2,2,2))

As variáveis que mais impactam são imobpl e gaf.

1;2 Regressão com dados empilhados

Primeiramente, instale e carregue o pacote pml.

library("plm")
## Loading required package: Formula

Existem três formas básicas de trabalhar com regressão em painel:

  • pooled (pooling)
  • efeitos fixos (thin)
  • efeitos aleatórios (random)

Criemos primeiro um data frame do tipo painel de dados a partir do arquivo lido.

# cria painel de dados
pdf<-pdata.frame(p)
pooled<-plm(eq,data=pdf,model="pooling")
summary(pooled)
## Pooling Model
## 
## Call:
## plm(formula = eq, data = pdf, model = "pooling")
## 
## Balanced Panel: n = 10, T = 4, N = 40
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -29566622  -5228631   1526712   6010800  24830720 
## 
## Coefficients:
##             Estimate Std. Error t-value  Pr(>|t|)    
## (Intercept) 17023516    9005231  1.8904 0.0675129 .  
## roa          -185924     144154 -1.2898 0.2060976    
## roe           736644     772179  0.9540 0.3470290    
## ml            114613     308623  0.3714 0.7127376    
## gaf         -6366668    3510100 -1.8138 0.0788088 .  
## imobpl       -251698      66875 -3.7637 0.0006545 ***
## ccapter      4074153   14852336  0.2743 0.7855559    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    8.1328e+15
## Residual Sum of Squares: 4.4335e+15
## R-Squared:      0.45487
## Adj. R-Squared: 0.35575
## F-statistic: 4.58928 on 6 and 33 DF, p-value: 0.0017233

O resultado gerado pelo método pooling é rigorosamente o mesmo daquele obtido pela regressão múltipla acima.

1.3 Regressão múltipla com efeito idiossincrático

Considere os dados relativos à Petrobras.

head(subset(p,p$empresa=="Petrobras"))
##     empresa  ano       eva  roa roe   ml  gaf imobpl ccapter
## 1 Petrobras 2010 -26944548 15.2 6.9 16.8 1.77   91.2    0.06
## 2 Petrobras 2011 -29751999 10.3 5.5 13.6 1.63  103.0    0.02
## 3 Petrobras 2012 -51029924  6.2 3.1  7.4 1.59  121.2    0.04
## 4 Petrobras 2013 -58811845  6.6 3.1  7.5 1.78  152.8    0.05

Questão
Note que o eva da Petrobras muda com o passar dos anos. Será que existe algum elemento que altera o eva de uma empresa ao longo dos anos? Esse é o assim chamado efeito não observado ou efeito idiossincrático.

Iremos agora incluir esse efeito a nossa análise, uma constante para cada empresa. O nome da empresa varia de acordo com ela. Esse valor será incluído na equação.

Vejamos quantas empresas fazem parte do banco de dados.

table(p$empresa)
## 
##     Ambev       Brf       Cbd       Ccr     Cielo       Jbs    Kroton 
##         4         4         4         4         4         4         4 
## Petrobras  Ultrapar      Vale 
##         4         4         4

Dez empresas. Em seguida faremos nova regressão múltipla a partir dos dados originais, mas desta vez incluindo a “constante” empresa.

eq2=eva ~ roa + roe + ml + gaf + imobpl + ccapter + empresa
mult_fixa=lm(eq2,data=p)
summary(mult_fixa)
## 
## Call:
## lm(formula = eq2, data = p)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -5126686 -1166281   -37188   894198  7556015 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -9211547    6185097  -1.489  0.14943    
## roa                  -7305     174208  -0.042  0.96690    
## roe                 387193     514947   0.752  0.45942    
## ml                  613992     180364   3.404  0.00233 ** 
## gaf               -2056804    1977846  -1.040  0.30874    
## imobpl             -345850      70388  -4.913 5.18e-05 ***
## ccapter           23515871    7321729   3.212  0.00373 ** 
## empresaBrf        28800685    4774458   6.032 3.15e-06 ***
## empresaCbd        29729832    4735833   6.278 1.73e-06 ***
## empresaCcr         1274578    5989102   0.213  0.83327    
## empresaCielo     -18553296   12155287  -1.526  0.13999    
## empresaJbs        30487213    4876130   6.252 1.84e-06 ***
## empresaKroton       811342    4051238   0.200  0.84296    
## empresaPetrobras   1856497    6476095   0.287  0.77683    
## empresaUltrapar   33606307    5188643   6.477 1.07e-06 ***
## empresaVale       31445378    6074130   5.177 2.66e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3382000 on 24 degrees of freedom
## Multiple R-squared:  0.9662, Adjusted R-squared:  0.9451 
## F-statistic: 45.79 on 15 and 24 DF,  p-value: 5.522e-14

O resultado mudou. A variável imobpl continua significante, mas agora a margem líquida (ml) passou a ser significante.

Além disso, um coeficiente foi criado para cada uma das dez empresas. E R^2 ficou inflacionado, quase igual a 100%, após a adição da nova variável dummy.

1.4 Regressão de dados em painel com efeitos fixos

Tentaremos agora remover essas constantes a fim de não alterar o valor de R^2.

efeitos_fixos=plm(eq,data=pdf,model="within")
summary(efeitos_fixos)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = eq, data = pdf, model = "within")
## 
## Balanced Panel: n = 10, T = 4, N = 40
## 
## Residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -5126686 -1166281   -37188        0   894198  7556015 
## 
## Coefficients:
##           Estimate Std. Error t-value  Pr(>|t|)    
## roa        -7305.2   174207.7 -0.0419  0.966899    
## roe       387192.9   514947.0  0.7519  0.459417    
## ml        613992.1   180363.9  3.4042  0.002333 ** 
## gaf     -2056804.1  1977846.4 -1.0399  0.308740    
## imobpl   -345850.3    70388.3 -4.9135 5.177e-05 ***
## ccapter 23515871.4  7321728.9  3.2118  0.003733 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1.5653e+15
## Residual Sum of Squares: 2.7457e+14
## R-Squared:      0.82459
## Adj. R-Squared: 0.71496
## F-statistic: 18.8037 on 6 and 24 DF, p-value: 5.4245e-08

O resultado mostra que as empresas com melhor margem líquida (ml) ou se endivida mais a partir de capital de terceiros (ccapter) têm melhor desempenho econômico.

Note como agora o valor de R^2 está menor. O efeito fixo consegue, portanto, não influenciar tão fortemente R^2.

Questão
Ainda é possível obter os valores para cada empresa? Sim, usando a função fixef.

fixef(efeitos_fixos)
##     Ambev       Brf       Cbd       Ccr     Cielo       Jbs    Kroton 
##  -9211547  19589138  20518285  -7936970 -27764843  21275666  -8400205 
## Petrobras  Ultrapar      Vale 
##  -7355050  24394759  22233830

E agora corrijamos a variável dummy do primeiro valor, ie, comparando todas as empresas com a primeira.

# note que o primeiro elemento corresponde á Ambev
fixef(efeitos_fixos) - fixef(efeitos_fixos)[1]
##       Ambev         Brf         Cbd         Ccr       Cielo         Jbs 
##         0.0  28800685.0  29729832.1   1274577.6 -18553295.9  30487213.2 
##      Kroton   Petrobras    Ultrapar        Vale 
##    811342.4   1856497.2  33606306.5  31445377.5

Veja que a Ambev tem valor zero, claro, se comparada a si mesma. Estes valores agora correspondem exatamente àqueles obtidos com a regressão fixa, com efeito idiossincrático.

1.5 Regressão de dados em painel com efeitos aleatórios

Se o efeito se manifesta de modo aleatório em relação ao tempo e ao indivíduo, não se pode aplicar a regressão com efeitos fixos já que esta parte do pressuposto de que os efeitos dos indivíduos não têm correlação.

efeitos_aleatorios=plm(eq,data=pdf,model="random")
summary(efeitos_aleatorios)
## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = eq, data = pdf, model = "random")
## 
## Balanced Panel: n = 10, T = 4, N = 40
## 
## Effects:
##                     var   std.dev share
## idiosyncratic 1.144e+13 3.382e+06 0.046
## individual    2.381e+14 1.543e+07 0.954
## theta: 0.891
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -7254557 -1209593   483685  2063074  6110377 
## 
## Coefficients:
##             Estimate Std. Error t-value  Pr(>|t|)    
## (Intercept)  4614026    7509996  0.6144  0.543174    
## roa          -167573     149863 -1.1182  0.271568    
## roe           714512     475351  1.5031  0.142317    
## ml            573257     176840  3.2417  0.002716 ** 
## gaf         -2127947    1974343 -1.0778  0.288939    
## imobpl       -304904      65282 -4.6705 4.854e-05 ***
## ccapter     21487748    7342876  2.9263  0.006165 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1.6433e+15
## Residual Sum of Squares: 3.956e+14
## R-Squared:      0.75926
## Adj. R-Squared: 0.71549
## F-statistic: 17.3464 on 6 and 33 DF, p-value: 6.0419e-09

No resultado mostrado acima, pode-se ver que tanto ml quando imobpl e ccapter continuam significantes. R^2 tampouco foi bastante influenciado.

1.6 Avaliação da Regressão de dados em painel

Vistos esses três modelos, fica a dúvida de qual usar. Faremos testes para comparar os modelos dois a dois, cada um com grau de complexidade, sendo o de efeitos aleatórios o mais complexo.

  • Para a decisão entre os modelos pooled e de efeitos fixos faremos um teste F para os efeitos dos indivíduos.
  • Para decidir entre os modelos pooled e de efeitos aleatórios usaremos o Teste Multiplicador de Lagrange de Breusch-Pagan.
  • E para escolher entre o efeitos fixos e o aleatório, faremos o teste de Hausman.
1.6.1 Teste entre pooled e efeitos fixos
  • H_0 : efeitos dos indivíduos não são significantes (pooled é o melhor modelo)
  • H_a : efeitos dos indivíduos são significantes (efeitos fixos é o melhor)
pFtest(efeitos_fixos,pooled)
## 
##  F test for individual effects
## 
## data:  eq
## F = 40.391, df1 = 9, df2 = 24, p-value = 2.414e-12
## alternative hypothesis: significant effects

Como esperado, o valor p é praticamente zero. Logo, a hipótese alternativa ganha: o método dos efeitos fixos é neste caso mais recomendável que uma regressão múltipla convencional.

1.6.2 Teste entre pooled e efeitos aleatórios
  • H_0 : efeito do painel não é significante (pooled é o melhor modelo)
  • H_a : efeito do painel é significante (efeitos aleatórios é o melhor)
plmtest(pooled,type="bp")
## 
##  Lagrange Multiplier Test - (Breusch-Pagan) for balanced panels
## 
## data:  eq
## chisq = 28.314, df = 1, p-value = 1.031e-07
## alternative hypothesis: significant effects

Novamente, o valor p é minúsculo. Logo, a hipótese alternativa ganha: o método dos efeitos aleatórios é neste caso mais indicado que uma regressão múltipla convencional.

1.6.3 Teste entre efeitos fixos e aleatórios
  • H_0 : a diferença entre os coeficientes não é sistemática (efeitos aleatórios é melhor)
  • H_a : a diferença entre os coeficientes é sistemática (efeitos fixos é melhor)
phtest(efeitos_fixos,efeitos_aleatorios)
## 
##  Hausman Test
## 
## data:  eq
## chisq = 7.146, df = 6, p-value = 0.3076
## alternative hypothesis: one model is inconsistent

Aqui o valor p = 0.30. Logo, para alpha=0.05, não se tem motivos para rejeitar a H_0. O modelo de efeitos aleatórios, mais complexo, é mais indicado neste caso.

1.7 Diagnóstico do modelo

Faremos agora a validação do modelo mais adequado ao banco de dados, ie, o modelo de efeitos aleatórios. A validação será feita de acordo com:

  • R^2
  • Teste F
  • normalidade dos resíduos
  • Homocedasticidade dos resíduos
  • independência dos resíduos
  • colinearidade

Repetindo, este é o resumo do modelo:

summary(efeitos_aleatorios)
## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = eq, data = pdf, model = "random")
## 
## Balanced Panel: n = 10, T = 4, N = 40
## 
## Effects:
##                     var   std.dev share
## idiosyncratic 1.144e+13 3.382e+06 0.046
## individual    2.381e+14 1.543e+07 0.954
## theta: 0.891
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -7254557 -1209593   483685  2063074  6110377 
## 
## Coefficients:
##             Estimate Std. Error t-value  Pr(>|t|)    
## (Intercept)  4614026    7509996  0.6144  0.543174    
## roa          -167573     149863 -1.1182  0.271568    
## roe           714512     475351  1.5031  0.142317    
## ml            573257     176840  3.2417  0.002716 ** 
## gaf         -2127947    1974343 -1.0778  0.288939    
## imobpl       -304904      65282 -4.6705 4.854e-05 ***
## ccapter     21487748    7342876  2.9263  0.006165 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1.6433e+15
## Residual Sum of Squares: 3.956e+14
## R-Squared:      0.75926
## Adj. R-Squared: 0.71549
## F-statistic: 17.3464 on 6 and 33 DF, p-value: 6.0419e-09

Nota-se que R^2 = 73,95%, logo o modelo consegue explicar bem os dados. O teste F retornou valor p = 6.0e-09. Logo, o modelo é válido.

1.7.1 normalidade dos resíduos
residuos<-efeitos_aleatorios$residuals
hist(residuos,col="orange",main="Histograma dos resíduos")

shapiro.test(residuos)$p.value
## [1] 0.2018046

Aqui, o valor p=0.2018. Passou!

1.7.2 Homocedasticidade dos resíduos
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(efeitos_aleatorios)$p.value
##           BP 
## 0.0002138632

Como o valor p, 0.0002139, é menor que alpha=0.05, o modelo tem problemas de heterocedasticidade (os resíduos estão muito altos para algumas observações). Nesse caso, pode-se tentar transformar a variável ou aplicar correção ao teste, chamada de correção de White.

1.7.3 Independência dos resíduos

Para uma análise de dados em painel, usamos o teste de Breusch-Godfrey/Wooldridge já que cada indivíduo é observado mais de uma vez.

pbgtest(efeitos_aleatorios)$p.value
## [1] 0.7453253

Aqui o valor p é alto, então o modelo passa neste teste.

1.7.4 Multicolinearidade das variáveis independentes

Procura-se aqui vif´s com valores < 5.

library(car)
## Loading required package: carData
vif(efeitos_aleatorios)
##       roa       roe        ml       gaf    imobpl   ccapter 
## 10.375198 16.818211  4.642237  1.928390  1.283805  2.294353

Os valores de retorno sobre o ativo (roa) e retorno sobre o retorno sobre o patrimônio líquido (roe ou return of equity) parecem estar superpostos. Logo, há multicolinearidade no modelo. Sugere-se reconstruir o modelo, removendo-se uma dessas variáveis.

1.7.5 Retirando o ROE

Fiquemos com o roa e refaçamos o modelo. Em seguida, verifiquemos novamente a colinearidade.

eq3=eva ~ roa + + ml + gaf + imobpl + ccapter 
efeitos_aleatorios2=plm(eq3,data=pdf,model="random")
vif(efeitos_aleatorios2)
##      roa       ml      gaf   imobpl  ccapter 
## 2.129892 3.387912 1.788065 1.272439 1.960844

Todos os valores abaixo de 5. Logo, o novo modelo passa no teste de colinearidade. Refaçamos os outros três testes.

residuos<-efeitos_aleatorios2$residuals
shapiro.test(residuos)$p.value
## [1] 0.3254487
bptest(efeitos_aleatorios2)$p.value
##           BP 
## 9.476019e-05
pbgtest(efeitos_aleatorios2)$p.value
## [1] 0.4015832

O novo modelo passa nos testes de normalidade e de independência dos resíduos, mas continua com problemas de heterocedasticidade.