Ngày 4: Mô hình hồi quy tuyến tính cho dữ liệu panel

Việc 1: Sử dụng mô hình hồi quy tuyến tính để đánh giá xem tỷ lệ tội phạm (biến số kết cục robbery) có liên quan như thế nào đến hiện diện của cảnh sát (biến số tiên lượng police) từ dữ liệu tội phạm (Crime dataset reduced.csv)
Đọc dữ liệu

crime= read.csv("D:\\OneDrive\\Statistical courses\\Ton Duc Thang Uni_April2023\\Datasets\\Crime dataset reduced.csv")
names (crime)
## [1] "city"          "year"          "police"        "robbery"      
## [5] "population"    "blackpct"      "femheadpct"    "welfarepercap"
## [9] "edexppercap"

1.1. Nhìn chung tỷ lệ tội phạm có liên quan đến tỷ lệ cảnh sát không? a. Vẽ biểu đồ mô tả mối liên quan giữa tỷ lệ tội phạm và tỷ lệ cảnh sát

library(ggplot2)
ggplot(data=crime, aes(x=police, y=robbery)) + geom_point() + geom_smooth(method="lm", se=F) + theme(legend.position="none")
## `geom_smooth()` using formula = 'y ~ x'

  1. Thực hiện mô hình tuyến tính robbery = alpha + beta*police
summary (lm (robbery ~ police, data= crime))
## 
## Call:
## lm(formula = robbery ~ police, data = crime)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2701 -1.6965 -0.3829  1.9435  5.9325 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.4750     1.0514   3.305 0.001289 ** 
## police        1.7632     0.5094   3.461 0.000771 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.3 on 108 degrees of freedom
## Multiple R-squared:  0.09985,    Adjusted R-squared:  0.09152 
## F-statistic: 11.98 on 1 and 108 DF,  p-value: 0.0007714
  1. Bạn có nhận xét gì về mối liên quan này? Nhận xét này có vấn đề gì không?

1.2 Mối liên quan giữa tỷ lệ tội phạm và tỷ lệ cảnh sát có khác nhau ở từng thành phố không?
a. Vẽ biểu đồ mô tả mối liên quan giữa tỷ lệ tội phạm và tỷ lệ cảnh sát ở từng thành phố. Bạn có nhận xét gì khi so sánh với kết quả ở mục 1.1?

ggplot(data=crime, aes(x=police, y=robbery, col=city)) + geom_point(na.rm=T) + geom_smooth(method="lm", se=F) + theme(legend.position="none")
## `geom_smooth()` using formula = 'y ~ x'

  1. Đánh giá mối liên quan giữa tỷ lệ tội phạm và tỷ lệ cảnh sát cho từng thành phố
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(magrittr)
crime %>% group_by(city) %>% mutate(beta = lm(robbery ~ police)$coefficients[2]) %>% mutate(alpha = lm(robbery ~ police)$coefficients[1]) %>% summarize(mbeta=mean(beta), malpha=mean(alpha))
## # A tibble: 5 × 3
##   city     mbeta malpha
##   <chr>    <dbl>  <dbl>
## 1 fresno   -2.75   8.79
## 2 losangel -8.79  28.5 
## 3 oakland  -1.76  12.5 
## 4 sacramen -5.69  15.3 
## 5 sanfran  -1.83  13.2
  1. Bạn có nhận xét như thế nào về sử dụng mô hình hồi quy tuyến tính để đánh giá mối liên quan giữa tỷ lệ tội phạm và hiện diện của cảnh sát.

Việc 2: Nhận thấy rằng dữ liệu tội phạm là dữ liệu panel, bạn quyết định sử dụng mô hình hồi quy cho dữ liệu panel để đánh giá xem tỷ lệ tội phạm (biến số kết cục robbery) có liên quan như thế nào đến hiện diện của cảnh sát (biến số tiên lượng police)

2.1 Thực hiện mô hình naive

library (plm)
## Warning: package 'plm' was built under R version 4.2.3
## 
## Attaching package: 'plm'
## The following objects are masked from 'package:dplyr':
## 
##     between, lag, lead
summary (plm(robbery ~ police, model="pooling", data=crime))
## Pooling Model
## 
## Call:
## plm(formula = robbery ~ police, data = crime, model = "pooling")
## 
## Balanced Panel: n = 5, T = 22, N = 110
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -4.27006 -1.69655 -0.38288  1.94346  5.93250 
## 
## Coefficients:
##             Estimate Std. Error t-value  Pr(>|t|)    
## (Intercept)  3.47505    1.05140  3.3052 0.0012887 ** 
## police       1.76323    0.50943  3.4612 0.0007714 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    634.49
## Residual Sum of Squares: 571.13
## R-Squared:      0.09985
## Adj. R-Squared: 0.091515
## F-statistic: 11.98 on 1 and 108 DF, p-value: 0.00077136
  1. Diễn giải ý nghĩa của mô hình
  2. Mô hình naive có phù hợp không? Tại sao?

