---
title: "gee tables"
output:
  html_document: 
    toc: true
    toc_float: true
    code_folding: hide
---

```{r setup, message=FALSE, warning=FALSE, echo = FALSE}
library(tidyverse)
library(broom)
library(emmeans)
library(ggpubr)
library(patchwork)
library(scales)
library(readxl)
library(kableExtra)
library(gee)
library(geepack)
library(pracma)


knitr::opts_chunk$set(
  warning = FALSE,
  message = FALSE,
  fig.width = 8,
  fig.height = 6
)
theme_set(theme_bw() + theme(legend.position = "bottom"))
```


This file contains code to reproduce generalized estimating equation (GEE) models and tables 2 and 3 for the paper.

# Read Data 

```{r}
## load scalar features data
scalars = readRDS(here::here("data", "scalar_metrics_pupils.rds"))

```

# Gee models

```{r}
## Percent Constriction
min_gee <- geeglm(min_constriction ~ eye + prepost*user_type,
               data = scalars, 
               id = factor(subject_id), 
               family = gaussian,
               corstr = "independence")

## AUC
auc_gee <- geeglm(AUC ~ eye + prepost*user_type,
               data = scalars, 
               id = factor(subject_id), 
               family = gaussian,
               corstr = "independence")
```


# Group means 

```{R}
# Int + right + post + occasional + daily + post*occasional + post*daily

pre_non = c(1, 0, 0, 0, 0, 0, 0)
pre_occ = c(1, 0, 0, 1, 0, 0, 0)
pre_daily = c(1, 0, 0, 0, 1, 0, 0)

# Int + right + post + occasional + daily + post*occasional + post*daily
post_non = c(1, 0, 1, 0, 0, 0, 0)
post_occ = c(1, 0, 1, 1, 0, 1, 0)
post_daily = c(1, 0, 1, 0, 1, 0, 1)

## min. constriction means
coef_min <- summary(min_gee)$coefficients[,1] # coefficients from gee model
se_min <- summary(min_gee)$coefficients[,2] # Std. Error from gee model


## Group Means at pre event for table 3 - min. constriction
mean_pre_non_min = sum(pre_non*coef_min)
mean_pre_occ_min = sum(pre_occ*coef_min)
mean_pre_daily_min = sum(pre_daily*coef_min)

## SE at pre event for table 3 - min. constriction
se_pre_non_min = sum(pre_non*se_min)
se_pre_occ_min = sum(pre_occ*se_min)
se_pre_daily_min = sum(pre_daily*se_min)

mean_pre_non_min = str_c(round(mean_pre_non_min,2), 
                         " (", round(se_pre_non_min,2), ")")
mean_pre_occ_min = str_c(round(mean_pre_occ_min,2), 
                         " (", round(se_pre_occ_min,2), ")")
mean_pre_daily_min = str_c(round(mean_pre_daily_min,2), 
                         " (", round(se_pre_daily_min,2), ")")

## Group Means at post event for table 3 - min. constriction
mean_post_non_min = sum(post_non*coef_min)
mean_post_occ_min = sum(post_occ*coef_min)
mean_post_daily_min = sum(post_daily*coef_min)

## SE at post event for table 3 - min. constriction

se_post_non_min = sum(post_non*se_min)
se_post_occ_min = sum(post_occ*se_min)
se_post_daily_min = sum(post_daily*se_min)


mean_post_non_min = str_c(round(mean_post_non_min,2), 
                         " (", round(se_post_non_min,2), ")")
mean_post_occ_min = str_c(round(mean_post_occ_min,2), 
                         " (", round(se_post_occ_min,2), ")")
mean_post_daily_min = str_c(round(mean_post_daily_min,2), 
                         " (", round(se_post_daily_min,2), ")")
#######################################
############## auc means ################
#######################################

coef_auc <- summary(auc_gee)$coefficients[,1] # coefficients from gee model
se_auc <- summary(auc_gee)$coefficients[,2] # coefficients from gee model

## Group Means at pre event for table 3 - AUC
mean_pre_non_auc = sum(pre_non*coef_auc)
mean_pre_occ_auc = sum(pre_occ*coef_auc)
mean_pre_daily_auc = sum(pre_daily*coef_auc)

## Group SE at pre event for table 3 - AUC
se_pre_non_auc = sum(pre_non*se_auc)
se_pre_occ_auc = sum(pre_occ*se_auc)
se_pre_daily_auc = sum(pre_daily*se_auc)

mean_pre_non_auc = str_c(round(mean_pre_non_auc,2), 
                         " (", round(se_pre_non_auc,2), ")")
mean_pre_occ_auc = str_c(round(mean_pre_occ_auc,2), 
                         " (", round(se_pre_occ_auc,2), ")")
mean_pre_daily_auc = str_c(round(mean_pre_daily_auc,2), 
                         " (", round(se_pre_daily_auc,2), ")")

## Group Means at post event for table 3 -AUC
mean_post_non_auc = sum(post_non*coef_auc)
mean_post_occ_auc = sum(post_occ*coef_auc)
mean_post_daily_auc= sum(post_daily*coef_auc)

## Group SE at post event for table 3 - AUC
se_post_non_auc = sum(post_non*se_auc)
se_post_occ_auc = sum(post_occ*se_auc)
se_post_daily_auc = sum(post_daily*se_auc)

mean_post_non_auc = str_c(round(mean_post_non_auc,2), 
                         " (", round(se_post_non_auc,2), ")")
mean_post_occ_auc = str_c(round(mean_post_occ_auc,2), 
                         " (", round(se_post_occ_auc,2), ")")
mean_post_daily_auc = str_c(round(mean_post_daily_auc,2), 
                         " (", round(se_post_daily_auc,2), ")")

###Contrast statements for differences Differences; Post - pre###
non_dif = (post_non - pre_non)
occ_dif = (post_occ - pre_occ)
daily_dif = (post_daily - pre_daily)
se_non_dif = sum(non_dif*se_min)

non_occ_dif = occ_dif - non_dif 

# Non vs daily

non_daily_dif = daily_dif - non_dif 

## Occ vs Daily
occ_daily_dif =  daily_dif -occ_dif 
```

