L6_PanelLogit

library(plm)
library(pglm)
Warning: package 'pglm' was built under R version 4.5.1
Loading required package: maxLik
Loading required package: miscTools

Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.

If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
library(tidyverse)
Warning: package 'stringr' was built under R version 4.5.1
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.4     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::between() masks plm::between()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks plm::lag(), stats::lag()
✖ dplyr::lead()    masks plm::lead()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(texreg) 
Warning: package 'texreg' was built under R version 4.5.1
Version:  1.39.4
Date:     2024-07-23
Author:   Philip Leifeld (University of Manchester)

Consider submitting praise using the praise or praise_interactive functions.
Please cite the JSS article in your publications -- see citation("texreg").

Attaching package: 'texreg'

The following object is masked from 'package:tidyr':

    extract

Panel Logit

Datplogit <- haven::read_dta("C:\\Users\\Huynh Chuong\\Desktop\\University\\UEL\\Class_QuantMethods\\Giáo trình\\Data\\mus\\mus18data.dta")
names(Datplogit)
 [1] "plan"     "site"     "coins"    "tookphys" "year"     "id"      
 [7] "black"    "income"   "age"      "female"   "educdec"  "time"    
[13] "outpdol"  "drugdol"  "suppdol"  "mentdol"  "inpdol"   "med"     
[19] "totadm"   "inpmis"   "mentvis"  "mdu"      "notmdvis" "num"     
[25] "mhi"      "ndisease" "physlim"  "ghindx"   "mdeoff"   "pioff"   
[31] "child"    "femchild" "lfam"     "lpi"      "idp"      "lcoins"  
[37] "fmde"     "hlthg"    "hlthf"    "hlthp"    "xghindx"  "linc"    
[43] "lnum"     "lnmed"    "dmed"     "dmdu"    
table(Datplogit$dmdu)

    0     1 
 6308 13878 
table(Datplogit$dmdu, Datplogit$year)
   
       1    2    3    4    5
  0 1729 1832 1790  480  477
  1 3909 3742 3755 1235 1237
Datplogit %>% group_by(id) %>%
  summarise(trungbinh=mean(dmdu)) %>%
  head(n=20)
# A tibble: 20 × 2
       id trungbinh
    <dbl>     <dbl>
 1 125024     0.2  
 2 125025     0.2  
 3 125026     0.2  
 4 125027     0.4  
 5 125057     0.2  
 6 125058     0.8  
 7 125059     0    
 8 125075     0    
 9 125126     0    
10 125127     0.333
11 125128     0    
12 125129     0.667
13 125133     0.667
14 125134     1    
15 125135     0.667
16 125152     1    
17 125153     1    
18 125159     0.6  
19 125160     0.4  
20 125161     1    
Datplogit <- pdata.frame(Datplogit, index=c("id", "year"))
Datplogit <- Datplogit %>%
  mutate(Y=as_factor(dmdu))
pdim(Datplogit)
Unbalanced Panel: n = 5908, T = 1-5, N = 20186

Thiết lập mô hình logit

model<- Y~lcoins+ ndisease

Ước lượng mô hình Logit, Probit (POOLED)

M1_logit <- glm(model, Datplogit,family = binomial(link = 'logit'))
M1_probit <- glm(model, Datplogit,family = binomial(link = 'probit'))

Ước lượng mô hình Logit, Probit (Random effect)

plogit<-pglm(Y~lcoins+ ndisease , Datplogit, family = binomial(link = 'logit'), model="random" )
pprobit<- update(plogit, family = binomial(link = 'probit'))
screenreg(list(logit = M1_logit , probit = M1_probit, plogit = plogit, pprobit = pprobit ),digits = 3)

==============================================================================
                logit           probit          plogit          pprobit       
------------------------------------------------------------------------------
(Intercept)          0.557 ***       0.349 ***       0.845 ***       0.492 ***
                    (0.037)         (0.022)         (0.078)         (0.045)   
