ICAR_abbreviations_reliability.R

### ESTIMATE TEST-RETEST RELIABILITY OF ICAR5

library(pacman)
p_load(kirkegaard, weights, psych, magrittr, ggplot2)
#install kirkegaard from https://github.com/Deleetdk/kirkegaard if necessary

# data --------------------------------------------------------------------

#SAPA data
load("data/sapa/sapaICARData18aug2010thru20may2013.rdata")
d_sapa = sapaICARData18aug2010thru20may2013[-c(1:16)] #rename and subset ICAR data
d_sapa_all = sapaICARData18aug2010thru20may2013
rm(sapaICARData18aug2010thru20may2013) #remove the copy


# scales ------------------------------------------------------------------

#item vectors
l_items = list(ICAR5 = c("VR.04", "VR.19", "LN.58", "MR.46", "R3D.04"),
               ICAR16 = c("VR." + c("04", "16", "17", "19"),
                          "LN." + c("07", "33", "34", "58"),
                          "MR." + c("45", "46", "47", "55"),
                          "R3D." + c("03", "04", "06", "08")),
               ICAR60 = colnames(d_sapa))

# cors --------------------------------------------------------------------

d_sapa_cors = wtd.cors(d_sapa) #get cors
diag(d_sapa_cors) = 1
m_sapa_cors = as.matrix(d_sapa_cors) #matrix version

#latent cors
d_sapa_cors_latent_5 = tetrachoric(d_sapa[l_items$ICAR5])
d_sapa_cors_latent_16 = tetrachoric(d_sapa[l_items$ICAR16])
d_sapa_cors_latent_60 = tetrachoric(d_sapa)


# make keys ---------------------------------------------------------- 

#make keys
keys_list = list(
  ICAR5 = l_items$ICAR5,
  ICAR16 = l_items$ICAR16,
  ICAR60 = l_items$ICAR60)
sapa_keys = make.keys(d_sapa_cors, keys_list)


# get scores --------------------------------------------------------------

# score the items using the keys
d_ICAR_scored = scoreOverlap(sapa_keys, m_sapa_cors, impute="none")

# see the correlations between the items and loadings based scales
print(d_ICAR_scored, short=FALSE)


# IRT ---------------------------------------------------------------------

#this takes time
fa_5 = irt.fa(d_sapa_cors_latent_5)
fa_16 = irt.fa(d_sapa_cors_latent_16)
fa_60 = irt.fa(d_sapa_cors_latent_60)

#get scores
d_scores_5 = score.irt(fa_5, items = d_sapa[l_items$ICAR5])
d_scores_16 = score.irt(fa_16, items = d_sapa[l_items$ICAR16])
d_scores_60 = score.irt(fa_60, items = d_sapa)

#save scores
d_ICAR_scores_IRT = data.frame(ICAR5 = d_scores_5$theta1 %>% standardize(),
                               ICAR16 = d_scores_16$theta1 %>% standardize(),
                               ICAR60 = d_scores_60$theta1 %>% standardize())

#cors
wtd.cors(d_ICAR_scores_IRT)

#cors with criteria vars
d_ICAR_scores_IRT$Edu = d_sapa_all$p2edu %>% as.numeric()

#reliability from criteria cors
wtd.cors(d_ICAR_scores_IRT)[4, 1] / wtd.cors(d_ICAR_scores_IRT)[4, 3]

#BMI?
d_ICAR_scores_IRT$BMI = d_sapa_all$BMI
ggplot(d_ICAR_scores_IRT, aes(ICAR60, BMI)) + geom_smooth()