#パッケージ呼び出し

library(AER)
##  要求されたパッケージ car をロード中です
##  要求されたパッケージ carData をロード中です
##  要求されたパッケージ lmtest をロード中です
##  要求されたパッケージ zoo をロード中です
## 
##  次のパッケージを付け加えます: 'zoo'
##  以下のオブジェクトは 'package:base' からマスクされています:
## 
##     as.Date, as.Date.numeric
##  要求されたパッケージ sandwich をロード中です
##  要求されたパッケージ survival をロード中です
library(AlgDesign)
library(car)
library(censReg)
##  要求されたパッケージ maxLik をロード中です
##  要求されたパッケージ 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/
## 
## 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/
library(DescTools)
## 
##  次のパッケージを付け加えます: 'DescTools'
##  以下のオブジェクトは 'package:car' からマスクされています:
## 
##     Recode
library(estimatr)
library(fixest)
library(ggplot2)
library(haven)
library(marginaleffects)
library(MatchIt)
library(mfx)
##  要求されたパッケージ MASS をロード中です
##  要求されたパッケージ betareg をロード中です
library(modelsummary)
## 
##  次のパッケージを付け加えます: 'modelsummary'
##  以下のオブジェクトは 'package:DescTools' からマスクされています:
## 
##     Format, Mean, Median, N, SD, Var
library(openxlsx)
library(plm)
library(psych)
## 
##  次のパッケージを付け加えます: 'psych'
##  以下のオブジェクトは 'package:modelsummary' からマスクされています:
## 
##     SD
##  以下のオブジェクトは 'package:ggplot2' からマスクされています:
## 
##     %+%, alpha
##  以下のオブジェクトは 'package:DescTools' からマスクされています:
## 
##     AUC, ICC, SD
##  以下のオブジェクトは 'package:car' からマスクされています:
## 
##     logit
library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
library(survival)
library(tableone)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.2.1
## ✔ purrr     1.0.4     ✔ tidyr     1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ psych::%+%()     masks ggplot2::%+%()
## ✖ psych::alpha()   masks ggplot2::alpha()
## ✖ dplyr::between() masks plm::between()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks plm::lag(), stats::lag()
## ✖ dplyr::lead()    masks plm::lead()
## ✖ dplyr::recode()  masks car::recode()
## ✖ dplyr::select()  masks MASS::select()
## ✖ purrr::some()    masks car::some()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(wooldridge)
## 
##  次のパッケージを付け加えます: 'wooldridge' 
## 
##  以下のオブジェクトは 'package:MASS' からマスクされています:
## 
##     cement
library(dplyr)

#データ読み込み

dataf <- 
  openxlsx::read.xlsx("仮想回答データ1群300人_神橋編集.xlsx")

#変数作成

dataf <- dataf %>% dplyr::mutate(
  respondent = case_when(group == 0 ~ id,
                         group == 1 ~ 300 + id,
                         group == 2 ~ 600 + id,
                         group == 3 ~ 900 + id),
  STR = respondent*100 + set_id,
  message = case_when(group == 0 ~ "control",
                      group == 1 ~ "selfish",
                      group == 2 ~ "altruistic",
                      group == 3 ~ "norm"),
message = factor(message,
                   levels = c("control", "selfish", "altruistic", "norm")))

#clogit実行

clogout1 <- survival::clogit(
  choice ~ price + crowd + review + time + strata(STR),
  data = dataf
)

clogout1
## Call:
## survival::clogit(choice ~ price + crowd + review + time + strata(STR), 
##     data = dataf)
## 
##              coef  exp(coef)   se(coef)      z      p
## price  -1.772e-03  9.982e-01  3.528e-05 -50.22 <2e-16
## crowd  -1.446e+00  2.355e-01  3.305e-02 -43.75 <2e-16
## review  9.617e-01  2.616e+00  2.126e-02  45.24 <2e-16
## time   -1.865e-02  9.815e-01  8.268e-04 -22.55 <2e-16
## 
## Likelihood ratio test=10248  on 4 df, p=< 2.2e-16
## n= 28800, number of events= 14400
clogout1$loglik
## [1] -9981.319 -4857.356

#メッセージ交差項

clogout2 <- survival::clogit(
  choice ~ price + crowd + crowd:factor(message) + review + time + strata(STR),
  data = dataf)
clogout2
## Call:
## survival::clogit(choice ~ price + crowd + crowd:factor(message) + 
##     review + time + strata(STR), data = dataf)
## 
##                                       coef  exp(coef)   se(coef)       z
## price                           -1.795e-03  9.982e-01  3.566e-05 -50.352
## crowd                           -9.414e-01  3.901e-01  5.396e-02 -17.444
## review                           9.775e-01  2.658e+00  2.161e-02  45.223
## time                            -1.888e-02  9.813e-01  8.329e-04 -22.662
## crowd:factor(message)selfish    -7.281e-01  4.828e-01  8.129e-02  -8.958
## crowd:factor(message)altruistic -9.709e-01  3.788e-01  8.505e-02 -11.415
## crowd:factor(message)norm       -4.777e-01  6.202e-01  7.819e-02  -6.110
##                                        p
## price                            < 2e-16
## crowd                            < 2e-16
## review                           < 2e-16
## time                             < 2e-16
## crowd:factor(message)selfish     < 2e-16
## crowd:factor(message)altruistic  < 2e-16
## crowd:factor(message)norm       9.99e-10
## 
## Likelihood ratio test=10403  on 7 df, p=< 2.2e-16
## n= 28800, number of events= 14400
clogout2$loglik
## [1] -9981.319 -4779.833