## Vigotsky et al., Methods Matter, PeerJ, 2018
## By Andrew Vigotsky, avigotsky@gmail.com
## Last edited: April 10, 2018

library(here)
setwd(here())
library(lme4)
library(dplyr)
library(foreign)
library(lmerTest)
library(sjstats)
library(car)
library(psychometric)
data0 <- read.csv('strength_hypertrophy_dataset.csv')


# Measurement correlation ##############################
# Calculates correlation matrix of different measurements 
measurements <- c("us_30","us_50", "us_70","us_avg","csa")
round(cor(data0[,measurements],use="pairwise.complete.obs"),3) # between subject
# within subject, weighted by # of points
output <- list()
for(i in 1:length(measurements)) {
  measure_cor <- c()
  for(j in 1:length(measurements)) {
    points <- c()
    sub_cor <- c()
    for(k in 1:19) {
      df.current <- data0[data0$subject == k & !is.na(data0[,measurements[j]]) & !is.na(data0[,measurements[i]]),]
      kpts <- nrow(df.current)
      if(kpts > 2) { # avoid correlations w/ <= 2pts
        points[k] <- kpts
        sub_cor[k] <- r2z(cor(df.current[,measurements[i]],df.current[,measurements[j]],use="pairwise.complete.obs")) # converts to z
      }
    }
    measure_cor[j] <- sum(sub_cor*(points-3))/sum(points-3) # weight by df
  }
  output[[i]] <- z2r(measure_cor) # back to r
}
df <- do.call(rbind,output)
rownames(df) <- measurements
colnames(df) <- measurements
print(round(df,3))


# Between-subject analyses ##############################
# These analyses are the traditional "between-subject" analyses. They correlate each all subjects' (Δsize,Δstrength) with one another. This is split into two periods - training and detraining.
# Note that all correlations are squared to obtain R^2 estimates (VAF).
phase <- c("training","detraining")
output <- list()
for(k in 1:length(phase)) { # let 1=training; 2=detraining
  pre <- data0[data0$week == (k-1)*8,]
  mid <- data0[data0$week == k*4,] # needed for CSA
  post <- data0[data0$week == k*8,]
  change <- post - pre
  mid_change <- post - mid
  
  # BOOTSTRAP 95% CIs
  us_30 <- c()
  us_50 <- c()
  us_70 <- c()
  us_avg <- c()
  csa <- c()
  set.seed(368271)
  for(i in 1:2000) {
    subsamp <- sample(1:19,19,replace=T)
    temp <- data.frame()
    mid_temp <- data.frame()
    for(j in 1:19) {
      temp <- rbind(temp,change[subsamp[j],])
      mid_temp <- rbind(mid_temp,mid_change[subsamp[j],])
    }
    us_30[i] <- (cor(temp$mvc,temp$us_30,use="pairwise.complete.obs")^2*100)
    us_50[i] <- (cor(temp$mvc,temp$us_50,use="pairwise.complete.obs")^2*100)
    us_70[i] <- (cor(temp$mvc,temp$us_70,use="pairwise.complete.obs")^2*100)
    us_avg[i] <- (cor(temp$mvc,temp$us_avg,use="pairwise.complete.obs")^2*100)
    csa[i] <- (cor(mid_temp$mvc,mid_temp$csa,use="pairwise.complete.obs")^2*100)
    if(i %% 100 == 0)
      print(i)
  }
  df = data.frame("phase" = phase[k],
                  "us_30" = round(cor(change$mvc,change$us_30,use="pairwise.complete.obs")^2*100,1),
                  "us_30_lower" = round(quantile(us_30,prob=0.025),1),
                  "us_30_upper" = round(quantile(us_30,prob=0.975),1),
                  "us_50" = round(cor(change$mvc,change$us_50,use="pairwise.complete.obs")^2*100,1),
                  "us_50_lower" = round(quantile(us_50,prob=0.025),1),
                  "us_50_upper" = round(quantile(us_50,prob=0.975),1),
                  "us_70" = round(cor(change$mvc,change$us_70,use="pairwise.complete.obs")^2*100,1),
                  "us_70_lower" = round(quantile(us_70,prob=0.025),1),
                  "us_70_upper" = round(quantile(us_70,prob=0.975),1),
                  "us_avg" = round(cor(change$mvc,change$us_avg,use="pairwise.complete.obs")^2*100,1),
                  "us_avg_lower" = round(quantile(us_avg,prob=0.025),1),
                  "us_avg_upper" = round(quantile(us_avg,prob=0.975),1),
                  "csa" = round(cor(mid_change$mvc,mid_change$csa,use="pairwise.complete.obs")^2*100,1),
                  "csa_lower" = round(quantile(csa,prob=0.025),1),
                  "csa_upper" = round(quantile(csa,prob=0.975),1)
                )
  output[[k]] <- df
}

