Estatística Descritiva
Adhmir
02/04/2022
O intuito é abordar as principais medidas de tendência e dispersão, gráficos e modelos de regressão com dados em painel
1 Estatística descritiva
Primeiro vamos configurar a pasta de trabalho
setwd("C:/Users/adh_r/OneDrive - FURB/PHD/A parte/Ideia oficinas/Software R/markdowns")Importar as bibliotecas
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.3 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(readxl)Vamos ler o último arquivo criado
df <- read_excel("hhi_winsor.xlsx")
head(df,3)## # A tibble: 3 × 10
## CompanyName TRBCSector Ano AtivoTotal EBITDA tamanho ROA marketShare
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Kepler Weber SA Industrials 2019 949. 113. 6.86 0.119 0.00188
## 2 Kepler Weber SA Industrials 2018 676. 96.8 6.52 0.143 0.00149
## 3 Kepler Weber SA Industrials 2017 704. 52.6 6.56 0.0746 0.00178
## # ℹ 2 more variables: ihh <dbl>, concentracao <dbl>Podemos observar que o ano está em ordem decrescente, vamos organizar em ordem crescente pelo ano e nome da empresa.
df <- df %>% arrange(Ano, CompanyName)
head(df)## # A tibble: 6 × 10
## CompanyName TRBCSector Ano AtivoTotal EBITDA tamanho ROA marketShare
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3R Petroleum O… Energy 2014 846. 22.6 6.74 0.0267 0.000826
## 2 ATMA Participa… Industria… 2014 2633. 25.0 7.88 0.00949 0.00848
## 3 Aeris Industri… Energy 2014 1233. 166. 7.12 0.134 0.00120
## 4 Aes Brasil Ene… Utilities 2014 4637. 1390. 8.44 0.300 0.00800
## 5 Afluente Trans… Utilities 2014 89.4 11.2 4.49 0.126 0.000154
## 6 AgroGalaxy Par… Consumer … 2014 2236. 103. 7.71 0.0462 0.0142
## # ℹ 2 more variables: ihh <dbl>, concentracao <dbl>
Pode organizar de forma fácil quantas colunas quiser pelo comando
arrange()
Vamos observar a estrutura do conjunto de dados
str(df)## tibble [1,848 × 10] (S3: tbl_df/tbl/data.frame)
## $ CompanyName : chr [1:1848] "3R Petroleum Oleo e Gas SA" "ATMA Participacoes SA" "Aeris Industria e Comercio de Equipamentos para Geracao de Energia SA" "Aes Brasil Energia SA" ...
## $ TRBCSector : chr [1:1848] "Energy" "Industrials" "Energy" "Utilities" ...
## $ Ano : chr [1:1848] "2014" "2014" "2014" "2014" ...
## $ AtivoTotal : num [1:1848] 846.1 2633.2 1233.4 4636.8 89.4 ...
## $ EBITDA : num [1:1848] 22.6 25 165.6 1390.2 11.2 ...
## $ tamanho : num [1:1848] 6.74 7.88 7.12 8.44 4.49 ...
## $ ROA : num [1:1848] 0.02674 0.00949 0.1343 0.29982 0.12554 ...
## $ marketShare : num [1:1848] 0.000826 0.00848 0.001204 0.008 0.000154 ...
## $ ihh : num [1:1848] 7754 673 7754 926 926 ...
## $ concentracao: num [1:1848] 1 0 1 0 0 0 1 0 0 0 ...As três primeiras, inclusive o ano, colunas são caracteres e a demais numéricos.
1.1 Medidas de tendência e dispersão
Podemos criar um resumo estatístico do conjunto de dados inteiro
summary(df)## CompanyName TRBCSector Ano AtivoTotal
## Length:1848 Length:1848 Length:1848 Min. : 0.0
## Class :character Class :character Class :character 1st Qu.: 433.2
## Mode :character Mode :character Mode :character Median : 1849.2
## Mean : 13044.9
## 3rd Qu.: 7924.2
## Max. :987419.0
## EBITDA tamanho ROA marketShare
## Min. : -5504.00 Min. : 4.146 Min. :-19.33029 Min. :0.000000
## 1st Qu.: 21.17 1st Qu.: 6.071 1st Qu.: 0.04550 1st Qu.:0.001362
## Median : 164.16 Median : 7.523 Median : 0.10059 Median :0.005650
## Mean : 1525.55 Mean : 7.530 Mean : 0.06596 Mean :0.029221
## 3rd Qu.: 882.16 3rd Qu.: 8.978 3rd Qu.: 0.14990 3rd Qu.:0.022057
## Max. :145036.00 Max. :10.748 Max. : 1.13394 Max. :0.878841
## ihh concentracao
## Min. : 497.0 Min. :0.0000
## 1st Qu.: 672.4 1st Qu.:0.0000
## Median : 936.8 Median :0.0000
## Mean :1660.4 Mean :0.2798
## 3rd Qu.:2292.1 3rd Qu.:1.0000
## Max. :7753.5 Max. :1.0000Nesta saída temos mínimo, primeiro quartil, média, mediana, terceiro quartil e máximo.
Talvez, seja necessário ter a estatística por grupo, pode usar o
seguinte comando: by(df, df$Ano, summary)
No entanto, estamos acostumados com saídas mais elaboradas, que abordam
desde medidas de tendência até de dispersão. Ao invés de calcularmos,
podemos utilizar o pacote pastecs
#install.packages("pastecs")
library(pastecs)##
## Attaching package: 'pastecs'## The following objects are masked from 'package:dplyr':
##
## first, last## The following object is masked from 'package:tidyr':
##
## extractstat.desc(df)## CompanyName TRBCSector Ano AtivoTotal EBITDA tamanho
## nbr.val NA NA NA 1.848000e+03 1.848000e+03 1.848000e+03
## nbr.null NA NA NA 0.000000e+00 0.000000e+00 0.000000e+00
## nbr.na NA NA NA 0.000000e+00 0.000000e+00 0.000000e+00
## min NA NA NA 7.660000e-03 -5.503999e+03 4.146265e+00
## max NA NA NA 9.874190e+05 1.450360e+05 1.074779e+01
## range NA NA NA 9.874190e+05 1.505400e+05 6.601526e+00
## sum NA NA NA 2.410697e+07 2.819216e+06 1.391454e+04
## median NA NA NA 1.849221e+03 1.641615e+02 7.522519e+00
## mean NA NA NA 1.304490e+04 1.525550e+03 7.529511e+00
## SE.mean NA NA NA 1.333584e+03 1.792978e+02 4.362867e-02
## CI.mean NA NA NA 2.615491e+03 3.516477e+02 8.556670e-02
## var NA NA NA 3.286569e+09 5.940895e+07 3.517596e+00
## std.dev NA NA NA 5.732861e+04 7.707720e+03 1.875526e+00
## coef.var NA NA NA 4.394715e+00 5.052422e+00 2.490899e-01
## ROA marketShare ihh concentracao
## nbr.val 1848.00000000 1.848000e+03 1.848000e+03 1.848000e+03
## nbr.null 0.00000000 0.000000e+00 0.000000e+00 1.331000e+03
## nbr.na 0.00000000 0.000000e+00 0.000000e+00 0.000000e+00
## min -19.33028721 1.321555e-08 4.969510e+02 0.000000e+00
## max 1.13393645 8.788413e-01 7.753539e+03 1.000000e+00
## range 20.46422366 8.788413e-01 7.256588e+03 1.000000e+00
## sum 121.89008445 5.400000e+01 3.068358e+06 5.170000e+02
## median 0.10058617 5.649682e-03 9.368430e+02 0.000000e+00
## mean 0.06595784 2.922078e-02 1.660367e+03 2.797619e-01
## SE.mean 0.01316428 1.902335e-03 3.805870e+01 1.044477e-02
## CI.mean 0.02581843 3.730953e-03 7.464259e+01 2.048480e-02
## var 0.32025516 6.687688e-03 2.676762e+06 2.016043e-01
## std.dev 0.56591091 8.177828e-02 1.636081e+03 4.490036e-01
## coef.var 8.57988878 2.798635e+00 9.853735e-01 1.604949e+001.2 Gráficos
Apesar do R já tem um comando para gráficos plot(), podemos
usar o pacote ggplt2 para obter alguns mais elaborados
#install.packages("ggplot2")
library(ggplot2)
gf1 <- ggplot(df, aes(x=factor(TRBCSector), y=ihh, fill=factor(TRBCSector))) +
geom_boxplot(alpha=0.3) +
theme(legend.position="none",
axis.text.x=element_text(angle = -10, hjust = 0)) +
labs( x = "Setores", y = "concentracao")
gf1
Outra possibilidade é fazer um gráfico de densidade, pois a distribuição do dados fica bem visível
gf2 <- ggplot(df, aes(x = tamanho)) +
geom_density(colour = 1, fill = "blue")
gf2 +facet_grid(Ano ~ factor(TRBCSector))
Podemos criar gráficos de dispersão
gf3 <-ggplot(df, aes(x=marketShare,
y=ihh))+
geom_point()+
geom_smooth()+
#facet_grid(Ano ~ factor(TRBCSector))+
labs( x = "Participação de mercado", y = "concentracao")
gf3## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
Podemos analisar a correlação também. Mas, para isso precisamos criar um conjunto de dados com a correlação
library(ggcorrplot)
correlacao <- round(cor(df[4:9]), 2)
head(correlacao)## AtivoTotal EBITDA tamanho ROA marketShare ihh
## AtivoTotal 1.00 0.95 0.33 0.02 0.78 0.20
## EBITDA 0.95 1.00 0.29 0.03 0.72 0.18
## tamanho 0.33 0.29 1.00 0.12 0.45 0.00
## ROA 0.02 0.03 0.12 1.00 0.03 0.00
## marketShare 0.78 0.72 0.45 0.03 1.00 0.20
## ihh 0.20 0.18 0.00 0.00 0.20 1.00pvalor <- round(cor_pmat(df[4:9]),3)
head(pvalor)## AtivoTotal EBITDA tamanho ROA marketShare ihh
## AtivoTotal 0.000 0.000 0.000 0.378 0.000 0.000
## EBITDA 0.000 0.000 0.000 0.215 0.000 0.000
## tamanho 0.000 0.000 0.000 0.000 0.000 0.841
## ROA 0.378 0.215 0.000 0.000 0.199 0.864
## marketShare 0.000 0.000 0.000 0.199 0.000 0.000
## ihh 0.000 0.000 0.841 0.864 0.000 0.000E então criar um gráfico
#install.packages("ggcorrplot")
g4<- ggcorrplot(correlacao,
method = "circle",
hc.order = TRUE,
outline.col = "white",
type = "upper",
colors = c("darkblue", "white", "darkgreen"),
lab = TRUE,
p.mat = pvalor,
#insig = "blank"
)
g4
2 Análise de regressão dados em painel
Para termos como exemplo, vamos assumir o seguinte modelo:
\[ Ebitda= \alpha + \beta_1marketshare_{it} + \beta_2hhi_{it} + \beta_3tam_{it} + \varepsilon_{it}\]
Aqui vamos abordar três tópicos
- Dados empilhados (Pooled)
- Pressupostos
* Normalidade dos resíduos
* Vif
* heterocedasticidade
* autocorrelação
- Pressupostos
- Efeito fixo
- Teste F de chow
- Efeito Aleatório
- Lagrange Multiplier of Breusch Pagan
- Hausmann
2.1 Dados empilhados
Considerando o modelo proposto, vamos analisar o modelo de dados empilhado. Não será testado o pressuposto da normalidade das variáveis.
pooled <- lm(EBITDA~marketShare+ihh+tamanho ,data = df)
summary(pooled)##
## Call:
## lm(formula = EBITDA ~ marketShare + ihh + tamanho, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -50773 -133 267 502 85577
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.853e+02 5.782e+02 0.493 0.6218
## marketShare 6.878e+04 1.740e+03 39.530 <2e-16 ***
## ihh 1.712e-01 7.776e-02 2.202 0.0278 *
## tamanho -1.400e+02 7.433e+01 -1.883 0.0598 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5325 on 1844 degrees of freedom
## Multiple R-squared: 0.5235, Adjusted R-squared: 0.5227
## F-statistic: 675.2 on 3 and 1844 DF, p-value: < 2.2e-16Podemos observar o seguinte; (I) ao menos uma das variáveis é estatisticamente significante, pois o teste F obteve um p valor abaixo de 0.01. Os todos os coeficientes foram estatísticamente significantes, com níveis diferentes.O modelo explica aproximadamente 52% das variações do ebitida.
2.1.1 Normalidade dos resíduos
Observar a comportamento dos resíduos, primeiro podemos inspecionar por meio de um gráfico.
Residuals vs Fitted
Os resíduos versos os estimados (fitted), pemite observar se a lineariedade se mantém, isso pode ser observado pela linha vermelhar próxima da linha tracejada. Se for homocedástico, a distribuição dever ser aproximadamente a mesma ao longo do eixo x. Também, é possível observar se existem resíduos outliers.,/p.
Exemplo de não lineariedade
Normal Q-Q plot
O gráfico QQ, ou gráfico quantil-quantil, é uma ferramenta gráfica para nos ajudar a avaliar se um conjunto de dados veio plausivelmente de alguma distribuição teórica, como Normal ou exponencial. Por exemplo, se executarmos uma análise estatística que assume que nossa variável dependente é normalmente distribuída, podemos usar um gráfico QQ Normal para verificar essa suposição. É apenas uma verificação visual, não uma prova hermética, por isso é um pouco subjetivo. Mas nos permite ver rapidamente se nossa suposição é plausível e, se não, como a suposição é violada e quais pontos de dados contribuem para a violação.
Note que na padronizada se tem uma distribuição com “caudas pesadas” frente uma distribuição Normal:
Observe que os pontos caem ao longo de uma linha no meio do gráfico, mas se curvam nas extremidades. Gráficos QQ normais que exibem esse comportamento geralmente significam que seus dados têm valores mais extremos do que seria esperado se eles realmente viessem de uma distribuição Normal.
Scale-Location
Um gráfico de localização de escala é um tipo de gráfico que exibe os valores ajustados de um modelo de regressão ao longo do eixo x e a raiz quadrada dos resíduos padronizados ao longo do eixo y.
Verifique se a linha vermelha é aproximadamente horizontal ao longo do gráfico. Se for, então a suposição de homocedasticidade é provavelmente satisfeita para um determinado modelo de regressão. Ou seja, a dispersão dos resíduos é aproximadamente igual em todos os valores ajustados.
Verifique se não há um padrão claro entre os resíduos. Em outras palavras, os resíduos devem ser espalhados aleatoriamente ao redor da linha vermelha com variabilidade aproximadamente igual em todos os valores ajustados.
Mais em statology
Residual vs leverage
Não há uma única definição aceita para o que constitui um outlier. Uma definição possível é que um outlier é qualquer ponto que não é bem aproximado pelo modelo (tem um grande resíduo) e que influencia significativamente o ajuste do modelo (tem grande alavancagem). É aqui que entra o gráfico Residuais vs Alavancagem.
cook distance
a distância de Cook ou Cook’s D é uma estimativa comumente usada da influência de um ponto de dados ao realizar uma análise de regressão de mínimos quadrados . [1] Em uma análise prática de mínimos quadrados ordinários , a distância de Cook pode ser usada de várias maneiras: para indicar pontos de dados influentes que valem particularmente a pena verificar quanto à validade; ou para indicar regiões do espaço de projeto onde seria bom poder obter mais pontos de dados. É nomeado após o estatístico americano R. Dennis Cook , que introduziu o conceito em 1977.
par(mfrow = c(2, 2)) # Split the plotting panel into a 2 x 2 grid
plot(pooled)
cooksD <- cooks.distance(pooled)
pontoinfluentes <- cooksD[(cooksD > (3 * mean(cooksD, na.rm = TRUE)))]
pontoinfluentes## 54 214 225 362 533 587 670
## 0.57272991 0.03818052 0.28590258 0.70051657 0.65904161 0.04488703 0.69400371
## 830 841 895 911 978 1149 1203
## 0.03301919 0.25901454 0.04338981 0.04519834 1.22911695 2.28549946 0.04966824
## 1219 1286 1457 1527 1594 1765 1781
## 0.12952816 0.92470467 4.02425182 0.20402473 0.39977257 4.64943346 0.03349287
## 1835
## 0.91357610Podemos decidir remover os pontos influentes considerando a Cook Distance
#install.packages("dplyr")
library(dplyr)
ids_influentes <- names(pontoinfluentes)
outliers <- df[ids_influentes,]
outliers <- outliers %>% arrange(CompanyName,Ano,TRBCSector)
outliers## # A tibble: 22 × 10
## CompanyName TRBCSector Ano AtivoTotal EBITDA tamanho ROA marketShare
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Cogna Educaca… Academic … 2014 16639. 2039. 9.72 0.123 0.644
## 2 Cogna Educaca… Academic … 2015 17601. 2551. 9.78 0.145 0.677
## 3 Cogna Educaca… Academic … 2016 18220. 2105. 9.81 0.116 0.673
## 4 Cogna Educaca… Academic … 2017 31949. 1756. 10.4 0.0550 0.766
## 5 Cogna Educaca… Academic … 2018 34118. 2702. 10.4 0.0792 0.726
## 6 Cogna Educaca… Academic … 2019 30784. 1030. 10.3 0.0335 0.592
## 7 Oi SA em Recu… Technology 2014 103008. 7617. 10.7 0.0739 0.398
## 8 Oi SA em Recu… Technology 2016 68639. -1739. 10.7 -0.0253 0.314
## 9 Petroleo Bras… Energy 2014 900135 81302 10.7 0.0903 0.879
## 10 Petroleo Bras… Energy 2015 804945 91719 10.7 0.114 0.865
## # ℹ 12 more rows
## # ℹ 2 more variables: ihh <dbl>, concentracao <dbl>Podemos rodar um modelo sem os pontos influentes
df_cookDistance <- df %>% anti_join(outliers)## Joining with `by = join_by(CompanyName, TRBCSector, Ano, AtivoTotal, EBITDA,
## tamanho, ROA, marketShare, ihh, concentracao)`modelo_cookdistance <- lm(EBITDA~marketShare+ihh+tamanho ,data = df_cookDistance)
summary(modelo_cookdistance)##
## Call:
## lm(formula = EBITDA ~ marketShare + ihh + tamanho, data = df_cookDistance)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15193.1 -387.6 -21.2 338.0 25896.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.874e+03 2.017e+02 -9.288 <2e-16 ***
## marketShare 3.312e+04 1.088e+03 30.434 <2e-16 ***
## ihh 3.062e-02 2.610e-02 1.173 0.241
## tamanho 2.777e+02 2.707e+01 10.258 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1775 on 1822 degrees of freedom
## Multiple R-squared: 0.5255, Adjusted R-squared: 0.5247
## F-statistic: 672.7 on 3 and 1822 DF, p-value: < 2.2e-16Para evitar que os dados fiquem desbalanceados podemos optar por retirar as empresas que constam nos pontos influentes
valores_unicos <-outliers %>% distinct(factor(CompanyName))
valores_unicos## # A tibble: 6 × 1
## `factor(CompanyName)`
## <fct>
## 1 Cogna Educacao SA
## 2 Oi SA em Recuperacao Judicial
## 3 Petroleo Brasileiro SA Petrobras
## 4 Rede D'Or Sao Luiz SA
## 5 Telefonica Brasil SA
## 6 Vale SAdf_cookDistance_balanceado <- df %>%
filter(CompanyName != 'Cogna Educacao SA'&
CompanyName != 'Oi SA em Recuperacao Judicial'&
CompanyName != 'Petroleo Brasileiro SA Petrobras' &
CompanyName != "Rede D'Or Sao Luiz SA" &
CompanyName != "Telefonica Brasil SA" &
CompanyName != "Vale SA" )
modelo_cookdistance_bln <- lm(EBITDA~marketShare+ihh+tamanho,
data = df_cookDistance_balanceado)
summary(modelo_cookdistance_bln)##
## Call:
## lm(formula = EBITDA ~ marketShare + ihh + tamanho, data = df_cookDistance_balanceado)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15363.5 -392.8 -27.3 333.5 18845.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.894e+03 1.863e+02 -10.165 <2e-16 ***
## marketShare 3.374e+04 1.289e+03 26.169 <2e-16 ***
## ihh 3.402e-02 2.343e-02 1.452 0.147
## tamanho 2.802e+02 2.558e+01 10.955 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1591 on 1808 degrees of freedom
## Multiple R-squared: 0.5037, Adjusted R-squared: 0.5029
## F-statistic: 611.8 on 3 and 1808 DF, p-value: < 2.2e-16par(mfrow = c(2, 2)) # Split the plotting panel into a 2 x 2 grid
plot(modelo_cookdistance_bln)
Podemos observar que não há valores influentes acima de 1.
Vamos observar os pressupostos considerrando o modelo pooled
Para realizar o teste de Kolmogorov-Smirnov pode se usar o seguinte comando
residuos <- resid(pooled)
ks.test(residuos, "pnorm",mean(residuos), sd(residuos))##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: residuos
## D = 0.32305, p-value < 2.2e-16
## alternative hypothesis: two-sidedPara realizar o teste de Shapiro-Wilk pode se usar o seguinte comando
shapiro.test(residuos)##
## Shapiro-Wilk normality test
##
## data: residuos
## W = 0.30721, p-value < 2.2e-16Como podemos ver pelos gráficos e pelos testes de normalidades, se rejeita a hipótese de que os erros possuam distribuição normal. Vale ressaltar que o teste de Shapiro-Wilk é para amostras pequenas. Vamos chegar a heterocedasticidade
2.1.2 Teste heterocedasticidade
#===Teste heterocedasticidade
#install.packages("skedastic")
#library(skedastic)
#white_lm(pooled, interactions = TRUE)
#função white_lm não encontrada.
library(lmtest)## Carregando pacotes exigidos: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericbptest(pooled, ~marketShare+ihh+tamanho + marketShare*ihh +marketShare*tamanho + tamanho*ihh + I(marketShare^2)+I(ihh^2)+I(tamanho^2),data = df)##
## studentized Breusch-Pagan test
##
## data: pooled
## BP = 946.04, df = 9, p-value < 2.2e-16#lm(mpg~disp+hp, data=mtcars)
#bptest(pooled, ~ disp*hp + I(disp^2) + I(hp^2), data = mtcars)O resultado aponta que a covariancia dos resíduos não é constante, isto é, se tem a presença de heterocedasticidade
2.1.3 Teste Autocorrelação
#===Teste Autocorrelação durbin watson
#install.packages("lmtest")
library(lmtest)
dwtest(pooled) ##
## Durbin-Watson test
##
## data: pooled
## DW = 1.9809, p-value = 0.3402
## alternative hypothesis: true autocorrelation is greater than 0Ou um pouco mais detalhado
#install.packages("car")
library(car)## Carregando pacotes exigidos: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## somedurbinWatsonTest(pooled)## lag Autocorrelation D-W Statistic p-value
## 1 0.009552761 1.980881 0.548
## Alternative hypothesis: rho != 0Por outro lado, não tem auto correlação residual de primeira ordem
2.1.4 Teste de multicolineariedade
#===Teste multicolineariedade - Variance Inflation Factor VIF
car::vif(pooled)## marketShare ihh tamanho
## 1.318928 1.054356 1.265735
Ou então, pelo pacote plm()
vif(pooled)## marketShare ihh tamanho
## 1.318928 1.054356 1.265735Os valores de vif também estão dentro do esperado entre 1 e 10. Precisamos, corrigir os resíduos pela correção de White, erros padrão robustos
2.1.4.1 Erros padrão robustos
#install.packages("sandwich")
#===correção dos resíduos
library(sandwich)
coeftest(pooled, vcov = vcovHC(pooled, type = "HC1"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.8527e+02 8.9887e+02 0.3174 0.75100
## marketShare 6.8784e+04 1.3180e+04 5.2188 2.004e-07 ***
## ihh 1.7122e-01 7.0127e-02 2.4416 0.01472 *
## tamanho -1.3997e+02 1.6655e+02 -0.8404 0.40077
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Podemos observar que após a correção a variável tamanho da empresa deixou de ser estatisticamente significante a 10%
2.2 Efeito Fixo
2.2.1 Matriz dummy
#install.packages("plm")
library(plm)##
## Attaching package: 'plm'## The following objects are masked from 'package:dplyr':
##
## between, lag, leadDeclarar dados em painel
dfpanel <- pdata.frame(df, index=c("CompanyName", "Ano"),
drop.index=FALSE, row.names=TRUE)
pdim(dfpanel)## Balanced Panel: n = 308, T = 6, N = 1848Construir o modelo e ver os resultados
efeitofixo <- plm(EBITDA~marketShare+ihh+tamanho,
data = dfpanel,
effect = "individual",
index = c("CompanyName","Ano"),
model = "within")
summary(efeitofixo)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = EBITDA ~ marketShare + ihh + tamanho, data = dfpanel,
## effect = "individual", model = "within", index = c("CompanyName",
## "Ano"))
##
## Balanced Panel: n = 308, T = 6, N = 1848
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -37125.9187 -126.6958 0.7935 122.6558 52976.2212
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## marketShare -4656.00577 6778.03689 -0.6869 0.49223
## ihh -0.87885 0.42177 -2.0837 0.03735 *
## tamanho 479.99744 207.57941 2.3124 0.02089 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1.0242e+10
## Residual Sum of Squares: 1.0165e+10
## R-Squared: 0.0075172
## Adj. R-Squared: -0.19266
## F-statistic: 3.88049 on 3 and 1537 DF, p-value: 0.00888122.2.1.1 Obter a constante
A constante muda para cada empresa, para obter a constante média é
necessário utilizar o comando within_intercept() e para
observar todas as constantes fixef
#fx_level <- fixef(efeitofixo)
overallint <- within_intercept(efeitofixo,
vcov = NULL,
return.model = TRUE,)
summary(overallint)## Pooling Model
##
## Call:
## plm(formula = form, data = data, model = "pooling")
##
## Balanced Panel: n = 308, T = 6, N = 1848
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -37125.9187 -126.6958 0.7935 122.6558 52976.2212
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) -493.33927 1772.48814 -0.2783 0.78079
## marketShare -4656.00577 6778.03689 -0.6869 0.49222
## ihh -0.87885 0.42177 -2.0837 0.03732 *
## tamanho 479.99744 207.57941 2.3124 0.02087 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1.0242e+10
## Residual Sum of Squares: 1.0165e+10
## R-Squared: 0.0075172
## Adj. R-Squared: 0.0059026
## F-statistic: 3.88049 on 3 and 1844 DF, p-value: 0.00885392.2.1.2 Last Squared Dummy Variables
Model_lsdv <- lm(data = df, EBITDA ~ marketShare + ihh + tamanho+ factor(CompanyName))
resumo_lsdv <- summary(Model_lsdv)
resumo_lsdv$coefficients[1:4]## [1] 3576.4701688 -4656.0057735 -0.8788451 479.9974407Para evitar um resultado muito longo, este comando traz apenas os quatro primeiros coeficientes. Se pode observar que o valor da constante mudou, isso ocorreu porque não se trata mais de um beta médio. Agora caso seja necessário saber os valores das constantes de cada empresa, é necessário somar o coeficiente da respectiva dummy e somar ao intercepto
Seguindo com os pressupostos do modelo com a matriz dummy
2.2.1.3 Vif
vif(lm(EBITDA~marketShare+ihh+tamanho+factor(Ano)+factor(CompanyName),data = df))## GVIF Df GVIF^(1/(2*Df))
## marketShare 8.845914e+01 1 9.405272
## ihh 1.453853e+02 1 12.057582
## tamanho 5.640478e+01 1 7.510312
## factor(Ano) 1.496877e+00 5 1.041163
## factor(CompanyName) 4.790083e+05 307 1.0215312.2.1.4 Testando efeito do tempo
efeitoFixo_tempo <- plm(EBITDA~marketShare+ihh+tamanho+factor(Ano),
data = dfpanel,
index = c("CompanyName","Ano"),
model = "within")
summary(efeitoFixo_tempo)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = EBITDA ~ marketShare + ihh + tamanho + factor(Ano),
## data = dfpanel, model = "within", index = c("CompanyName",
## "Ano"))
##
## Balanced Panel: n = 308, T = 6, N = 1848
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -36730.416 -317.404 31.621 298.661 52390.597
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## marketShare 1403.91832 6824.78086 0.2057 0.8370456
## ihh -0.50670 0.43733 -1.1586 0.2467937
## tamanho -128.96718 237.62401 -0.5427 0.5873903
## factor(Ano)2015 181.65559 205.74472 0.8829 0.3774194
## factor(Ano)2016 235.03847 207.03643 1.1353 0.2564474
## factor(Ano)2017 506.70042 206.46709 2.4541 0.0142326 *
## factor(Ano)2018 788.97484 212.46778 3.7134 0.0002118 ***
## factor(Ano)2019 1114.42589 227.75001 4.8932 1.097e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1.0242e+10
## Residual Sum of Squares: 9964100000
## R-Squared: 0.027094
## Adj. R-Squared: -0.17295
## F-statistic: 5.33306 on 8 and 1532 DF, p-value: 1.2619e-062.2.2 Efeito tempo, teste F
Pode fazer o teste F para testar a regressão com o efeito fixo tempo e sem
pFtest(efeitoFixo_tempo, efeitofixo)##
## F test for individual effects
##
## data: EBITDA ~ marketShare + ihh + tamanho + factor(Ano)
## F = 6.1655, df1 = 5, df2 = 1532, p-value = 1.149e-05
## alternative hypothesis: significant effects2.2.3 Teste F de Chow
pFtest(efeitofixo, pooled)##
## F test for individual effects
##
## data: EBITDA ~ marketShare + ihh + tamanho
## F = 20.748, df1 = 307, df2 = 1537, p-value < 2.2e-16
## alternative hypothesis: significant effectsSe o valor de p for < 0,05, então o modelo de efeitos fixos é a melhor escolha
2.2.4 Breusch Pagan Heterocedasticidade
bptest(EBITDA~marketShare+ihh+tamanho, data = dfpanel, studentize=F)##
## Breusch-Pagan test
##
## data: EBITDA ~ marketShare + ihh + tamanho
## BP = 48841, df = 3, p-value < 2.2e-16Como o valor p < 0,05, detectamos heterocedasticidade
2.2.5 matriz de covariância robusta - White
Teste de heterocedasticidade, correção de White. Os da função vcovHC: Robust Covariance Matrix Estimators argumentos são os seguintes
vcovHC(
x,
method = c("arellano", "white1", "white2"),
type = c("HC0", "sss", "HC1", "HC2", "HC3", "HC4"),
cluster = c("group", "time"))
No entanto, para replicação do STATA é necessário escolher
type = "sss", pois faz uma correção para pequenas amostras
efeitofixoRobusto <- coeftest(efeitofixo,
vcovHC(efeitofixo, type = "sss"))
efeitofixoRobusto##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## marketShare -4656.00577 12900.76357 -0.3609 0.71822
## ihh -0.87885 0.37843 -2.3223 0.02035 *
## tamanho 479.99744 115.73741 4.1473 3.548e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1efeitofixoRobusto <- coeftest(efeitofixo,
vcovHC(efeitofixo, type = "HC1"))
efeitofixoRobusto##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## marketShare -4656.00577 12883.28996 -0.3614 0.71785
## ihh -0.87885 0.37792 -2.3255 0.02018 *
## tamanho 479.99744 115.58065 4.1529 3.463e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.2.6 correlação serial
pbgtest(efeitofixo)##
## Breusch-Godfrey/Wooldridge test for serial correlation in panel models
##
## data: EBITDA ~ marketShare + ihh + tamanho
## chisq = 467.59, df = 6, p-value < 2.2e-16
## alternative hypothesis: serial correlation in idiosyncratic errorsComo valor de p > 0,05, concluímos que há correlação serial
2.3 Efeito Aleatório
random <- plm(EBITDA~marketShare+ihh+tamanho,
data = dfpanel,
index = c("CompanyName","Ano"),
model = "random")
summary(random)## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = EBITDA ~ marketShare + ihh + tamanho, data = dfpanel,
## model = "random", index = c("CompanyName", "Ano"))
##
## Balanced Panel: n = 308, T = 6, N = 1848
##
## Effects:
## var std.dev share
## idiosyncratic 6613247 2572 0.235
## individual 21526062 4640 0.765
## theta: 0.7793
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -29875.394 -135.747 24.621 101.705 61218.090
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) -1247.18445 1059.28247 -1.1774 0.2390
## marketShare 52582.09897 3393.58575 15.4946 <2e-16 ***
## ihh 0.16015 0.16190 0.9891 0.3226
## tamanho 128.87182 132.34450 0.9738 0.3302
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1.5087e+10
## Residual Sum of Squares: 1.2847e+10
## R-Squared: 0.14849
## Adj. R-Squared: 0.1471
## Chisq: 321.557 on 3 DF, p-value: < 2.22e-162.3.1 LM multiplier de Breusch-Pagan
plmPooled <- plm( EBITDA ~ marketShare + ihh + tamanho,
data=dfpanel, model="pooling")
plmtest(plmPooled, type=c("bp"))##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: EBITDA ~ marketShare + ihh + tamanho
## chisq = 2581.4, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effectsO teste LM ajuda você a decidir entre uma regressão de efeitos aleatórios e uma regressão OLS simples.
A hipótese nula no teste LM é que as variâncias entre as entidades são zero. Isto é, nenhuma diferença significativa entre as unidades (ou seja, nenhum efeito de painel).
Neste caso a hipótese nula foi rejeitada, entre o modelo pooled e o aleatório, a melhor escolha é o aleatório
2.3.2 Teste de Hausmann
phtest(efeitofixo, random)##
## Hausman Test
##
## data: EBITDA ~ marketShare + ihh + tamanho
## chisq = 103.21, df = 3, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistentP valor < 0.05 deve se escolher o modelo de efeito aleatório
3 Bônus
As vezes é necesário calcular o mesmo modelo em um sub conjunto de dados, e coletar ao menos os p valores e os respectivos coeficientes
O pacote broom resolve isso.
#install.packages("broom")
library(broom)Da seguinte forma.
df_subset <- dfpanel %>% select(Ano,EBITDA,marketShare,ihh,tamanho)
coef_pvlr <- df_subset %>%
nest(-Ano) %>%
mutate(fit = map(data, ~ lm(EBITDA ~ marketShare + ihh + tamanho, data = .)),
results = map(fit, tidy)) %>%
unnest(results) %>%
select(Ano,term, estimate, p.value) ## Warning: Supplying `...` without names was deprecated in tidyr 1.0.0.
## ℹ Please specify a name for each selection.
## ℹ Did you want `data = -Ano`?
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.coef_pvlr <- as.data.frame(coef_pvlr)
names(coef_pvlr) <- c("Ano", "Variavel", "Coeficientes","P valor")
coef_pvlr## Ano Variavel Coeficientes P valor
## 1 2014 (Intercept) -1.379113e+02 8.821893e-01
## 2 2014 marketShare 4.710207e+04 1.031696e-46
## 3 2014 ihh 2.082800e-01 8.425411e-02
## 4 2014 tamanho -7.045401e+01 5.644277e-01
## 5 2015 (Intercept) 2.335039e+02 8.179371e-01
## 6 2015 marketShare 5.868888e+04 1.555128e-54
## 7 2015 ihh 1.391598e-01 3.055366e-01
## 8 2015 tamanho -1.242277e+02 3.516653e-01
## 9 2016 (Intercept) -4.130209e+02 6.659368e-01
## 10 2016 marketShare 5.197037e+04 2.115394e-49
## 11 2016 ihh 1.624163e-01 2.198383e-01
## 12 2016 tamanho -1.556936e+01 9.003387e-01
## 13 2017 (Intercept) 4.681825e+02 7.382735e-01
## 14 2017 marketShare 6.913845e+04 3.258891e-44
## 15 2017 ihh 1.281742e-01 4.905155e-01
## 16 2017 tamanho -1.558273e+02 3.865405e-01
## 17 2018 (Intercept) 9.687039e+02 5.504810e-01
## 18 2018 marketShare 8.486290e+04 3.228978e-47
## 19 2018 ihh 1.672719e-01 4.428174e-01
## 20 2018 tamanho -2.441422e+02 2.353623e-01
## 21 2019 (Intercept) 1.840288e+03 3.517708e-01
## 22 2019 marketShare 1.052264e+05 6.051731e-49
## 23 2019 ihh 2.337586e-01 3.873521e-01
## 24 2019 tamanho -3.993276e+02 1.048198e-01Gráfico dos coeficientes
coefmarketshare <- coef_pvlr %>% filter( Variavel =="marketShare")
ggplot(coefmarketshare,aes(Ano, Coeficientes, group=1)) +
geom_line(color="Black")+
geom_point(color="Blue")
Podemos identificar que o coeficiente aumenta de acordo com o tempo, no entanto, é possível realizar inúmeras análises
4 Referência
Croissant, Y., & Millo, G. (2008). Panel Data Econometrics in R: The plm Package. Journal of Statistical Software, 27(2), 1–43. https://doi.org/10.18637/jss.v027.i02
Leppert, P. (2021).R Tutorial: Panel Data Analysis 1. RPubs