# Functions for estimates and p-values of difference of differences using contrast statements 

```{r}
### Function for pairwise difference of differences
contrast_pair = function(gee_model, contrast_statment, df, lower_tail = FALSE, log_p = FALSE){
  
  coef2 <- summary(gee_model)$coefficients[,1] # coefficients from gee model
  cov2 <- summary(gee_model)$cov.scaled # scaled covariance from gee model
  contrast_statment.m <- matrix(contrast_statment, byrow = TRUE, 1, 7) # convert contrast statment to 1 by 7 matrix

  contrast_statment_est <- contrast_statment.m%*%coef2 # get estimate of difference from contrast statement * coefficients

  chisq.v <- t(contrast_statment_est)%*%solve(contrast_statment.m%*% cov2  %*%t(contrast_statment.m)) %*% contrast_statment_est # calculate chi-squared test statistic
  
  cont_p_value  <- pchisq(chisq.v, df, ncp=0, lower.tail=lower_tail, log.p=log_p) # calculate p value
  cont_p_value = format.pval(cont_p_value, digits = 2)
  
  contrast_statment_est = round(contrast_statment_est, 3)
  
  return (list(contrast_statment_est, cont_p_value))
}


#### Function for Overall difference of differences
contrast_overal = function(gee_model, dif1, dif2, dif3, df, lower_tail = FALSE, log_p = FALSE){
  
  coef2 <- summary(gee_model)$coefficients[,1] # coefficients from gee model
  cov2 <- summary(gee_model)$cov.scaled # scaled covariance from gee model
  contrast_statment.m <- matrix(c(dif1, dif2, dif3), byrow=T, 3, 7) # convert contrast statment to 1 by 7 matrix

  contrast_statment_est <- contrast_statment.m%*%coef2 # get estimate of difference from contrast statement * coefficients

  chisq.v <- t(contrast_statment_est)%*%solve(contrast_statment.m%*% cov2  %*%t(contrast_statment.m)) %*% contrast_statment_est # calculate chi-squared test statistic
  
  overall_p_value  <- pchisq(chisq.v, df, ncp=0, lower.tail=lower_tail, log.p=log_p) # calculate p value
  overall_p_value = format.pval(overall_p_value, digits = 2)
  return (overall_p_value)
}
```



## Min. Constriction Results

```{r}

## Pairwise estimates for differences and p-values for table 4
non_daily_results = contrast_pair(min_gee, non_daily_dif, 1) # dif_daily - dif_non

non_occ_results = contrast_pair(min_gee, non_occ_dif, 1) #  dif_occ - Dif_non 

occ_daily_results = contrast_pair(min_gee, occ_daily_dif, 1)#  dif_daily - Dif_occ 

overall_min_constrict = contrast_overal(min_gee, non_dif, occ_dif, daily_dif, df = 3 )



### Creating row for table 

row_min_con = t(as.matrix(c(unlist(non_daily_results), unlist(non_occ_results), unlist(occ_daily_results), overall_min_constrict)))
colnames(row_min_con) = c("estimate", "p-value", "estimate", "p-value", "estimate", "p-value", "overall")

```


