new_sc.R

library(sjPlot)

library(lme4)
library(dplyr)
library(nlme)
# library(ggplot)
library(ggplot2)
library(lmerTest)
# library(EMAtools)

lme.dscore <- function(mod, data, type) {
  if (type == "lme4") {
    mod1 <- lmerTest::lmer(mod, data = data)
    eff <- cbind(summary(mod1)$coefficients[, 4], summary(mod1)$coefficients[, 3])
  }

  if (type == "nlme") {
    eff = cbind(summary(mod)$tTable[, 4], summary(mod)$fixDF$terms)
  }

  colnames(eff) <- c("t", "df")
  eff <- as.data.frame(eff)
  eff$d <- (2 * eff$t) / sqrt(eff$df)
  eff <- eff[-1,]
  
  # Format estimates and errors in scientific notation
  eff$t <- formatC(eff$t, format = "e", digits = 2)
  eff$d <- formatC(eff$d, format = "e", digits = 2)
  
  return(eff)
}

####################################### Group #####################################################
df_i <- read.csv("D:/Kuake_Backup/icdl_ff_homo/Results_2223/fformations_each_version/version4+no bad time/4_P_group_sclag.csv")
# df_i <- read.csv("D:/icdl_ff_homo/icdl_ff_homo/Intermediate_results/3_P_group_sclag.csv")
# df_i <- read.csv("D:/icdl_ff_homo/icdl_ff_homo/Intermediate_results/2_P_group_sclag.csv")
# df_i <- read.csv("D:/icdl_ff_homo/icdl_ff_homo/PO_individual_group_sclag_3.csv")
# df_i <- read.csv("D:/Kuake_Backup/icdl_ff_homo/Results_2324/persecond+nobadtime/4_P_group_sclag.csv")

pa <- df_i %>% mutate(
  diagnosisPerson = factor(diagnosisPerson, levels = c("TH", "HL")),
)
head(pa)

# model <- lmer(Social_contact_ratio ~ diagnosisPerson * Group_Size+ (1 | person) + (1|group), data = pa)
# model <- lmer(Social_contact_ratio ~ diagnosisPerson * HL_ratio  + (1 | person) + (1|group), data = pa)
# model <- lmer(Social_contact_ratio ~ diagnosisPerson * Homophily_degree  + (1 | person) + (1|group), data = pa)
# model <- lmer(Social_contact_ratio ~ diagnosisPerson * HL_ratio * Group_Size + (1 | person) + (1|group), data = pa)
# model <- lmer(Social_contact_ratio ~ diagnosisPerson * Homophily_degree * Group_Size + (1 | person) + (1|group), data = pa)
# VarCorr(model)

# To check whether or not to include group as a random factor
model1 <- lmer(Social_contact_ratio ~ diagnosisPerson * Homophily_degree + (1 | person), data = pa) # Without group
model2 <- lmer(Social_contact_ratio ~ diagnosisPerson * Homophily_degree + (1 | group) + (1 | person), data = pa) # With group

anova(model1, model2)
VarCorr(model1)

# # Extract coefficients and format
# coef_summary <- summary(model)$coefficients
# coef_table <- data.frame(
#   Predictors = rownames(coef_summary),
#   Estimates = formatC(coef_summary[, 1], format = "e", digits = 2),
#   Std_Error = formatC(coef_summary[, 2], format = "e", digits = 2),
#   CI = apply(coef_summary[, c(1, 2)], 1, function(x) formatC(x[1] + c(-1, 1) * 1.96 * x[2], format = "e", digits = 2)),
#   t_Value = formatC(coef_summary[, 3], format = "e", digits = 2),
#   p_Value = formatC(coef_summary[, 4], format = "e", digits = 2)
# )

# # Clean up CI column
# coef_table$CI <- sapply(coef_table$CI, function(x) paste0("[", x[1], ", ", x[2], "]"))

# # Print the custom summary table
# print(coef_table)

# Displaying the model using sjPlot
lme.dscore(model1,data = pa,type='lme4')
tab_model(model, show.stat = TRUE, show.se = TRUE, show.r2 = TRUE, dv.labels = "Time in Social Contact")



######################################################################################################
#Dyadic Social Contact Analysis

