govdoc-cr-analysis/mlm/readme_did_model_fit.R

library(dplyr)
library(lubridate)
library(rdd)

readme_df_filepath <- "/mmfs1/gscratch/comdata/users/mjilg/govdoc-cr-data/final_data/README_weekly_count_data.csv"
readme_df = read.csv(readme_df_filepath, header = TRUE)

window_num <- 5
readme_df <- readme_df |>
  filter(week_index >= (- window_num) & week_index <= (window_num)) |>
  mutate(scaled_age = scale(age)) |>
  mutate(scaled_age_at_commit = scale(age_at_commit))|>
  mutate(log1p_count = log1p(commit_count))

library(lme4)
library(optimx)
library(lattice)

all_gmodel <- glmer.nb(log1p_count ~ before_after * week_index + scaled_age + (before_after * week_index | project_id),
                       control=glmerControl(optimizer="bobyqa",
                                            optCtrl=list(maxfun=2e5)), nAGQ=0,
                       data=readme_df)

summary(all_gmodel)
#saveRDS(all_gmodel, "020125_readme_model.rda")

model_residuals <- residuals(all_gmodel)
acf(model_residuals)
vif(all_gmodel)

#getting between-group variance 
variance_components <- as.data.frame(VarCorr(all_gmodel))

#getting the BLUPs
library(broom.mixed)
library(ggplot2)
condvals <- broom.mixed::tidy(all_gmodel, effects = "ran_vals", conf.int = TRUE)
glmer_ranef_Dweek <- test_condvals [which(test_condvals $term == "before_after:week_index"),]
has_zero <- function(estimate, low, high){
  return(ifelse((low < 0),ifelse((high > 0), 1, 0), 2))
}
glmer_ranef_Dweek <- glmer_ranef_Dweek |>
  mutate(ranef_grouping = has_zero(estimate, conf.low, conf.high)) |>
  mutate(rank = rank(estimate))
g <- glmer_ranef_Dweek |> 
  ggplot(aes(x=rank, y=estimate, col = as.factor(ranef_grouping))) +
  geom_linerange(aes(ymin= conf.low, ymax= conf.high)) +
  theme_bw()
g

write.csv(glmer_ranef_Dweek, "0201_readme_dweek_ranefs.csv")