# Description: Frequency vs. Intensity study 
#              (1) Data: OIDP
#              (2) Method: correlation using factor analysis model 
# Written by Seungwon Chung 

# install packages
library(tidyverse)
library(psych)
library(psychometric)
library(lavaan)

# set directories
folder <- ".../Impact frequency and intensity/"
setwd(folder) 
code.folder <- paste0(folder,"Code/")
data.folder <- paste0(folder,"Data/")
result.folder <- paste0(folder,"Result/")

# read in data (csv file)
OIDP_clean <- read.csv(paste0(data.folder,"OIDP data_clean.csv"))

############################
## Descriptive Statistics ##
############################

# subset item responses (freq/int) 
frequency <- OIDP_clean %>% 
  select(ends_with("freq_i"))
intensity <- OIDP_clean %>% 
  select(ends_with("diff_i"))
fre_int <- cbind(frequency, intensity)

# get sum scores
freq_sum <- apply(frequency,1,sum)
int_sum <- apply(intensity,1,sum)
fre_int_sum <- cbind(freq_sum, int_sum)

# subset group variables
OIDP_clean$group.age <- ifelse(OIDP_clean$age<=median(OIDP_clean$age),0,1) # define age group (0: young, 1: old)
OIDP_scores_group <- cbind(fre_int_sum,OIDP_clean$language,OIDP_clean$group.age,OIDP_clean$female) 
OIDP_scores_group <- as.data.frame(OIDP_scores_group)

names(OIDP_scores_group) <- c("fre_sum","int_sum","language","age","gender")

# summary (frequency)
OIDP_scores_group %>%                               
  summarize(min = min(fre_sum,na.rm=TRUE),
            q1 = quantile(fre_sum, 0.25,na.rm=TRUE),
            median = median(fre_sum,na.rm=TRUE),
            mean = mean(fre_sum,na.rm=TRUE),
            q3 = quantile(fre_sum, 0.75,na.rm=TRUE),
            max = max(fre_sum,na.rm=TRUE))

# summary (severity)
OIDP_scores_group %>%                               
  summarize(min = min(int_sum,na.rm=TRUE),
            q1 = quantile(int_sum, 0.25,na.rm=TRUE),
            median = median(int_sum,na.rm=TRUE),
            mean = mean(int_sum,na.rm=TRUE),
            q3 = quantile(int_sum, 0.75,na.rm=TRUE),
            max = max(int_sum,na.rm=TRUE))

# summary by language (frequency)
OIDP_scores_group %>%                               
  group_by(language) %>% 
  summarize(min = min(fre_sum,na.rm=TRUE),
            q1 = quantile(fre_sum, 0.25,na.rm=TRUE),
            median = median(fre_sum,na.rm=TRUE),
            mean = mean(fre_sum,na.rm=TRUE),
            q3 = quantile(fre_sum, 0.75,na.rm=TRUE),
            max = max(fre_sum,na.rm=TRUE))

# summary by language (severity)
OIDP_scores_group %>%                               
  group_by(language) %>% 
  summarize(min = min(int_sum,na.rm=TRUE),
            q1 = quantile(int_sum, 0.25,na.rm=TRUE),
            median = median(int_sum,na.rm=TRUE),
            mean = mean(int_sum,na.rm=TRUE),
            q3 = quantile(int_sum, 0.75,na.rm=TRUE),
            max = max(int_sum,na.rm=TRUE))

# summary by age (frequency)
OIDP_scores_group %>%                               
  group_by(age) %>% 
  summarize(min = min(fre_sum,na.rm=TRUE),
            q1 = quantile(fre_sum, 0.25,na.rm=TRUE),
            median = median(fre_sum,na.rm=TRUE),
            mean = mean(fre_sum,na.rm=TRUE),
            q3 = quantile(fre_sum, 0.75,na.rm=TRUE),
            max = max(fre_sum,na.rm=TRUE))

# summary by age (severity)
OIDP_scores_group %>%                               
  group_by(age) %>% 
  summarize(min = min(int_sum,na.rm=TRUE),
            q1 = quantile(int_sum, 0.25,na.rm=TRUE),
            median = median(int_sum,na.rm=TRUE),
            mean = mean(int_sum,na.rm=TRUE),
            q3 = quantile(int_sum, 0.75,na.rm=TRUE),
            max = max(int_sum,na.rm=TRUE))

