zoe_fergusion_stat_consult_oct_25_2022.R

library(geepack)    # Run install.packages("geepack") if necessary
library(tidyr)      # Run install.packages("tidyr") if necessary



### Generate synthetic data
### For demonstration purposes only.
### You do not need to understand how this section works.
### Please see line 43 for format of data frame and line 58 for fitting model.

set.seed(3)

participants_per_group <- 50
num_participants <- 2 * participants_per_group
id <- seq(1:num_participants)
GP <- runif(num_participants, 1, 7)
Treatment <- c(rep(0, participants_per_group), rep(1, participants_per_group))

# Perceived credibility for control group
PC_Defendant0 <- runif(participants_per_group, 4, 7)
PC_Plaintiff0 <- 8 - PC_Defendant0 + rnorm(participants_per_group, sd=0.5)
PC_Plaintiff0[PC_Plaintiff0 < 1] <- 1
PC_Plaintiff0[PC_Plaintiff0 > 7] <- 7

# Perceived credibility for treatment group
PC_Defendant1 <- runif(participants_per_group, 1, 4)
PC_Plaintiff1 <- 8 - PC_Defendant1 + rnorm(participants_per_group, sd=0.5)
PC_Plaintiff1[PC_Plaintiff1 < 1] <- 1
PC_Plaintiff1[PC_Plaintiff1 > 7] <- 7

PC_Defendant <- c(PC_Defendant0, PC_Defendant1)
PC_Plaintiff <- c(PC_Plaintiff0, PC_Plaintiff1)

# Create dataframe in desired format
df_wide <- data.frame(id, GP, Treatment, PC_Plaintiff, PC_Defendant)
df <- df_wide %>% pivot_longer(c(PC_Defendant, PC_Plaintiff),
                               values_to="PC",
                               names_to="IsPlaintiff")
df[which(df$IsPlaintiff == "PC_Defendant"), ]$IsPlaintiff <- "0"
df[which(df$IsPlaintiff == "PC_Plaintiff"), ]$IsPlaintiff <- "1"
df$IsPlaintiff <- as.numeric(df$IsPlaintiff)

head(df)
# > head(df)
# # A tibble: 6 × 5
#      id    GP Treatment IsPlaintiff    PC
#   <int> <dbl>     <dbl>       <dbl> <dbl>
# 1     1  2.01         0           0  6.30
# 2     1  2.01         0           1  2.07
# 3     2  5.85         0           0  6.05
# 4     2  5.85         0           1  2.56
# 5     3  3.31         0           0  4.63
# 6     3  3.31         0           1  3.56


### Run GEE

model <- geeglm(PC ~ GP * (Treatment + IsPlaintiff), id=as.factor(id), corstr="exchangeable", data=df)


coef(summary(model))
# > coef(summary(model))
#                   Estimate    Std.err        Wald   Pr(>|W|)
# (Intercept)     4.56603474 0.43199386 111.7180310 0.00000000
# GP             -0.15730560 0.09715288   2.6216638 0.10541386
# Treatment       0.06748147 0.12014930   0.3154472 0.57435665
# IsPlaintiff    -1.35533447 0.85983554   2.4846326 0.11496319
# GP:Treatment   -0.02169462 0.02627627   0.6816743 0.40901077
# GP:IsPlaintiff  0.38709105 0.19454803   3.9588821 0.04662466