## AUC results

```{r}

## Pairwise estimates for differences and p-values for table 4
non_daily_results_auc = contrast_pair(auc_gee, non_daily_dif, 1) # dif_daily - dif_non

non_occ_results_auc = contrast_pair(auc_gee, non_occ_dif, 1) #dif_occ - Dif_non 

occ_daily_results_auc = contrast_pair(auc_gee, occ_daily_dif, 1)#  dif_daily - Dif_occ 


overall_auc = contrast_overal(auc_gee,non_dif, occ_dif, daily_dif, df = 3 )



### Creating row for table 

row_auc = t(as.matrix(c(unlist(non_daily_results_auc), unlist(non_occ_results_auc), unlist(occ_daily_results_auc), overall_auc)))
colnames(row_auc) = c("estimate", "p-value", "estimate", "p-value", "estimate", "p-value", "overall")

```



## Table 3 ; Estimated Differences ; Non vs Occ vs Daily 

```{r, eval = FALSE}
#tab_differences = round(rbind(row_min_con, row_auc),4)

tab_differences =rbind(row_min_con, row_auc)

rownames(tab_differences) = c("Percent Constriction", "AUC")

kable(tab_differences,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Table 3") %>%
  kable_classic("striped", full_width = F,
                position = "center", font_size = 12) %>%
  add_header_above(c("Scalar Feature"=1, "Daily vs. Non" = 2, "Occ. vs. Non" = 2, "Daily vs. Occ" =2, " "=1))

```



# Table 4

## Gee model with non vs user (general)

```{r}

#### Taking the easy way out and making a new variable 

scalars = scalars %>% mutate(group_general = case_when(
  user_type == "non-user" ~ "non-user",
  TRUE ~ "user"
))

## Percent constriction
min_gee_general <- geeglm(min_constriction ~ eye + prepost*group_general,
               data = scalars, 
               id = factor(subject_id), 
               family = gaussian,
               corstr = "independence")

## AUC
auc_gee_general <- geeglm(AUC ~ eye + prepost*group_general,
               data = scalars, 
               id = factor(subject_id), 
               family = gaussian,
               corstr = "independence")

```

# Group means 

```{R}

coef_min_gen <- summary(min_gee_general)$coefficients[,1] # coefficients from gee model for min constriction - general user

se_min_gen <- summary(min_gee_general)$coefficients[,2] # SE from gee model for min constriction - general user


coef_auc_gen <- summary(auc_gee_general)$coefficients[,1] # coefficients from gee model for min auc - general user
se_auc_gen <- summary(auc_gee_general)$coefficients[,2] # SE from gee model for min auc - general user


# Int + right + post + general + post*general 

## Pre contrast statements
pre_non_gen = c(1, 0, 0, 0, 0)
pre_use_gen = c(1, 0, 0, 1, 0)

## post contrast statements
post_non_gen = c(1, 0, 1, 0, 0)
post_use_gen = c(1, 0, 1, 1, 1)

general_dif = post_use_gen - pre_use_gen

### Means for general user - min. constriction
mean_pre_use_gen_min = sum(pre_use_gen*coef_min_gen)
mean_post_use_gen_min = sum(post_use_gen*coef_min_gen)

### SE for general user - min. constriction
se_pre_use_gen_min = sum(pre_use_gen*se_min_gen)
se_post_use_gen_min = sum(post_use_gen*se_min_gen)


mean_pre_use_gen_min = str_c(round(mean_pre_use_gen_min, 2), 
                             " (", round(se_pre_use_gen_min,2),
                             ")")
mean_post_use_gen_min = str_c(round(mean_post_use_gen_min, 2), 
                             " (", round(se_post_use_gen_min,2),
                             ")")



## Means for general user - auc
mean_pre_use_gen_auc = sum(pre_use_gen*coef_auc_gen)
mean_post_use_gen_auc = sum(post_use_gen*coef_auc_gen)

## Means for general user - auc
se_pre_use_gen_auc = sum(pre_use_gen*se_auc_gen)
se_post_use_gen_auc = sum(post_use_gen*se_auc_gen)

mean_pre_use_gen_auc = str_c(round(mean_pre_use_gen_auc, 2), 
                             " (", round(se_pre_use_gen_auc,2),
                             ")")
mean_post_use_gen_auc = str_c(round(mean_post_use_gen_auc, 2), 
                             " (", round(se_post_use_gen_auc,2),
                             ")")

# Non vs general

non_general_dif = (post_use_gen - pre_use_gen) - (post_non_gen - pre_non_gen) 

```