lcoins              -0.153 ***      -0.091 ***      -0.236 ***      -0.136 ***
                    (0.008)         (0.005)         (0.016)         (0.009)   
ndisease             0.057 ***       0.034 ***       0.090 ***       0.052 ***
                    (0.003)         (0.002)         (0.005)         (0.003)   
sigma                                                1.872 ***       1.090 ***
                                                    (0.046)         (0.026)   
------------------------------------------------------------------------------
AIC              24110.136       24117.746       21847.146       21842.321    
BIC              24133.874       24141.484                                    
Log Likelihood  -12052.068      -12055.873      -10919.573      -10917.161    
Deviance         24104.136       24111.746                                    
Num. obs.        20186           20186                                        
==============================================================================
*** p < 0.001; ** p < 0.01; * p < 0.05
# Newton-Raphson
plogit_NR<- pglm(Y~lcoins+ ndisease , Datplogit, family = binomial(link = 'logit'), control = list(method = "nr"))  
# Phương pháp ước lượng hỗn hợp
logit_mixed <- lme4::glmer(Y~lcoins+ ndisease + (1 | id), 
                     data = Datplogit, 
                 family = binomial)
screenreg(list(logit = M1_logit , probit = M1_probit, 
               plogit = plogit, pprobit = pprobit, NewtonRaphson=plogit_NR, mixed_model=logit_mixed ),
          digits = 3)

===================================================================================================================
                     logit           probit          plogit          pprobit         NewtonRaphson   mixed_model   
-------------------------------------------------------------------------------------------------------------------
(Intercept)               0.557 ***       0.349 ***       0.845 ***       0.492 ***       0.845 ***       0.848 ***
                         (0.037)         (0.022)         (0.078)         (0.045)         (0.078)         (0.075)   
lcoins                   -0.153 ***      -0.091 ***      -0.236 ***      -0.136 ***      -0.236 ***      -0.234 ***
                         (0.008)         (0.005)         (0.016)         (0.009)         (0.016)         (0.016)   
ndisease                  0.057 ***       0.034 ***       0.090 ***       0.052 ***       0.090 ***       0.089 ***
                         (0.003)         (0.002)         (0.005)         (0.003)         (0.005)         (0.005)   
sigma                                                     1.872 ***       1.090 ***       1.872 ***                
                                                         (0.046)         (0.026)         (0.046)                   
-------------------------------------------------------------------------------------------------------------------
AIC                   24110.136       24117.746       21847.146       21842.321       21847.146       22051.476    
BIC                   24133.874       24141.484                                                       22083.127    
Log Likelihood       -12052.068      -12055.873      -10919.573      -10917.161      -10919.573      -11021.738    
Deviance              24104.136       24111.746                                                                    
Num. obs.             20186           20186                                                           20186        
Num. groups: id                                                                                        5908        
Var: id (Intercept)                                                                                       3.030    
===================================================================================================================
*** p < 0.001; ** p < 0.01; * p < 0.05

Tobit panel

library(censReg)
Warning: package 'censReg' was built under R version 4.5.1

Please cite the 'censReg' package as:
Henningsen, Arne (2017). censReg: Censored Regression (Tobit) Models. R package version 0.5. http://CRAN.R-Project.org/package=censReg.

If you have questions, suggestions, or comments regarding the 'censReg' package, please use a forum or 'tracker' at the R-Forge site of the 'sampleSelection' project:
https://r-forge.r-project.org/projects/sampleselection/
# Tải dữ liệu ví dụ
data( "Affairs", package = "AER" )
panelResult <- censReg( affairs ~ age + yearsmarried + religiousness +  occupation + rating, left  = 0, right=Inf  , data = Affairs,  )
summary( panelResult )

Call:
censReg(formula = affairs ~ age + yearsmarried + religiousness + 
    occupation + rating, left = 0, right = Inf, data = Affairs)

