Re-analysis if volunteer data

volunteer_reanalysis.R

# Code from https://osf.io/aqi5j/ made by original authors, adapted by Daniel Lakens
# Download the csv file from https://osf.io/aqi5j/
# Paper: http://www.tandfonline.com/doi/full/10.1080/23743603.2016.1273647

require(car)

# IMPORTANT: change this to match your path/file
fileLocation <- 'C:\\Users\\Daniel\\surfdrive\\Data\\blog volunteer study\\VolunteeringDataFile_05.09.16.csv'

# Dienes functions
# Author: Zoltan Dienes

#Bf function calculates the Bayes factor as in Dienes calculator
#uniform can take values 0 or 1. 
#If it is 1, specify lower and upper values.
#If it is 0, specify tail=2 (two-tailed) or tail=1 (one-tailed, half normal)
Bf<-function(sd, obtained, uniform, lower=0, upper=1, meanoftheory=0, sdtheory=1, tail=2)
{
#Authors Danny Kaye & Thom Baguley
#Version 1.0
#19/10/2009
#test data can be found starting at p100
#
    area <- 0
    if(identical(uniform, 1)){
        theta <- lower
        range <- upper - lower
        incr <- range / 2000
        for (A in -1000:1000){
            theta <- theta + incr
            dist_theta <- 1 / range
            height <- dist_theta * dnorm(obtained, theta, sd)
            area <- area + height * incr
        }
    }else{
        theta <- meanoftheory - 5 * sdtheory
        incr <- sdtheory / 200
        for (A in -1000:1000){
            theta <- theta + incr
            dist_theta <- dnorm(theta, meanoftheory, sdtheory)
            if(identical(tail, 1)){
                if (theta <= 0){
                    dist_theta <- 0
                } else {
                    dist_theta <- dist_theta * 2
                }
            }
            height <- dist_theta * dnorm(obtained, theta, sd)
            area <- area + height * incr
        }
    }
    LikelihoodTheory <- area
    Likelihoodnull <- dnorm(obtained, 0, sd)
    BayesFactor <- LikelihoodTheory / Likelihoodnull
    ret <- list("LikelihoodTheory" = LikelihoodTheory, "Likelihoodnull" = Likelihoodnull, "BayesFactor" = BayesFactor)
    ret
}

# study analysis starts here

# load file
srcData <- read.csv(fileLocation)

# double checking missingness in data file
apply(srcData, 2, function(x){all(complete.cases(x))})

# changing categorical iv to factor due to ocd
srcData$Condition <- as.factor(srcData$Condition)

# reweighting the composite score
# this doesn't really matter for the final result since we're using standardized effects
srcData$SWBx2_T2<-srcData$SWB_T2
srcData$SWBx2_T1<-srcData$SWB_T1
srcData$SWB_T2<-srcData$SWB_T2/2
srcData$SWB_T1<-srcData$SWB_T1/2

# clean out the existing t2t1 columns since we're recreating them
varlist <- names(srcData)
varlist <- varlist[-grep('T2T1', varlist)]

# get all the names for t2 vars and t1 vars
t2vars <- varlist[grep('T2_|_T2', varlist)]
t1vars <- varlist[grep('T1_|_T1', varlist)]

# get the variable names and creating new names for t2t1
varnames <- gsub('T2_', '', gsub('_T2', '', t2vars))
t2t1names <- paste0('T2T1_', varnames)

# subtract t1 from t2 for all variables
for(i in 1:length(t2vars)) {
    srcData[,t2t1names[i]] <- srcData[,t2vars[i]] - srcData[,t1vars[i]]
}

results <- list()

# loop through t2t1 variables
for(var0 in t2t1names){

    # checking if we got all the variables we want
    cat(paste0("Processing ", var0, "...\n"))

    # hist(data[data$condition==1,2])
    # hist(data[data$condition==2,2])

    # testing for hov
    hov <- leveneTest(as.formula(paste0(var0, "~Condition")), srcData)['group', 'Pr(>F)']

    # group n
    n1 <- summary(srcData$Condition)['1']
    n2 <- summary(srcData$Condition)['2']

    # group means
    m1 <- mean(srcData[srcData$Condition==1,var0])
    m2 <- mean(srcData[srcData$Condition==2,var0])

    # group sd
    sd1 <- sd(srcData[srcData$Condition==1,var0])
    sd2 <- sd(srcData[srcData$Condition==2,var0])

    # pooled sd
    #NOTE DL: Original file had typo, n2-2 instead of n2-1 in sd formula
    sdp <- sqrt((((n1-1)*sd1^2 + (n2-1)*sd2^2)/(n1+n2-2)))

    # sem
    semp <- sdp*sqrt(1/n1+1/n2)

    # mean diff
    meandiff <- m2 - m1

    # numbers we need for the dienes app
    results[[var0]] <- c()
    
    # uncomment to include hov test results
    # results[[var0]]['hov'] <- hov

    # the "raw" numbers. these are probably not what we want
    # to compare to the N(.5, .15^2) priors
    # but they're included for completeness
    results[[var0]]['raw mean.'] <- meandiff
    results[[var0]]['pooled sd.'] <- sdp
    results[[var0]]['standard error.'] <- semp

    results[[var0]]['standardized mean.'] <- meandiff/sdp
    results[[var0]]['se of standardized mean.'] <- semp/sdp

    results[[var0]]['bf.'] <- Bf(semp/sdp, meandiff/sdp, F, meanoftheory=0.5, sdtheory=.15, tail=2)$BayesFactor
#NOTE DL: I'm also storing means, sd, and n for equivalence test. 
    results[[var0]]['m1'] <- m1
    results[[var0]]['m2'] <- m2
    results[[var0]]['sd1'] <- sd1
    results[[var0]]['sd2'] <- sd2
    results[[var0]]['n1'] <- n1
    results[[var0]]['n2'] <- n2
}

