#パッケージ呼び出し
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