Observations:
         Total  Left-censored     Uncensored Right-censored 
           601            451            150              0 

Coefficients:
              Estimate Std. error t value  Pr(> t)    
(Intercept)    8.17420    2.74145   2.982  0.00287 ** 
age           -0.17933    0.07909  -2.267  0.02337 *  
yearsmarried   0.55414    0.13452   4.119 3.80e-05 ***
religiousness -1.68622    0.40375  -4.176 2.96e-05 ***
occupation     0.32605    0.25442   1.282  0.20001    
rating        -2.28497    0.40783  -5.603 2.11e-08 ***
logSigma       2.10986    0.06710  31.444  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Newton-Raphson maximisation, 7 iterations
Return code 1: gradient close to zero (gradtol)
Log-likelihood: -705.5762 on 7 Df
# Tải gói cần thiết
library(AER)
Warning: package 'AER' was built under R version 4.5.1
Loading required package: car
Warning: package 'car' was built under R version 4.5.1
Loading required package: carData
Warning: package 'carData' was built under R version 4.5.1

Attaching package: 'car'
The following object is masked from 'package:dplyr':

    recode
The following object is masked from 'package:purrr':

    some
Loading required package: lmtest
Loading required package: zoo

Attaching package: 'zoo'
The following objects are masked from 'package:base':

    as.Date, as.Date.numeric
Loading required package: sandwich
Loading required package: survival
library(plm)

# Giả sử bạn có dữ liệu panel
set.seed(123)
n <- 100  # Số lượng cá nhân
T <- 5    # Số lượng thời gian
id <- rep(1:n, each = T)
time <- rep(1:T, n)
X1 <- rnorm(n * T)
X2 <- rnorm(n * T)
epsilon <- rnorm(n * T, mean = 0, sd = 1)
Y <- pmax(0, 5 + 2 * X1 - 3 * X2 + epsilon)  # Biến phụ thuộc với giới hạn từ 0

# Tạo dataframe
Datapanel <- data.frame(id, time, Y, X1, X2)
# Chuyển đổi dữ liệu thành dạng panel
Datapanel <- pdata.frame(Datapanel, index = c("id", "time"))

# Mô hình Tobit
tobit_model <- tobit(Y ~ X1 + X2, data = Datapanel, left = 0)
# Tóm tắt kết quả
summary(tobit_model)

Call:
tobit(formula = Y ~ X1 + X2, left = 0, data = Datapanel)

Observations:
         Total  Left-censored     Uncensored Right-censored 
           500             50            450              0 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  5.00857    0.04604 108.792   <2e-16 ***
X1           2.07645    0.04766  43.566   <2e-16 ***
X2          -2.97534    0.04843 -61.433   <2e-16 ***
Log(scale)  -0.01364    0.03331  -0.409    0.682    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Scale: 0.9865 

Gaussian distribution
Number of Newton-Raphson Iterations: 9 
Log-likelihood: -655.3 on 4 Df
Wald-statistic:  5525 on 2 Df, p-value: < 2.22e-16 
tobit_model2 <- censReg(Y ~ X1 + X2, data = Datapanel, left = 0)
# Tóm tắt kết quả
summary(tobit_model2)

Call:
censReg(formula = Y ~ X1 + X2, left = 0, data = Datapanel)

Observations:
         Total  Left-censored     Uncensored Right-censored 
           500             50            450              0 

Coefficients:
            Estimate Std. error t value Pr(> t)    
(Intercept)  5.00897    0.04750 105.449  <2e-16 ***
X1           2.07470    0.04768  43.511  <2e-16 ***
X2          -2.97683    0.04850 -61.372  <2e-16 ***
logSigmaMu  -2.03724    0.99891  -2.039  0.0414 *  
logSigmaNu  -0.02240    0.03718  -0.602  0.5470    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Newton-Raphson maximisation, 8 iterations
Return code 8: successive function values within relative tolerance limit (reltol)
Log-likelihood: -655.1903 on 5 Df