Load the neccessary packages

library(plyr)
library(stringr)

Set the working directory and import the data files And naming the columns when group is equal to 1 = OCPD, 3 = HC and 4 = OCD

setwd("C:/Users/World/Desktop/Joe/Research/CDS - Columbia/ITC/DD")
ITC = ldply(list.files(pattern = "[0-9]+__[0-9]+_[0-9]+_[0-9]+_[0-9]+_[0-9]+_[0-9]+_[A-Z]+_TMS_Outcomes.txt"), 
    function(filename) {
        dum = read.table(filename)
        dum$filename = filename
        dum$ID = str_extract(dum$filename, "^[0-9]{4}")
        dum$group = str_extract(dum$filename, "^[0-9]")
        dum <- rename(dum, c(V1 = "TrialN", V2 = "SSd", V3 = "LLd", V4 = "SSa", 
            V5 = "LLa", V6 = "Choice"))
        return(dum)
    })

ITC$check <- ifelse(str_extract(ITC$filename, "^[0-9]{1}") %in% 1 | str_extract(ITC$filename, 
    "^[0-9]{1}") %in% 3:4, "1", "2")

ITC <- ITC[ITC[, 10] != 2, ]

head(ITC)
##   TrialN SSd LLd  SSa  LLa Choice
## 1      1   0   2 44.0 57.2      2
## 2      2   0   4 24.7 27.2      1
## 3      3   2   6 43.9 44.3      1
## 4      4   0   2 18.5 22.2      2
## 5      5   2   6 15.0 17.3      2
## 6      6   0   2 46.1 50.7      2
##                                     filename   ID group check
## 1 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1
## 2 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1
## 3 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1
## 4 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1
## 5 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1
## 6 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1

To check the Ns in each group divide the following by 58 and then by 2 (or 116)

N <- table(ITC$group)
N/116
## 
##  1  3  4 
## 22 23 22

Seperating between Accelerate and Delay. If the file is completed first it is accelerate if later it is delay. Make time column and data column

ITC$time <- str_replace(ITC$filename, "[0-9]+__[0-9]+_[0-9]+_[0-9]+_[0-9]+_([0-9]+)_([0-9]+)_[A-Z]+_TMS_Outcomes.txt", 
    "\\1\\2")
ITC$hour <- as.numeric(str_replace(ITC$filename, "[0-9]+__[0-9]+_[0-9]+_[0-9]+_([0-9]+)_[0-9]+_[0-9]+_[A-Z]+_TMS_Outcomes.txt", 
    "\\1"))
ITC$hour <- ITC$hour + 12
ITC$time <- paste(ITC$hour, ITC$time, sep = "")
ITC$hour <- NULL

ITC$date <- str_replace(ITC$filename, "[0-9]+__([0-9]+)_([0-9]+)_([0-9]+)_[0-9]+_[0-9]+_[0-9]+_[A-Z]+_TMS_Outcomes.txt", 
    "\\1\\2\\3")

If the date and id are equal and the time is less than the next i+58th question it is accelerate(1), if not it is delay(2). Column 11 is equal to ITC$time. Column 8 is equal to ITC$ID and column 12 is equal to ITC$date. Renaming the column

for (i in 1:nrow(ITC)) {
    a = ITC[i, 8]
    b = ITC[i + 58, 8]
    x = ITC[i, 11]
    y = ITC[i + 58, 11]
    w = ITC[i, 12]
    z = ITC[i + 58, 12]
    ITC[i, 13] <- ifelse(a == b & x < y & w == z, 1, 2)
}

colnames(ITC)[13] <- "Acc.D"

Now that we have ID, time, date, group and accelerate vs delay, we can start the analysis Setting the delays in years Starting with the titrator

ITC$SSd <- ifelse(ITC[, 1] %in% 1:36, ITC$SSd/52, ITC$SSd/12)
ITC$LLd <- ifelse(ITC[, 1] %in% 1:36, ITC$LLd/52, ITC$LLd/12)

head(ITC)
##   TrialN     SSd     LLd  SSa  LLa Choice
## 1      1 0.00000 0.03846 44.0 57.2      2
## 2      2 0.00000 0.07692 24.7 27.2      1
## 3      3 0.03846 0.11538 43.9 44.3      1
## 4      4 0.00000 0.03846 18.5 22.2      2
## 5      5 0.03846 0.11538 15.0 17.3      2
## 6      6 0.00000 0.03846 46.1 50.7      2
##                                     filename   ID group check   time
## 1 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 2 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 3 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 4 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 5 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 6 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
##     date Acc.D
## 1 182013     1
## 2 182013     1
## 3 182013     1
## 4 182013     1
## 5 182013     1
## 6 182013     1

dummy.T <- subset(ITC, !(TrialN %in% 1:36))

for (i in 1:nrow(dummy.T)) {
    x = dummy.T[i, 6]
    y = dummy.T[i + 1, 6]
    dummy.T[i, 14] <- ifelse(x < y, (dummy.T[i, 5] + dummy.T[i + 1, 5])/2, "NA")
}

