Exemplo de Regressão de dados em painel (livro análise de dados - Fávero) - 2021

Regressão com Dados em Painel

library(readxl)
DesempenhoPainelCurto <- read_excel("DesempenhoPainelCurto.xls", 
    col_types = c("text", "numeric", "text", 
        "text", "numeric", "numeric", "numeric"))
View(DesempenhoPainelCurto)

require(plm)

Carregando pacotes exigidos: plm

data=pdata.frame(DesempenhoPainelCurto, index=c("id","tempo"))
  
head(data)

summary(data)

  estudante               id         classe              tempo    
 Length:720         1      : 24   Length:720         1      : 30  
 Class :character   2      : 24   Class :character   10     : 30  
 Mode  :character   3      : 24   Mode  :character   11     : 30  
                    4      : 24                      12     : 30  
                    5      : 24                      13     : 30  
                    6      : 24                      14     : 30  
                    (Other):576                      (Other):540  
     nota_Y          horas_X1       faltas_X2    
 Min.   : 28.50   Min.   :12.80   Min.   : 0.00  
 1st Qu.: 41.97   1st Qu.:21.30   1st Qu.: 9.00  
 Median : 62.18   Median :24.10   Median :14.00  
 Mean   : 61.96   Mean   :24.02   Mean   :14.17  
 3rd Qu.: 81.87   3rd Qu.:26.60   3rd Qu.:19.00  
 Max.   :100.00   Max.   :37.30   Max.   :28.00  
                                                 

coplot(nota_Y ~ tempo|id, type="b", data=data)

Modelo Pooled

reg.pooled=plm(nota_Y~horas_X1+faltas_X2, 
               data=data, model="pooling")
summary(reg.pooled)

Pooling Model

Call:
plm(formula = nota_Y ~ horas_X1 + faltas_X2, data = data, model = "pooling")

Balanced Panel: n = 30, T = 24, N = 720

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-30.500493  -6.923118  -0.094251   7.188799  29.736521 

Coefficients:
             Estimate Std. Error  t-value Pr(>|t|)    
(Intercept) 96.796643   2.472157  39.1547   <2e-16 ***
horas_X1     0.057336   0.094426   0.6072   0.5439    
faltas_X2   -2.556470   0.056921 -44.9124   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    275280
Residual Sum of Squares: 71988
R-Squared:      0.73849
Adj. R-Squared: 0.73776
F-statistic: 1012.4 on 2 and 717 DF, p-value: < 2.22e-16

Modelo Efeitos Fixos

reg.ef=plm(nota_Y~horas_X1+faltas_X2, 
           data=data, model="within")
summary(reg.ef)

Oneway (individual) effect Within Model

Call:
plm(formula = nota_Y ~ horas_X1 + faltas_X2, data = data, model = "within")

Balanced Panel: n = 30, T = 24, N = 720

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-10.102434  -3.086259  -0.063843   2.777009  13.762533 

Coefficients:
           Estimate Std. Error t-value  Pr(>|t|)    
horas_X1   0.302644   0.041098   7.364 5.116e-13 ***
faltas_X2 -0.354595   0.108605  -3.265  0.001149 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    13620
Residual Sum of Squares: 12167
R-Squared:      0.10662
Adj. R-Squared: 0.06637
F-statistic: 41.0561 on 2 and 688 DF, p-value: < 2.22e-16

summary(fixef(reg.ef))

   Estimate Std. Error t-value  Pr(>|t|)    
1   80.0153     1.6330  48.998 < 2.2e-16 ***
2   80.2816     1.5512  51.754 < 2.2e-16 ***
3   80.7463     1.6722  48.288 < 2.2e-16 ***
4   79.8659     1.4909  53.570 < 2.2e-16 ***
5   81.3204     1.8304  44.429 < 2.2e-16 ***
6   79.2648     1.4376  55.135 < 2.2e-16 ***
7   81.4106     1.9055  42.724 < 2.2e-16 ***
8   79.1101     1.3806  57.303 < 2.2e-16 ***
9   82.0856     1.9758  41.545 < 2.2e-16 ***
10  80.9245     1.7466  46.334 < 2.2e-16 ***
11  59.0288     2.1572  27.363 < 2.2e-16 ***
12  59.6457     2.0710  28.800 < 2.2e-16 ***
13  60.0113     2.2380  26.815 < 2.2e-16 ***
14  59.4963     1.9803  30.044 < 2.2e-16 ***
15  61.4886     2.4307  25.297 < 2.2e-16 ***
16  58.0637     1.9182  30.270 < 2.2e-16 ***
17  60.7477     2.5144  24.160 < 2.2e-16 ***
18  58.5425     1.8105  32.335 < 2.2e-16 ***
19  61.8489     2.6075  23.720 < 2.2e-16 ***
20  60.0639     2.3157  25.938 < 2.2e-16 ***
21  38.4986     2.8207  13.648 < 2.2e-16 ***
22  38.1981     2.7191  14.048 < 2.2e-16 ***
23  39.0781     2.9184  13.390 < 2.2e-16 ***
24  37.9335     2.6222  14.466 < 2.2e-16 ***
25  39.5431     3.0851  12.817 < 2.2e-16 ***
26  37.2788     2.5094  14.855 < 2.2e-16 ***
27  40.3043     3.1862  12.650 < 2.2e-16 ***
28  37.3468     2.4277  15.384 < 2.2e-16 ***
29  39.6169     3.3045  11.989 < 2.2e-16 ***
30  39.5464     2.9965  13.198 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(lm(nota_Y~horas_X1+faltas_X2+as.factor(id), data=data))

