Panel 資料分析

資料觀察

rm(list=ls())
library(readr)
fatality <- read_csv("https://raw.githubusercontent.com/tpemartin/Econometric-Analysis/master/Part%20II/fatality.csv")

載入Panel套件

library(plm)

Loading required package: Formula

宣告資料為Panel data frame

fatality<-pdata.frame(fatality,c("state","year"))

X軸為beertax Y軸為mrall

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:plm’:

    between, lag, lead

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(ggplot2)
fatality %>% ggplot(aes(x=beertax,y=mrall))+
  geom_point()

X軸為mrall Y軸為beertax

fatality %>% ggplot(aes(y=beertax,x=mrall))+
  geom_point()

不同州用不同顏色畫離散圖

library(ggplot2)
ggplot(data=fatality,aes(x=beertax,y=mrall,color=state))+
  geom_point()->f1
f1

不同年用不同顏色畫離散圖

library(ggplot2)
ggplot(data=fatality,aes(x=beertax,y=mrall,color=year))+
  geom_point()

去除每個州的中間點，即進行Demean

fatality$mrall_demean<-Within(fatality$mrall,effect=c('individual'))
fatality$beertax_demean<-Within(fatality$beertax,effect=c('individual'))

Demean 之後再畫一次離散圖

ggplot(data=fatality,aes(x=beertax_demean,y=mrall_demean,color=year))+
  geom_point()->f2
f2

模型估計

迴歸模型設定

model<-mrall~beertax

使用Pooled OLS

pool1<-plm(model, data=fatality, model='pooling')
summary(pool1)

Pooling Model

Call:
plm(formula = model, data = fatality, model = "pooling")

Balanced Panel: n=48, T=7, N=336

Residuals :
     Min.   1st Qu.    Median   3rd Qu.      Max. 
-1.09e-04 -3.78e-05 -9.44e-06  2.85e-05  2.28e-04 

Coefficients :
              Estimate Std. Error t-value  Pr(>|t|)    
(Intercept) 1.8533e-04 4.3567e-06 42.5391 < 2.2e-16 ***
beertax     3.6461e-05 6.2170e-06  5.8647 1.082e-08 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    1.0892e-06
Residual Sum of Squares: 9.8747e-07
R-Squared:      0.093363
Adj. R-Squared: 0.090648
F-statistic: 34.3943 on 1 and 334 DF, p-value: 1.0822e-08

f1

str(pool1)