colnames(dummy.T)[14] <- "IP"

head(dummy.T)
##    TrialN SSd  LLd SSa LLa Choice
## 37     37   0 0.25  50  55      2
## 38     38   0 0.25  50  60      2
## 39     39   0 0.25  50  65      2
## 40     40   0 0.25  50  70      2
## 41     41   0 0.25  50  75      2
## 42     42   0 0.25  50  80      2
##                                      filename   ID group check   time
## 37 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 38 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 39 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 40 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 41 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 42 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
##      date Acc.D IP
## 37 182013     1 NA
## 38 182013     1 NA
## 39 182013     1 NA
## 40 182013     1 NA
## 41 182013     1 NA
## 42 182013     1 NA

Calculating the k and delta

dummy.T[, 14] <- as.numeric(dummy.T[, 14])
## Warning: NAs introduced by coercion

for (i in 1:nrow(dummy.T)) {
    x = dummy.T[i, 6]
    y = dummy.T[i + 1, 6]
    dummy.T[i, 15] <- ifelse(x < y, (dummy.T[i, 4]/dummy.T[i, 14])^(dummy.T[i, 
        3] - dummy.T[i, 2]), "NA")
    dummy.T[i, 16] <- ifelse(x < y, ((dummy.T[i, 14]/dummy.T[i, 4]) - 1)^4, 
        "NA")
}

dummy.T <- rename(dummy.T, c(V15 = "T.d", V16 = "T.k"))
head(dummy.T)
##    TrialN SSd  LLd SSa LLa Choice
## 37     37   0 0.25  50  55      2
## 38     38   0 0.25  50  60      2
## 39     39   0 0.25  50  65      2
## 40     40   0 0.25  50  70      2
## 41     41   0 0.25  50  75      2
## 42     42   0 0.25  50  80      2
##                                      filename   ID group check   time
## 37 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 38 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 39 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 40 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 41 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
## 42 1002__1_8_2013_1_22_15_PM_TMS_Outcomes.txt 1002     1     1 132215
##      date Acc.D IP T.d T.k
## 37 182013     1 NA  NA  NA
## 38 182013     1 NA  NA  NA
## 39 182013     1 NA  NA  NA
## 40 182013     1 NA  NA  NA
## 41 182013     1 NA  NA  NA
## 42 182013     1 NA  NA  NA

Calculating group differences

dummy.T <- na.omit(dummy.T)
summary(aov(T.k ~ group, dummy.T))
##              Df Sum Sq Mean Sq F value Pr(>F)
## group         2      9    4.38    2.14   0.12
## Residuals   187    382    2.04

The non titrator Calculating the implied k Recode 0 as impatient and 1 as patient

dummy.NT <- subset(ITC, TrialN %in% 1:36)
dummy.NT$NT.k <- (dummy.NT$LLa - dummy.NT$SSa)/((dummy.NT$SSa * dummy.NT$LLd) - 
    (dummy.NT$LLa * dummy.NT$SSd))
dummy.NT$Choice <- ifelse(dummy.NT$Choice == 2, 1, 0)

dummy.NT.glm <- function(dummy.NT) {
    glm(dummy.NT$Choice ~ dummy.NT$NT.k, family = binomial)
}
dummy.NT.AIC <- function(dummy.NT) {
    extractAIC(glm(dummy.NT$Choice ~ dummy.NT$NT.k, family = binomial))
}

dummy.NT.glm <- dlply(dummy.NT, .(ID, Acc.D), dummy.NT.glm)
dummy.NT.AIC <- dlply(dummy.NT, .(ID, Acc.D), dummy.NT.AIC)

intercept <- sapply(dummy.NT.glm, `[[`, 1)

k <- data.frame((-intercept[1, ])/intercept[2, ])

After k is calculated we need to calculate the difference in k values between groups First I resutrcture a dataframe to fit the data

AIC <- data.frame(sapply(dummy.NT.AIC, `[[`, 2))
AIC$ID <- rownames(AIC)
AIC$ID <- str_extract(AIC$ID, "^[0-9]{4}")
AIC$group <- str_extract(AIC$ID, "^[0-9]{1}")

ITC.n <- cbind(AIC, k)
ITC.n <- rename(ITC.n, c(sapply.dummy.NT.AIC........2. = "AIC", X..intercept.1.....intercept.2... = "k"))

head(ITC.n)
##          AIC   ID group      k
## 1002.1 15.01 1002     1 1.6391
## 1002.2 15.29 1002     1 1.6372
## 1003.1 17.97 1003     1 3.2484
## 1003.2 31.13 1003     1 5.6103
## 1006.1  4.00 1006     1 0.7993
## 1006.2 25.82 1006     1 3.0610

Now I do an ANOVA to see if there are group differences

ITC.n$group <- as.factor(ITC.n$group)

summary(aov(k ~ group, ITC.n))
##              Df   Sum Sq  Mean Sq F value Pr(>F)
## group         2 4.02e+32 2.01e+32     0.3   0.74
## Residuals   130 8.60e+34 6.62e+32