#read in dataframe
# df<-read.csv("D:/icdl_ff_homo/icdl_ff_homo/PO_individual_group_sclag_3.csv")
# df<-read.csv("D:/icdl_ff_homo/icdl_ff_homo/final_social_contact_ratios_with_diagnosis.csv")
# df_ds_p<-read.csv("D:/icdl_ff_homo/icdl_ff_homo/grouped_final_social_contact_ratios_with_diagnosis_copy.csv")
# df_ds_po<-read.csv("D:/icdl_ff_homo/icdl_ff_homo/PO_grouped_final_social_contact_ratios_with_diagnosis_copy.csv")
# df_k<-read.csv("D:/icdl_ff_homo/icdl_ff_homo/kmeans_final_social_contact_ratios_with_diagnosis_copy.csv")
# df<-read.csv("D:/Kuake_Backup/icdl_ff_homo/Intermediate_results/2_sc.csv")
# df<-read.csv("D:/Kuake_Backup/icdl_ff_homo/Intermediate_results/3_sc.csv")
# df<-read.csv("D:/Kuake_Backup/icdl_ff_homo/Intermediate_results/4_sc.csv")
# df<-read.csv("D:/Kuake_Backup/icdl_ff_homo/Intermediate_results/4_45_sc.csv")
# df<-read.csv("D:/Kuake_Backup/icmi_zff-main/Intermediate_results/4_no_45_sc.csv")
# df<-read.csv("D:/Kuake_Backup/icdl_ff_homo/Results_2223/Intermediate_results/Total Time Updated/4_sc_count.csv")
head(df)

pa <- df %>% mutate(
  diagnosisPerson1 = factor(diagnosisPerson1, levels = c("HL", "TH")),
  diagnosisPerson2 = factor(diagnosisPerson2, levels = c("HL", "TH")),
  Pair = ifelse(Pair == "diff", "Discordant", "Concordant"),
  dyad = paste(Subject, Partner)
)
head(pa)


#aggregated dx*pair social contact analyses

# No grouping

sc <- lmer(scLag ~ diagnosisPerson1 * Pair + (1 | Subject), pa)
lme.dscore(sc,data = df,type='lme4')
tab_model(sc,show.stat = T,show.se=T,show.r2=F,dv.labels ="Shared Time in Social Contact")

# + date as a random factor
sc <- lmer(scLag ~ diagnosisPerson1 * Pair + (1 | Date) + (1 | Subject), pa)

# Add Group
# as fixed factor
# sc <- lmer(scLag ~ diagnosisPerson1 * Homophily * Group + (1 | Subject), pa)

# sc <- lmer(scLag ~ diagnosisPerson1 * Homophily * Group + (1 | Subject) + (1 | Date), pa)

# as random factor
# sc <- lmer(scLag ~ diagnosisPerson1 * Homophily + (1 | Subject) + (1 | Group), data = pa)

# how Homophily and group interact
# sc <- lmer(scLag ~ Homophily * Group + (1 | Subject), data = pa)


# sc <- lmer(scLag ~ Homophily * Group + (1 | Subject)  + (1 | Date) , data = pa)
# sc <- lmer(scLag ~ diagnosisPerson1 * Homophily + (Homophily | Group)  + (1 | Date) , data = pa)



# # with diagnosisPerson1
# # baseline
# model <- lmer(scLag ~ diagnosisPerson1 * Homophily + (1 | Subject), data = pa)

# # ds with p
# model1_ds_p <- lmer(scLag ~ diagnosisPerson1 * Homophily * Group + (1 | Subject), data = pa_ds_p)
# model2_ds_p <- lmer(scLag ~ diagnosisPerson1 * Homophily + (1 | Subject) + (1 | Group), data = pa_ds_p)

# # ds with po
# model1_ds_po <- lmer(scLag ~ diagnosisPerson1 * Homophily * Group + (1 | Subject), data = pa_ds_po)
# model2_ds_po <- lmer(scLag ~ diagnosisPerson1 * Homophily + (1 | Subject) + (1 | Group), data = pa_ds_po)

