Better_RIR_2.R

####WORKSPACE SETTING####

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
dati <- read.csv("data_careless.csv")

library(dplyr)

## 
## Caricamento pacchetto: 'dplyr'

## I seguenti oggetti sono mascherati da 'package:stats':
## 
##     filter, lag

## I seguenti oggetti sono mascherati da 'package:base':
## 
##     intersect, setdiff, setequal, union

####BETTER RIR FUNCTION####
rir <- function(data,factors,iterations=100){
  

#Calculating RIR    
  all_IR <- matrix(nrow=nrow(data),ncol=iterations) %>%  data.frame()
  d <- 1 #counter
  
  for (i in 1:iterations){
    
    #necessary in subsequent for cycle
    a <- 1
    b <- 0
    indexes_half1 <- list()
    indexes_half2 <- list()
    
    #Splitting the dataset in two, making sure that items of the same subscale have the same index
    for (length_factor in factors){
      
      #Sample randomly one half of the subscale, assignign those indeces to 
      #halfSubscale_1 and the others to half_subscale 2
      halfSubscale_1 <- sample(1:length_factor,floor(length_factor/2))
      halfSubscale_2 <- (1:length_factor)[!1:length_factor%in%halfSubscale_1]
      
      #If they are of different length, remove a random item from the longer one
      if (length(halfSubscale_2) != length(halfSubscale_1)){
        halfSubscale_2 <- halfSubscale_2[sample(halfSubscale_2,1)!=halfSubscale_2]
      }
      
      #Storing the indexes in two lists of all the indexes in the two halves
      indexes_half1[[a]] <- halfSubscale_1 + b
      indexes_half2[[a]] <- halfSubscale_2 + b
      
      #Counter of list position
      a <- a+1
      
      #Counter of total length of indeces
        b <- b + length_factor
    }
    
    #making the list into a vector
    indexes_half1 <- indexes_half1 %>% unlist
    indexes_half2 <- indexes_half2 %>% unlist
    
    #calculating IR for each participant
    IR <- apply(data,
          1,
          \(x){
            cor(x[indexes_half1] %>% as.numeric,
                x[indexes_half2] %>% as.numeric,
                use="pairwise.complete.obs") %>%
              suppressWarnings()
          })  
    
    #storing it into a dataframe (participant x iteration)
    all_IR[,d] <- IR
    d <- d+1
    
    #averaging the index of each participant
    RIR <- apply(all_IR,1,mean,na.rm=T)
    
  }

  return(list(RIR=RIR,all_IR=all_IR))
}





#### RRIC - Resampled Random Individual Correlation ####

rric <- function(data,iterations=100){

  #small function to transform the numbers to the smallest even number
  make_even <- function(x) {
    ifelse(x %% 2 == 0, x, x - 1)
  }
  
  #needed for for cycle below
  e <- 1
  all_RIC <- matrix(nrow=nrow(data),ncol=iterations)
  
  for (i in 1:iterations){
  
  #Giving the columns a random order
  col_random_order <- sample(1:ncol(data),make_even(ncol(data)))
  
  #splitting the columns in two halves
  indexes_half1 <- col_random_order[1:(length(col_random_order)/2)]
  indexes_half2 <- col_random_order[(1+length(col_random_order)/2):length(col_random_order)]
  
  #computing the random individual correlation between the two halves
  RIC <- apply(data,
              1,
              \(x){
                cor(x[indexes_half1] %>% as.numeric,
                    x[indexes_half2] %>% as.numeric,
                    use="pairwise.complete.obs") %>%
                  suppressWarnings()
              })  
  
  #storing results inside a data.frame case x iteration
  all_RIC[,e] <- RIC
  e <- e+1
  }
  
  #averaging all the RIC, making it RRIC
  RRIC <- apply(all_RIC,1,mean,na.rm=T)
  
  return(list(all_RIC = all_RIC, RRIC = RRIC))
}


#### RESAMPLED INDIVIDUAL CORRELATION DIFFERECE ####

ricd <- function(data,factors,iter=20){

  rir_results <- rir(dati,c(10,10,10),iter)
  
  rric_results <- rric(dati,iter)

  RICD <- abs(rir_results$RIR - rric_results$RRIC)
  
  return(RICD)
}


ricd_results <- ricd(
  data = dati,
  factors = c(rep(3,5),rep(5,5),rep(3,6),rep(7,3),rep(6,5)),
  iter = 100
)

plot(density(ricd_results%>%na.omit))

ricd_results <- ricd(
  data = dati[,sample(ncol(data))],
  factors = c(rep(3,5),rep(5,5),rep(3,6),rep(7,3),rep(6,5)),
  iter = 100
)

(rric(dati,100)$RRIC-rric(dati,100)$RRIC) %>% 
  abs %>% 
  density %>%
  plot