print(do.call(rbind,output))


# Within-subject analyses  ##############################
# The following code looks at the within-subject relationships between strength and hypertrophy.
# The bootstrap code will ocassionally throw convergence warnings. This does not seem to be too problematic, however, as different algorithms appear to produce similar estimates and we are not looking at *very* extremes (min and max).

measurements <- c("csa","us_avg","us_70","us_50","us_30")
output <- list()
for(j in 1:length(measurements)) {
  data <- data0[complete.cases(data0[,measurements[j]]),]
  grpmeans<-aggregate(data[,measurements[j]],list(data$subject),mean,na.rm=TRUE)
  names(grpmeans)<-c("subject","grpmn")
  data<-merge(data,grpmeans,by="subject")
  data$centered<-data[,measurements[j]]-data$grpmn
  full <- lmer(mvc ~ centered + (centered|subject), data=data)
  fixed <- lmer(mvc ~ 1 + (1|subject), data=data)
  # HLM R^2
  hlmr2 <- 1-var(resid(full))/var(resid(fixed))
  hlmicc <- icc(fixed)
  # ANCOVA
  form <- formula(paste("mvc ~ ",measurements[j],"+factor(subject)"))
  peta2 <- as.vector(eta_sq(Anova(lm(form,data),type=3),partial=T))[2]
  # BOOTSTRAP 95% CIs
  r2 <- c()
  etasquared <- c()
  set.seed(368271)
  for(i in 1:2000) {
    subsamp <- sample(1:19,19,replace=T)
    temp <- data.frame()
    for(k in 1:19) {
      temp <- rbind(temp,data[data$subject %in% subsamp[k],])
    }
    full <- lmer(mvc ~ centered + (centered|subject), data=temp)
    fixed <- lmer(mvc ~ 1 + (1|subject), data=temp)
    r2[i] <- 1-var(resid(full))/var(resid(fixed))
    etasquared[i] <- as.vector(eta_sq(Anova(lm(form,temp),type=3),partial=T))[2]
    if(i %% 100 == 0)
      print(i)
  }
  df <- data.frame("measurement" = measurements[j],
             "HLM_VAF" = round(hlmr2*100,1),
             "HLM_LOWER" = round(quantile(r2,prob=0.025)*100,1),
             "HLM_UPPER" = round(quantile(r2,prob=0.975)*100,1),
             "ANCOVA_VAF" = round(peta2*100,1),
             "ANCOVA_LOWER" = round(quantile(etasquared,prob=0.025)*100,1),
             "ANCOVA_UPPER" = round(quantile(etasquared,prob=0.975)*100,1),
             "ICC" = round(hlmicc,2)
             )
  output[[j]] <- df
}

print(do.call(rbind,output))


# Relationship between the strength-hypertrophy relationship and strength ##############
# This tests to see if those who have a stronger strength-hypertrophy correlation also tend to get stronger.
del_mvc <- c()
subject_csa <- c()
subject_us_avg <- c()
subject_us_30 <- c()
subject_us_50 <- c()
subject_us_70 <- c()
for(i in 1:19) {
  del_mvc[i] <- data0[data0$week == 8 & data0$subject == i,"mvc"] - data0[data0$week == 0 & data0$subject == i,"mvc"]
  if(length(na.omit(data0[data0$subject == i,"csa"])) > 2) # This stops cases with only 2 measures from inflating estimates
    subject_csa[i] <- cor(data0[data0$subject == i,"mvc"],data0[data0$subject == i,"csa"],use="pairwise.complete.obs")
  else
    subject_csa[i] <- NA
  subject_us_avg[i] <- cor(data0[data0$subject == i,"mvc"],data0[data0$subject == i,"us_avg"],use="pairwise.complete.obs")
  subject_us_30[i] <- cor(data0[data0$subject == i,"mvc"],data0[data0$subject == i,"us_30"],use="pairwise.complete.obs")
  subject_us_50[i] <- cor(data0[data0$subject == i,"mvc"],data0[data0$subject == i,"us_50"],use="pairwise.complete.obs")
  subject_us_70[i] <- cor(data0[data0$subject == i,"mvc"],data0[data0$subject == i,"us_70"],use="pairwise.complete.obs")
}
df.change <- data.frame(del_mvc,subject_csa,subject_us_avg,subject_us_30,subject_us_50,subject_us_70)
for(i in 1:5) {
  print(colnames(df.change)[i+1])
  c <- cor(df.change[i+1],df.change[,1],method="spearman",use="pairwise.complete.obs")
  print(c[[1]])
}