# # K means
# model1_k<- lmer(scLag ~ diagnosisPerson1 * Homophily * Group + (1 | Subject), data = pa_k)
# model2_k <- lmer(scLag ~ diagnosisPerson1 * Homophily + (1 | Subject) + (1 | Group), data = pa_k)



# # Baseline Model Without diagnosisPerson1
# model_no_diag <- lmer(scLag ~ Homophily + (1 | Subject), data = pa_ds_p)

# # ds with p Models Without diagnosisPerson1
# model1_ds_p_no_diag <- lmer(scLag ~ Homophily * Group + (1 | Subject), data = pa_ds_p)
# model2_ds_p_no_diag <- lmer(scLag ~ Homophily + (1 | Subject) + (1 | Group), data = pa_ds_p)

# # ds with po Models Without diagnosisPerson1
# model1_ds_po_no_diag <- lmer(scLag ~ Homophily * Group + (1 | Subject), data = pa_ds_po)
# model2_ds_po_no_diag <- lmer(scLag ~ Homophily + (1 | Subject) + (1 | Group), data = pa_ds_po)

# # K-means Models Without diagnosisPerson1
# model1_k_no_diag <- lmer(scLag ~ Homophily * Group + (1 | Subject), data = pa_k)
# model2_k_no_diag <- lmer(scLag ~ Homophily + (1 | Subject) + (1 | Group), data = pa_k)


# sc <- model


# summary(model2_ds_p_no_diag)
# summary(model2_ds_po_no_diag)
# summary(model2_k_no_diag)

# AIC(model, model1_ds_p, model2_ds_p, model1_ds_po, model2_ds_po, model1_k, model2_k)
# BIC(model, model1_ds_p, model2_ds_p, model1_ds_po, model2_ds_po, model1_k, model2_k)

# anova(model1, model2)
# anova(model1, model3)

# summary(model)$varcor
# summary(model1_ds_p)$varcor
# summary(model2_ds_p)$varcor
# summary(model1_ds_po)$varcor
# summary(model2_ds_po)$varcor
# summary(model1_k)$varcor
# summary(model2_k)$varcor


# aic_values <- c(
#   AIC(model_no_diag),
#   AIC(model1_ds_p_no_diag),
#   AIC(model2_ds_p_no_diag),
#   AIC(model1_ds_po_no_diag),
#   AIC(model2_ds_po_no_diag),
#   AIC(model1_k_no_diag),
#   AIC(model2_k_no_diag)
# )

# bic_values <- c(
#   BIC(model_no_diag),
#   BIC(model1_ds_p_no_diag),
#   BIC(model2_ds_p_no_diag),
#   BIC(model1_ds_po_no_diag),
#   BIC(model2_ds_po_no_diag),
#   BIC(model1_k_no_diag),
#   BIC(model2_k_no_diag)
# )

# # Create a data frame to summarize the results
# comparison_df <- data.frame(
#   Model = c("Baseline", "DS with p - Model 1", "DS with p - Model 2", 
#             "DS with po - Model 1", "DS with po - Model 2", 
#             "K Means - Model 1", "K Means - Model 2"),
#   AIC = aic_values,
#   BIC = bic_values
# )

# print(comparison_df)


# # Read the CSV file into a dataframe
# df_3 <-read.csv("D:/icdl_ff_homo/icdl_ff_homo/3_group.csv")

# # Rename the columns
# colnames(df_3)[colnames(df_3) == "homophily_degree"] <- "Homophily"
# colnames(df_3)[colnames(df_3) == "sclag"] <- "scLag"
# colnames(df_3)[colnames(df_3) == "group"] <- "Group"

# # head(df_3)
# sc <- lmer(scLag ~ Homophily + (1 | day) + (1 | Group), data = df_3)
lme.dscore(sc,data = df_3,type='lme4')
tab_model(sc,show.stat = T,show.se=T,show.r2=F,dv.labels ="Time in Social Contact")

# # Extract residuals
# residuals_model <- resid(model)

# # Q-Q plot
# ggplot(data = data.frame(residuals_model), aes(sample = residuals_model)) + 
#   stat_qq() + 
#   stat_qq_line() +
#   ggtitle("Q-Q Plot of Residuals for Baseline Model")

# summary(model2_ds_p_no_diag)
# summary(model2_ds_po_no_diag)