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.
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
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.
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:
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.
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.
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.
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.
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.
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.
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.
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.
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:
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.
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!
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.
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.
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.
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.