# Functions for estimates and p-values of difference of differecences using contrast statements 

```{r}
### Function for pairwise difference of differences
contrast_pair = function(gee_model, contrast_statment, df, lower_tail = FALSE, log_p = FALSE){
  
  coef2 <- summary(gee_model)$coefficients[,1] # coefficients from gee model
  cov2 <- summary(gee_model)$cov.scaled # scaled covariance from gee model
  contrast_statment.m <- matrix(contrast_statment, byrow=T, 1, 5) # convert contrast statment to 1 by 7 matrix

  contrast_statment_est <- contrast_statment.m%*%coef2 # get estimate of difference from contrast statement * coefficients

  chisq.v <- t(contrast_statment_est)%*%solve(contrast_statment.m%*% cov2  %*%t(contrast_statment.m)) %*% contrast_statment_est # calculate chi-squared test statistic
  
  cont_p_value  <- pchisq(chisq.v, df, ncp=0, lower.tail=lower_tail, log.p=log_p) # calculate p value
  cont_p_value = format.pval(cont_p_value, digits = 2)
  
  contrast_statment_est = round(contrast_statment_est, 3)
  
  return (list(contrast_statment_est, cont_p_value))
}

```




## Min. Constriction Results

```{r}

## Pairwise estimates for differences and p-values for table 4
non_general_results = contrast_pair(min_gee_general, non_general_dif, 1) # Dif_non - dif_daily

### Creating row for table 

row_min_con_gen = t(as.matrix(unlist(non_general_results)))
colnames(row_min_con_gen) = c("estimate", "p-value")

```


## AUC results

```{r}
## Pairwise estimates for differences and p-values for table 4
non_general_results_auc = contrast_pair(auc_gee_general, non_general_dif, 1) # Dif_non - dif_daily
 

### Creating row for table 

row_auc_gen = t(as.matrix(unlist(non_general_results_auc)))
colnames(row_auc_gen) = c("estimate", "p-value")

```


# Combined Table 3 + 4 

```{r}
row_auc_gen = t(as.matrix(unlist(non_general_results_auc)))
colnames(row_auc_gen) = c("estimate", "p-value")


row_min_con_gen = t(as.matrix(unlist(non_general_results)))
colnames(row_min_con_gen) = c("estimate", "p-value")

## full row for percent constriction
row_min_con = t(as.matrix(c(unlist(non_daily_results), unlist(non_occ_results), unlist(occ_daily_results), row_min_con_gen)))
colnames(row_min_con) = c("estimate", "p-value", "estimate", "p-value", "estimate", "p-value","estimate", "p-value")

## full row for AUC

row_auc = t(as.matrix(c(unlist(non_daily_results_auc), unlist(non_occ_results_auc), unlist(occ_daily_results_auc), row_auc_gen)))
colnames(row_auc) = c("estimate", "p-value", "estimate", "p-value", "estimate", "p-value",  "estimate", "p-value")


tab_differences = rbind(row_min_con, row_auc)

rownames(tab_differences) = c("Percent Constriction", "AUC")

kable(tab_differences,
      "html",
      booktabs = T,
      align = c("r")) %>%
  kable_classic("striped", full_width = F,
                position = "center", font_size = 12) %>%
  add_header_above(c("Scalar Feature"=1 ,"Daily vs. Non-User"=2 ,"Occasional vs. Non-User"=2, "Daily vs. Occasional"=2 ,"General vs. Non-User"=2 ))
```



# Table 2 - Adjusted means for pre and post 

```{r}
non_user_means = c(mean_pre_non_min, mean_pre_non_auc, mean_post_non_min, mean_post_non_auc)

occ_user_means = c(mean_pre_occ_min, mean_pre_occ_auc, mean_post_occ_min, mean_post_occ_auc)

daily_user_means = c(mean_pre_daily_min, mean_pre_daily_auc, mean_post_daily_min, mean_post_daily_auc)

gen_user_means = c(mean_pre_use_gen_min, mean_pre_use_gen_auc, mean_post_use_gen_min, mean_post_use_gen_auc)


means_table = as.matrix(rbind(non_user_means, occ_user_means, daily_user_means, gen_user_means))


colnames(means_table) = c("Percent Constriction", "AUC", "Percent Constriction", "AUC")

rownames(means_table) = c("non-user", 'occasional', 'daily', 'general user')

kable(means_table,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Table 2") %>%
  kable_classic("striped", full_width = F,
                position = "center", font_size = 12) %>%
  add_header_above(c("User Group"=1,"Baseline"=2, "Post-Smoking" = 2))
```