List of 10
 $ coefficients: Named num [1:2] 1.85e-04 3.65e-05
  ..- attr(*, "names")= chr [1:2] "(Intercept)" "beertax"
 $ vcov        : num [1:2, 1:2] 1.90e-11 -1.98e-11 -1.98e-11 3.87e-11
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2] "(Intercept)" "beertax"
  .. ..$ : chr [1:2] "(Intercept)" "beertax"
 $ residuals   : Named num [1:336] -2.86e-05 -1.57e-05 -1.42e-05 -2.62e-05 2.29e-05 ...
  ..- attr(*, "names")= chr [1:336] "1-1982" "1-1983" "1-1984" "1-1985" ...
 $ df.residual : int 334
 $ formula     :'Formula' with 1 left-hand and 1 right-hand side: mrall ~ beertax
  ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
 $ model       :Classes ‘pdata.frame’ and 'data.frame': 336 obs. of  2 variables:
  ..$ mrall  :Classes 'pseries', 'numeric'  atomic [1:336] 0.000213 0.000235 0.000234 0.000219 0.000267 ...
  .. .. ..- attr(*, "index")=Classes ‘pindex’ and 'data.frame': 336 obs. of  2 variables:
  .. .. .. ..$ state: Factor w/ 48 levels "1","4","5","6",..: 1 1 1 1 1 1 1 2 2 2 ...
  .. .. .. ..$ year : Factor w/ 7 levels "1982","1983",..: 1 2 3 4 5 6 7 1 2 3 ...
  ..$ beertax:Classes 'pseries', 'numeric'  atomic [1:336] 1.54 1.79 1.71 1.65 1.61 ...
  .. .. ..- attr(*, "index")=Classes ‘pindex’ and 'data.frame': 336 obs. of  2 variables:
  .. .. .. ..$ state: Factor w/ 48 levels "1","4","5","6",..: 1 1 1 1 1 1 1 2 2 2 ...
  .. .. .. ..$ year : Factor w/ 7 levels "1982","1983",..: 1 2 3 4 5 6 7 1 2 3 ...
  ..- attr(*, "terms")=Classes 'terms', 'formula'  language mrall ~ beertax
  .. .. ..- attr(*, "variables")= language list(mrall, beertax)
  .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "mrall" "beertax"
  .. .. .. .. ..$ : chr "beertax"
  .. .. ..- attr(*, "term.labels")= chr "beertax"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(mrall, beertax)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. .. ..- attr(*, "names")= chr [1:2] "mrall" "beertax"
  ..- attr(*, "index")=Classes ‘pindex’ and 'data.frame':   336 obs. of  2 variables:
  .. ..$ state: Factor w/ 48 levels "1","4","5","6",..: 1 1 1 1 1 1 1 2 2 2 ...
  .. ..$ year : Factor w/ 7 levels "1982","1983",..: 1 2 3 4 5 6 7 1 2 3 ...
 $ assign      : int [1:2] 0 1
 $ args        :List of 5
  ..$ model        : chr "pooling"
  ..$ effect       : chr "individual"
  ..$ random.method: chr "swar"
  ..$ random.dfcor : NULL
  ..$ inst.method  : chr "bvk"
 $ aliased     : Named logi [1:2] FALSE FALSE
  ..- attr(*, "names")= chr [1:2] "(Intercept)" "beertax"
 $ call        : language plm(formula = model, data = fatality, model = "pooling")
 - attr(*, "class")= chr [1:2] "plm" "panelmodel"

pool1$coefficients

 (Intercept)      beertax 
1.853308e-04 3.646054e-05 

f1+stat_function(fun=function(x){
  pool1$coefficients[1]+pool1$coefficients[2]*x
  },color="black")->f.ols
f.ols

使用random effect

re1<-plm(model, data=fatality, model='random')
summary(re1)

