variance decomposition for difference scores - inspired by Hedge 2018

diffScoreVarDecomp

diffScoreVarDecomp <- function(cA_m1, cB_m1, cA_m2, cB_m2, plot = "both", measurename = "x") {

# calculations are heavily based on the irr's package icc() code:

# bind the four vectors first for row-wise deletion:

alldat <- as.matrix(na.omit(cbind(cA_m1, cB_m1, cA_m2, cB_m2))) 

alldat <- cbind(alldat, "ds1" = alldat[,"cA_m1"] - alldat[,"cB_m1"], "ds2" = alldat[,"cA_m2"] - alldat[,"cB_m2"])

datvector <- c("componentscore A", "componentscore B", "differencescore A-B")  

    ns <- nrow(alldat)
    nr <- 2
    
# ns = number of subjects; 
# nr = number of raters/measurements - always two in this set-up
 

decompdat <- NULL

for (dat in 1:length(datvector)) {

  if(dat ==1){  usedat <- alldat[,c("cA_m1", "cA_m2")] }
  if(dat ==2){  usedat <- alldat[,c("cB_m1", "cB_m2")] }
  if(dat ==3){  usedat <- alldat[,c("ds1", "ds2")]}
  
    SStotal <- var(as.numeric(usedat))*(ns*nr-1)
    MSr <- var(apply(usedat,1,mean))*nr                    # MSR = mean square for rows - between subjects
    MSw <- sum(apply(usedat,1,var)/ns)                     # MSW = mean square for residual sources of variance - SPSS error between table 
    MSc <- var(apply(usedat,2,mean))*ns                    # MSC = mean square for columns - between sessions 
    MSe <- (SStotal-MSr*(ns-1)-MSc*(nr-1))/((ns-1)*(nr-1))  # MSE = mean square for error - SPSS error within table 

icc_A <- (MSr-MSe)/(MSr+(nr-1)*MSe+(nr/ns)*(MSc-MSe))
icc_C <- (MSr-MSe)/(MSr+(nr-1)*MSe)

## then for the first decomposing step, Hedge reduces the subject and session MS by the error MS (as is done in the ICC formulas) and divides again by the df: 

Betw_subject_var_abs <- (MSr - MSe)/nr
Betw_session_var_abs <- (MSc - MSe)/ns
Error_var_abs <- MSe

## and for relative variance: treat the sum of between subects, between sessions, and between error var as var total and compute proportion for each:
totalVar <- sum(Betw_subject_var_abs, Betw_session_var_abs, Error_var_abs)
Betw_subject_var_prc <- Betw_subject_var_abs / totalVar
Betw_session_var_prc <- Betw_session_var_abs / totalVar
Error_var_prc <- Error_var_abs / totalVar

mdat=data.frame(measure = measurename, score= datvector[dat], icc_Agreement = icc_A, icc_Consistency = icc_C, var_source = c("between subjects", "between sessions", "error"), var_abs = c(Betw_subject_var_abs, Betw_session_var_abs, Error_var_abs), var_prc = c(Betw_subject_var_prc, Betw_session_var_prc, Error_var_prc))

decompdat <- rbind(decompdat, mdat)
}

decompdat <- as.data.frame(decompdat)

# plotterdeplot: 
if(plot != "none") {
require(ggplot2)

prcplot <- ggplot(decompdat, aes(x = score, y=var_prc, fill=factor(var_source, levels=rev(c("between subjects", "between sessions", "error") )))) + 
    geom_bar( stat="identity") + 
    theme_minimal() + 
    theme(legend.title = element_blank(), panel.grid = element_blank(), axis.title.x = element_blank()) + 
    labs(title= decompdat$measure[1], y= "variance percent\n") + 
    scale_fill_manual(values=c("#cd1076", 	"#76cd10", 	"#1076cd"))

absplot <- ggplot(decompdat, aes(x = score, y=var_abs, fill=factor(var_source, levels=rev(c("between subjects", "between sessions", "error") )))) + 
    geom_bar( stat="identity") + 
    theme_minimal() + 
    theme(legend.title = element_blank(), panel.grid = element_blank(), axis.title.x = element_blank()) + 
    labs(title= decompdat$measure[1], y = "variance absolute\n") + 
    scale_fill_manual(values=c("#cd1076", 	"#76cd10", 	"#1076cd"))

if(plot == "both") {
  print(absplot)
  print(prcplot)
  }

if(plot == "absolute") {
  print(absplot)
}


if(plot == "percentage") {
  print(prcplot)
}
}

return(decompdat)
}