#Now we use this data to perform equivalence tests for all effects listed in Table 2

library("TOSTER")

#SWBX2
TOSTtwo(m1=results$T2T1_SWBx2[7],m2=results$T2T1_SWBx2[8],sd1=results$T2T1_SWBx2[9],sd2=results$T2T1_SWBx2[10],n1=results$T2T1_SWBx2[11],n2=results$T2T1_SWBx2[12],low_eqbound_d=-0.5,high_eqbound_d=0.5)
#Compare with BF
results$T2T1_SWBx2[6]

#SPANE_PA
TOSTtwo(m1=results$T2T1_SPANE_PA[7],m2=results$T2T1_SPANE_PA[8],sd1=results$T2T1_SPANE_PA[9],sd2=results$T2T1_SPANE_PA[10],n1=results$T2T1_SPANE_PA[11],n2=results$T2T1_SPANE_PA[12],low_eqbound_d=-0.5,high_eqbound_d=0.5)
#Compare with BF
results$T2T1_SPANE_PA[6]

#SPANE_NA
TOSTtwo(m1=results$T2T1_SPANE_NA[7],m2=results$T2T1_SPANE_NA[8],sd1=results$T2T1_SPANE_NA[9],sd2=results$T2T1_SPANE_NA[10],n1=results$T2T1_SPANE_NA[11],n2=results$T2T1_SPANE_NA[12],low_eqbound_d=-0.5,high_eqbound_d=0.5)
#Compare with BF
#NOTE DL: There is a Type in Table 2: Result for NA is 0.03, not 0.06
results$T2T1_SPANE_NA[6]

#SWLS
TOSTtwo(m1=results$T2T1_SWLS[7],m2=results$T2T1_SWLS[8],sd1=results$T2T1_SWLS[9],sd2=results$T2T1_SWLS[10],n1=results$T2T1_SWLS[11],n2=results$T2T1_SWLS[12],low_eqbound_d=-0.5,high_eqbound_d=0.5)
#Compare with BF
results$T2T1_SWLS[6]

#CESD
TOSTtwo(m1=results$T2T1_CESD[7],m2=results$T2T1_CESD[8],sd1=results$T2T1_CESD[9],sd2=results$T2T1_CESD[10],n1=results$T2T1_CESD[11],n2=results$T2T1_CESD[12],low_eqbound_d=-0.5,high_eqbound_d=0.5)
#Compare with BF
results$T2T1_CESD[6]

#PSS
TOSTtwo(m1=results$T2T1_PSS[7],m2=results$T2T1_PSS[8],sd1=results$T2T1_PSS[9],sd2=results$T2T1_PSS[10],n1=results$T2T1_PSS[11],n2=results$T2T1_PSS[12],low_eqbound_d=-0.5,high_eqbound_d=0.5)
#Compare with BF
results$T2T1_PSS[6]

#SCS
TOSTtwo(m1=results$T2T1_SCS[7],m2=results$T2T1_SCS[8],sd1=results$T2T1_SCS[9],sd2=results$T2T1_SCS[10],n1=results$T2T1_SCS[11],n2=results$T2T1_SCS[12],low_eqbound_d=-0.5,high_eqbound_d=0.5)
#Compare with BF
results$T2T1_SCS[6]

#UCLA
TOSTtwo(m1=results$T2T1_UCLA[7],m2=results$T2T1_UCLA[8],sd1=results$T2T1_UCLA[9],sd2=results$T2T1_UCLA[10],n1=results$T2T1_UCLA[11],n2=results$T2T1_UCLA[12],low_eqbound_d=-0.5,high_eqbound_d=0.5)
#Compare with BF
results$T2T1_UCLA[6]

#MLQ
TOSTtwo(m1=results$T2T1_MLQ[7],m2=results$T2T1_MLQ[8],sd1=results$T2T1_MLQ[9],sd2=results$T2T1_MLQ[10],n1=results$T2T1_MLQ[11],n2=results$T2T1_MLQ[12],low_eqbound_d=-0.5,high_eqbound_d=0.5)
#Compare with BF
results$T2T1_MLQ[6]