Call:
lm(formula = nota_Y ~ horas_X1 + faltas_X2 + as.factor(id), data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.1024  -3.0863  -0.0638   2.7770  13.7625 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      80.0153     1.6330  48.998  < 2e-16 ***
horas_X1          0.3026     0.0411   7.364 5.12e-13 ***
faltas_X2        -0.3546     0.1086  -3.265  0.00115 ** 
as.factor(id)2    0.2663     1.2197   0.218  0.82721    
as.factor(id)3    0.7310     1.2186   0.600  0.54880    
as.factor(id)4   -0.1494     1.2349  -0.121  0.90375    
as.factor(id)5    1.3051     1.2556   1.039  0.29898    
as.factor(id)6   -0.7504     1.2593  -0.596  0.55143    
as.factor(id)7    1.3953     1.2872   1.084  0.27878    
as.factor(id)8   -0.9052     1.2941  -0.699  0.48449    
as.factor(id)9    2.0703     1.3260   1.561  0.11891    
as.factor(id)10   0.9092     1.2322   0.738  0.46083    
as.factor(id)11 -20.9865     1.4288 -14.688  < 2e-16 ***
as.factor(id)12 -20.3696     1.3747 -14.818  < 2e-16 ***
as.factor(id)13 -20.0039     1.4874 -13.449  < 2e-16 ***
as.factor(id)14 -20.5190     1.3264 -15.469  < 2e-16 ***
as.factor(id)15 -18.5267     1.6253 -11.399  < 2e-16 ***
as.factor(id)16 -21.9516     1.2886 -17.036  < 2e-16 ***
as.factor(id)17 -19.2676     1.6970 -11.354  < 2e-16 ***
as.factor(id)18 -21.4728     1.2546 -17.115  < 2e-16 ***
as.factor(id)19 -18.1664     1.7743 -10.239  < 2e-16 ***
as.factor(id)20 -19.9514     1.5501 -12.871  < 2e-16 ***
as.factor(id)21 -41.5167     1.9457 -21.338  < 2e-16 ***
as.factor(id)22 -41.8172     1.8601 -22.481  < 2e-16 ***
as.factor(id)23 -40.9372     2.0320 -20.146  < 2e-16 ***
as.factor(id)24 -42.0818     1.7786 -23.660  < 2e-16 ***
as.factor(id)25 -40.4722     2.1991 -18.404  < 2e-16 ***
as.factor(id)26 -42.7365     1.6957 -25.203  < 2e-16 ***
as.factor(id)27 -39.7110     2.2911 -17.333  < 2e-16 ***
as.factor(id)28 -42.6685     1.6245 -26.265  < 2e-16 ***
as.factor(id)29 -40.3984     2.3912 -16.895  < 2e-16 ***
as.factor(id)30 -40.4689     2.1130 -19.153  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.205 on 688 degrees of freedom
Multiple R-squared:  0.9558,    Adjusted R-squared:  0.9538 
F-statistic: 479.9 on 31 and 688 DF,  p-value: < 2.2e-16

Modelo Efeitos Aleatórios

reg.ea=plm(nota_Y~horas_X1+faltas_X2,
           data=data, model="random", 
           random.method = "walhus")
summary(reg.ea)

Oneway (individual) effect Random Effect Model 
   (Wallace-Hussain's transformation)

Call:
plm(formula = nota_Y ~ horas_X1 + faltas_X2, data = data, model = "random", 
    random.method = "walhus")

Balanced Panel: n = 30, T = 24, N = 720

Effects:
                 var std.dev share
idiosyncratic 28.191   5.309 0.282
individual    71.793   8.473 0.718
theta: 0.8731

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-13.014599  -3.272224   0.066829   3.351574  16.476822 

Coefficients:
            Estimate Std. Error z-value  Pr(>|z|)    
(Intercept) 68.88629    2.39295 28.7872 < 2.2e-16 ***
horas_X1     0.24884    0.04426  5.6223 1.884e-08 ***
faltas_X2   -0.91108    0.10273 -8.8685 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    17832
Residual Sum of Squares: 14938
R-Squared:      0.16226
Adj. R-Squared: 0.15992
Chisq: 138.872 on 2 DF, p-value: < 2.22e-16

Comparação e escolha dos modelos

Modelo Pooled x Modelo de Efeitos Fixos

require(plm)
pFtest(reg.ef,reg.pooled)

    F test for individual effects

data:  nota_Y ~ horas_X1 + faltas_X2
F = 116.64, df1 = 29, df2 = 688, p-value < 2.2e-16
alternative hypothesis: significant effects

Como o valor p é inferior a 0,05, o modelo de Efeitos Fixos é melhor do que o modelo Pooled.

Modelo Pooled x Modelo de Efeitos Aleatórios

plmtest(reg.pooled, type="bp")

    Lagrange Multiplier Test - (Breusch-Pagan) for balanced panels

data:  nota_Y ~ horas_X1 + faltas_X2
chisq = 4269.1, df = 1, p-value < 2.2e-16
alternative hypothesis: significant effects

Modelo Efeitos Fixos x Modelo de Efeitos Aleatórios

phtest(reg.ef,reg.ea)

    Hausman Test

data:  nota_Y ~ horas_X1 + faltas_X2
chisq = 249.5, df = 2, p-value < 2.2e-16
alternative hypothesis: one model is inconsistent

Alguns testes para os modelos

Testando dependência transversal (cross-sectional)

pcdtest(reg.ea, test="cd")

    Pesaran CD test for cross-sectional dependence in panels

data:  nota_Y ~ horas_X1 + faltas_X2
z = 63.82, p-value < 2.2e-16
alternative hypothesis: cross-sectional dependence

Normalidade dos resíduos

shapiro.test(reg.ea$residuals)

    Shapiro-Wilk normality test

data:  reg.ea$residuals
W = 0.9974, p-value = 0.3148

Homocedasticidade dos resíduos

library(lmtest)

Carregando pacotes exigidos: zoo

Attaching package: 㤼㸱zoo㤼㸲

The following objects are masked from 㤼㸱package:base㤼㸲:

    as.Date, as.Date.numeric

bptest(reg.ea)

    studentized Breusch-Pagan test

data:  reg.ea
BP = 0.062415, df = 2, p-value = 0.9693

Testando correlação serial

# teste Breusch-Godfrey/Wooldridge - EFEITOS ALEATÓRIOS
pbgtest(reg.ea) 

    Breusch-Godfrey/Wooldridge test for serial correlation in panel
    models

data:  nota_Y ~ horas_X1 + faltas_X2
chisq = 455.38, df = 24, p-value < 2.2e-16
alternative hypothesis: serial correlation in idiosyncratic errors

Teste para efeitos individuais ou de tempo

# teste Wooldridge - POOLED
pwtest(reg.pooled) 

    Wooldridge's test for unobserved individual effects

data:  formula
z = 3.6195, p-value = 0.0002952
alternative hypothesis: unobserved effect

pwtest(reg.pooled, effect = "time") 

    Wooldridge's test for unobserved time effects

data:  formula
z = 2.9425, p-value = 0.003255
alternative hypothesis: unobserved effect

Testando raízes unitárias

require(tseries)

Carregando pacotes exigidos: tseries
Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

    㤼㸱tseries㤼㸲 version: 0.10-47

    㤼㸱tseries㤼㸲 is a package for time series analysis and
    computational finance.

    See 㤼㸱library(help="tseries")㤼㸲 for details.

adf.test(data$nota_Y, k=2)

p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  data$nota_Y
Dickey-Fuller = -5.8363, Lag order = 2, p-value = 0.01
alternative hypothesis: stationary

Referências

Breusch, Trevor S. 1978. «Testing for autocorrelation in dynamic linear models». Australian Economic Papers 17 (31): 334–55.

Breusch, Trevor S, e Adrian R Pagan. 1979. «A simple test for heteroscedasticity and random coefficient variation». Econometrica: Journal of the Econometric Society, 1287–94.

Breusch, Trevor Stanley, e Adrian Rodney Pagan. 1980. «The Lagrange multiplier test and its applications to model specification in econometrics». The Review of Economic Studies 47 (1): 239–53.

Gujarati, Damodar N., e Down C Porter. 2011. Econometria básica. 5a ed. New York: Mc Graw Hill. https://doi.org/10.1126/science.1186874.

Hausman, Jerry A. 1978. «Specification tests in econometrics». Econometrica: Journal of the econometric society, 1251–71.

Pesaran, M Hashem. 2015. «Testing weak cross-sectional dependence in large panels». Econometric Reviews 34 (6-10): 1089–1117.

Wooldridge, Jeffrey M. 2010. Econometric analysis of cross section and panel data. MIT Press.