2.2 Thực hiện mô hình fixed effects và diễn giải kết quả

fixed = plm(robbery ~ police, model="within", index=c("city", "year"), data=crime)
summary (fixed)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = robbery ~ police, data = crime, model = "within", 
##     index = c("city", "year"))
## 
## Balanced Panel: n = 5, T = 22, N = 110
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -3.062674 -0.921336 -0.097664  0.902659  3.234792 
## 
## Coefficients:
##        Estimate Std. Error t-value  Pr(>|t|)    
## police -4.16010    0.74123 -5.6124 1.666e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    252.49
## Residual Sum of Squares: 193.79
## R-Squared:      0.23247
## Adj. R-Squared: 0.19557
## F-statistic: 31.4991 on 1 and 104 DF, p-value: 1.6658e-07

2.3 Thực hiện mô hình random effects và diễn giải kết quả

random = plm(robbery ~ police, model="random", index=c("city", "year"), random.method="swar", data=crime)
summary (random)
## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = robbery ~ police, data = crime, model = "random", 
##     random.method = "swar", index = c("city", "year"))
## 
## Balanced Panel: n = 5, T = 22, N = 110
## 
## Effects:
##                 var std.dev share
## idiosyncratic 1.863   1.365 0.349
## individual    3.470   1.863 0.651
## theta: 0.8456
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -2.42398 -1.04861 -0.29648  0.94161  3.56167 
## 
## Coefficients:
##             Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept) 13.90197    1.71220  8.1194 4.686e-16 ***
## police      -3.40244    0.72859 -4.6699 3.013e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    261.59
## Residual Sum of Squares: 217.64
## R-Squared:      0.168
## Adj. R-Squared: 0.1603
## Chisq: 21.8081 on 1 DF, p-value: 3.0132e-06

2.4 Mô hình fixed effect hay random effects phù hợp hơn? Tại sao?

phtest(fixed, random) 
## 
##  Hausman Test
## 
## data:  robbery ~ police
## chisq = 30.88, df = 1, p-value = 2.745e-08
## alternative hypothesis: one model is inconsistent

2.5 Mối liên quan giữa tỷ lệ tội phạm và hiện diện của cảnh sát có thay đổi theo thời gian không?
a. Kiểm định ảnh hưởng của thời gian lên mối liên quan giữa tỷ lệ tội phạm và hiện diện của cảnh sát

library (lme4)
## Loading required package: Matrix
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
model_time = lmer (robbery ~ year + police + year:police + (1 + year |city), data= crime)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0661239 (tol = 0.002, component 1)
summary (model_time)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: robbery ~ year + police + year:police + (1 + year | city)
##    Data: crime
## 
## REML criterion at convergence: 376.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2145 -0.6406 -0.0806  0.6926  2.7839 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  city     (Intercept) 125.39162 11.1978       
##           year          0.01525  0.1235  -0.99
##  Residual               1.34821  1.1611       
## Number of obs: 110, groups:  city, 5
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)  46.94488   16.82587  23.59176   2.790  0.01026 * 
## year         -0.47227    0.19660  19.22022  -2.402  0.02655 * 
## police      -23.92478    7.55536  22.93109  -3.167  0.00432 **
## year:police   0.28647    0.09025  18.46483   3.174  0.00513 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) year   police
## year        -0.993              
## police      -0.946  0.956       
## year:police  0.926 -0.950 -0.992
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0661239 (tol = 0.002, component 1)

  1. Nhận xét về mối liên quan giữa tỷ lệ tội phạm và hiện diện cảnh sát về (i) mức độ liên quan, (ii) giả định của mô hình

  2. Mối liên quan giữa tỷ lệ tội phạm và hiện diện của cảnh sát vào năm 1972 khác năm 1980 như thế nào?

Việc 3: Bạn hãy ghi lại tất cả những hàm/ lệnh trên trong RMarkdown và share trên mạng rpubs.com/taikhoancuaban