Oneway (individual) effect Random Effect Model 
   (Swamy-Arora's transformation)

Call:
plm(formula = model, data = fatality, model = "random")

Balanced Panel: n=48, T=7, N=336

Effects:
                    var   std.dev share
idiosyncratic 3.605e-10 1.899e-05 0.119
individual    2.660e-09 5.158e-05 0.881
theta:  0.8622  

Residuals :
     Min.   1st Qu.    Median   3rd Qu.      Max. 
-4.71e-05 -1.20e-05 -2.15e-06  9.10e-06  9.64e-05 

Coefficients :
               Estimate  Std. Error t-value Pr(>|t|)    
(Intercept)  2.0671e-04  9.9971e-06 20.6773   <2e-16 ***
beertax     -5.2016e-06  1.2418e-05 -0.4189   0.6756    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    1.2648e-07
Residual Sum of Squares: 1.2642e-07
R-Squared:      0.00052508
Adj. R-Squared: -0.0024674
F-statistic: 0.175467 on 1 and 334 DF, p-value: 0.67557

f.ols+stat_function(fun=function(x){
  re1$coefficients[1]+re1$coefficients[2]*x
  },color="black",linetype="dashed")->f.ols.re
f.ols.re

使用fixed effects (or called within) models: Oneway (individual) effect

model

mrall ~ beertax

fe1<-plm(model, data=fatality, model='within', effect='individual')
summary(fe1)

Oneway (individual) effect Within Model

Call:
plm(formula = model, data = fatality, effect = "individual", 
    model = "within")

Balanced Panel: n=48, T=7, N=336

Residuals :
     Min.   1st Qu.    Median   3rd Qu.      Max. 
-5.87e-05 -8.28e-06 -1.27e-07  7.95e-06  8.98e-05 

Coefficients :
           Estimate  Std. Error t-value Pr(>|t|)    
beertax -6.5587e-05  1.8785e-05 -3.4915 0.000556 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    1.0785e-07
Residual Sum of Squares: 1.0345e-07
R-Squared:      0.040745
Adj. R-Squared: -0.11969
F-statistic: 12.1904 on 1 and 287 DF, p-value: 0.00055597

f2+stat_function(fun=function(x){
  fe1$coefficients*x
  },color="blue")->f.fe1
f.fe1

使用fixed effects (or called within) models: 用Twoways (individual and time) effect

fe2<-plm(model, data=fatality, model='within', effect='twoways')
summary(fe2)

Twoways effects Within Model

Call:
plm(formula = model, data = fatality, effect = "twoways", model = "within")

Balanced Panel: n=48, T=7, N=336

Residuals :
     Min.   1st Qu.    Median   3rd Qu.      Max. 
-5.96e-05 -8.10e-06  1.43e-07  8.23e-06  8.39e-05 

Coefficients :
           Estimate  Std. Error t-value Pr(>|t|)   
beertax -6.3998e-05  1.9738e-05 -3.2424 0.001328 **
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    1.029e-07
Residual Sum of Squares: 9.9193e-08
R-Squared:      0.036065
Adj. R-Squared: -0.14918
F-statistic: 10.5133 on 1 and 281 DF, p-value: 0.001328

f.fe1+stat_function(fun=function(x){
  fe2$coefficients*x
  },color="blue",linetype="dashed")->f.fe1.fe2
f.fe1.fe2

library(stargazer)

Please cite as: 

 Hlavac, Marek (2015). stargazer: Well-Formatted Regression and Summary Statistics Tables.
 R package version 5.2. http://CRAN.R-project.org/package=stargazer 

re1

Model Formula: mrall ~ beertax

Coefficients:
(Intercept)     beertax 
 2.0671e-04 -5.2016e-06 

stargazer(pool1,re1,fe1,fe2,type='text')

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

========================================================================================================
                                                 Dependent variable:                                    
             -------------------------------------------------------------------------------------------
                                                        mrall                                           
                       (1)                   (2)                   (3)                     (4)          
--------------------------------------------------------------------------------------------------------
beertax            0.00004***             -0.00001             -0.0001***              -0.0001***       
                    (0.00001)             (0.00001)             (0.00002)               (0.00002)       
                                                                                                        
Constant            0.0002***             0.0002***                                                     
                    (0.00000)             (0.00001)                                                     
                                                                                                        
--------------------------------------------------------------------------------------------------------
Observations           336                   336                   336                     336          
R2                    0.093                 0.001                 0.041                   0.036         
Adjusted R2           0.091                -0.002                -0.120                  -0.149         
F Statistic  34.394*** (df = 1; 334) 0.175 (df = 1; 334) 12.190*** (df = 1; 287) 10.513*** (df = 1; 281)
========================================================================================================
Note:                                                                        *p<0.1; **p<0.05; ***p<0.01

相關檢定

Hausman test：c_i是否與x_{it}有關

phtest(fe2,re1)

    Hausman Test

data:  model
chisq = 14.687, df = 1, p-value = 0.0001269
alternative hypothesis: one model is inconsistent

BP檢定：是否存在州的隨機效果項c_i

pwtest(pool1)

    Wooldridge's test for unobserved individual
    effects

data:  formula
z = 3.2727, p-value = 0.001065
alternative hypothesis: unobserved effect

White Test for heteroscedasticity

library(lmtest)
lm1<-lm(mrall~beertax,data=fatality)
bptest(lm1,varformula=~beertax+I(beertax^2), data=fatality)

    studentized Breusch-Pagan test

data:  lm1
BP = 15.895, df = 2, p-value = 0.0003536