# summary by gender (frequency)
OIDP_scores_group %>%                               
  group_by(gender) %>% 
  summarize(min = min(fre_sum,na.rm=TRUE),
            q1 = quantile(fre_sum, 0.25,na.rm=TRUE),
            median = median(fre_sum,na.rm=TRUE),
            mean = mean(fre_sum,na.rm=TRUE),
            q3 = quantile(fre_sum, 0.75,na.rm=TRUE),
            max = max(fre_sum,na.rm=TRUE))

# summary by gender (severity)
OIDP_scores_group %>%                               
  group_by(gender) %>% 
  summarize(min = min(int_sum,na.rm=TRUE),
            q1 = quantile(int_sum, 0.25,na.rm=TRUE),
            median = median(int_sum,na.rm=TRUE),
            mean = mean(int_sum,na.rm=TRUE),
            q3 = quantile(int_sum, 0.75,na.rm=TRUE),
            max = max(int_sum,na.rm=TRUE))

#############################
## Correlation among items ##
#############################

# polychoric correlations
Pcorr_fre_int <- lavCor(fre_int, ordered=colnames(c(frequency,intensity)),se="standard",output="est")

################################
## Correlation b/w sum scores ##
################################

freq_sum <- apply(frequency,1,sum)
int_sum <- apply(intensity,1,sum)
cor_sum <- cor(freq_sum,int_sum,use="pairwise.complete.obs")
CIr(cor_sum, n=873, level=.95) 

#####################################
## Correlation using latent scores ##
#####################################

# correlation without correlated residuals 
two_factor <- 'Frequency =~ oidp1_freq_i + oidp2_freq_i + oidp3_freq_i + oidp4_freq_i + oidp5_freq_i + oidp6_freq_i + oidp7_freq_i + oidp8_freq_i
               Intensity =~ oidp1_diff_i + oidp2_diff_i + oidp3_diff_i + oidp4_diff_i + oidp5_diff_i + oidp6_diff_i + oidp7_diff_i + oidp8_diff_i
               Frequency~~ 1*Frequency
               Intensity~~ 1*Intensity
               Frequency~~ Intensity'

#fit.model <- cfa(two_factor,data = fre_int, std.lv=TRUE, ordered=colnames(c(frequency,intensity)),estimator="ULS")
fit.model_two_factor <- cfa(two_factor,data = fre_int, std.lv=TRUE, ordered=colnames(c(frequency,intensity)),estimator="DWLS")
summary(fit.model_two_factor, fit.measures = TRUE, standardized=TRUE)
parameterEstimates(fit.model_two_factor,se=TRUE)

# correlation with correlated residuals 
two_factor_ecor <- 'Frequency =~ oidp1_freq_i + oidp2_freq_i + oidp3_freq_i + oidp4_freq_i + oidp5_freq_i + oidp6_freq_i + oidp7_freq_i + oidp8_freq_i
               	Intensity =~ oidp1_diff_i + oidp2_diff_i + oidp3_diff_i + oidp4_diff_i + oidp5_diff_i + oidp6_diff_i + oidp7_diff_i + oidp8_diff_i
               	Frequency~~ 1*Frequency
               	Intensity~~ 1*Intensity
               	Frequency~~ Intensity
		   	oidp1_freq_i~~oidp1_diff_i
			oidp2_freq_i~~oidp2_diff_i
			oidp3_freq_i~~oidp3_diff_i
			oidp4_freq_i~~oidp4_diff_i
			oidp5_freq_i~~oidp5_diff_i
			oidp6_freq_i~~oidp6_diff_i
			oidp7_freq_i~~oidp7_diff_i
			oidp8_freq_i~~oidp8_diff_i'
fit.model_two_factor_ecor <- cfa(two_factor_ecor,data = fre_int, std.lv=TRUE, ordered=colnames(c(frequency,intensity)),estimator="DWLS")
summary(fit.model_two_factor_ecor, fit.measures = TRUE, standardized=TRUE)
parameterEstimates(fit.model_two_factor_ecor,se=TRUE)