09 Feb Replication Wilke 2014

Replication of Wilke (2014)

Eight of 10 main effects of
domain and gender have the same direction as
in Wilke et al. (2014). Eight of nine effects had
the same direction in both studies, and the ninth
was small (and insignificant) in both studies

The tables can be found in Appendix A of the paper.

Setup

# Libraries & Directories
library(reshape2)   #To change from long to short format
library(data.table) #To handle data.frames more efficiently
library(lsr)        #For effect sizes
library(pwr)
library(xtable)     #For table outputs
library(scales)     #Percentage format
# Helper functions: download the utils folder and maybe change the path path argument
sapply(list.files(path="utils/", pattern="*.R", full.names = TRUE), source, .GlobalEnv)

# order of the risk domains
rskDom <- c("bgc","wgc","sta","env","fac","fse","par","kin","mat","mre")

Read data of current study

setwd("../4-Data") # Change this to your data directory
d   <- fread("risk.csv", key="i,fem") # risk likelihoods
dem <- fread("demographics.csv", drop = c("staDat","endDat","totTim"), key="i,fem") # demographics without timing variables
d <- d[dem,on = c("i","fem")] # Merge risk data and demographics data
rm(dem) # clean workspace
d[, source := "Present study"]

Read data by Wilke et al. (2014)

Note: the next data we will read in has to be obtained from the authors of
Wilke, A., Sherman, A., Curdt, B., Mondal, S., Fitzgerald, C., & Kruger, D. J. (2014). An evolutionary domain-specific risk scale. Evolutionary Behavioral Sciences, 8, 123â141. doi:10.1037/ebs0000011

dw2 <- load_and_preprocess_study2(filename = "wilkeStudy2.sav")
dw3 <- load_and_preprocess_study3(filename = "wilkeStudy3.sav")
# Make the order identical
setcolorder(dw2, names(d))
setcolorder(dw3, names(d))

d <- rbind(d,dw2,dw3) # Merge data
rm(dw2,dw3) # clean workspace

# Make the variables identical to the variables in d
d[, fem := factor(ifelse(fem=="Female",1,0), labels=c("Male","Female"))]
d[, mar := factor(ifelse((mar=="Single" |  mar == "No" | mar == 0), 0, 1), levels = 0:1, labels=c("Single","Relationship"))]
d[, bir := factor(bir, levels = c("First","Middle","Last","Other (step-sid or twins)"), ordered=F)]
d[, dom := factor(dom, levels=rskDom, labels=toupper(rskDom), ordered=F)]
d[source != "Study 2 by Wilke et al. (2014)", lik.f := factor(lik, levels=as.character(1:7), labels=c("Extremely Unlikely","Moderately Unlikely","Somewhat Unlikely","Not Sure","Somewhat Likely","Moderately Likely","Extremely Likely"), ordered=TRUE)]
d[source == "Study 2 by Wilke et al. (2014)", lik.f := factor(lik, levels=c("Extremely Unlikely","Moderately Unlikely","Somewhat Unlikely","Not sure","Somewhat Likely","Moderately Likely","Extremely Likely"), labels=c("Extremely Unlikely","Moderately Unlikely","Somewhat Unlikely","Not Sure","Somewhat Likely","Moderately Likely","Extremely Likely"), ordered=TRUE)]

setwd("../5-Code") # Reset the working directory

Domain effects

Run ordinal multinomial regression as mixed model with participants as random effects to test for effects of domains and/or gender on risk-taking likelihood. Intuitively, this model estimates the probability to check the next-higher response category (to surpass a hypothetical threshhold for the next-higher possible response). The model returns the log odds of responding with the next-higher response category given the predictors. Note, it does not return the probability directly.

# Mixed ordered multinomial logistic regression
library(ordinal)
# Fit a model with domain x gender as independent variables
res.d.dom.fem <- clmm(lik.f ~ dom * fem + (1 | i), data=d[source == "Present study"])
# Fit a model with  only ~ gender as independent variable
res.d.fem <- clmm(lik.f ~ fem + (1 | i), data=d[source == "Present study"])

# Compare the two models
anova(res.d.dom.fem, res.d.fem)

# Which effects are there
summary(res.d.dom.fem)

Comparison of Demographics

Compare the demographics between Wilke's studies and the present study

# Significance tests compared to study 3
# Change "Study 3 by Wilke et al. (2014)" to "Study 2 by Wilke et al. (2014)" if you want to see results compared to Study 2
d[!duplicated(i), t.test(age[source == "Present study"],
                         age[source == "Study 3 by Wilke et al. (2014)"])]
d[!duplicated(i), cohensD(age[source == "Present study"],
                          age[source == "Study 3 by Wilke et al. (2014)"], method="unequal")]
# Number of offspring
d[!duplicated(i), t.test(offNow[source == "Present study"],
                         offNow[source == "Study 3 by Wilke et al. (2014)"])]
d[!duplicated(i), cohensD(offNow[source == "Present study"],
                          offNow[source == "Study 3 by Wilke et al. (2014)"], method="unequal")]

# Marital status
tmp <- d[!duplicated(i), as.list(table(mar)), by = source][c(1,3), 3:2] # counts
prop.test(as.matrix(tmp))
ES.h(tmp[1,1]/sum(tmp[1,]), tmp[2,1]/sum(tmp[2,]))

Comparison of Effects

We compare the domain effects that resulted from the ordinal regression model above with the effects from study 3 and study 2 by Wilke and colleagues.

# Wilke et al (2014)'s data, fit the same model
# Study 3
res.dw3.dom.fem <- clmm(lik.f ~ dom * fem + (1 | i), data=d[source == "Study 3 by Wilke et al. (2014)"])
# Study 2
res.dw2.dom.fem <- clmm(lik.f ~ dom * fem + (1 | i), data=d[source == "Study 2 by Wilke et al. (2014)"])

# Table with Effect, CI, p-value
# The table output will be saved in LaTeX format
# Present Study
source("utils/expanddomains.r")
library(xtable)
# Store effects, cnofidence intervalls, p-values of current study
e1  <-  summary(res.d.dom.fem)$coefficients[7:25, "Estimate"]
ci1 <- confint(res.d.dom.fem)[7:25,]
p1  <- summary(res.d.dom.fem)$coefficients[7:25,"Pr(>|z|)"]

for (study in 2:3)
{
  # Store effects, cnofidence intervalls, p-values of current study
  res.dw.dom.fem <- get(paste0("res.dw",study,".dom.fem")) # get model for study x
  ew  <-  summary(res.dw.dom.fem)$coefficients[7:25, "Estimate"]
  ciw <- confint(res.dw.dom.fem)[7:25,]
  pw  <- summary(res.dw.dom.fem)$coefficients[7:25,"Pr(>|z|)"]




  # Make the table
  tab.rsk.dom <- data.frame(cbind(
      # Present study 
      sprintf("%.2f", e1),
      apply(ci1, 1, format.ci),
      format.p(p1),
      # Wilke et al (2014)
      sprintf("%.2f", ew),
      apply(ciw, 1, format.ci),
      format.p(pw),
      # Whether the effect was replicated
      ifelse((sign(e1)==sign(ew) & p1<.05 & pw<.05) | (p1>=.05 & pw>=.05),
          "Yes", ifelse(!(sign(e1)==sign(ew) & ! ((p1<.05 & pw<.05) | (p1>=.05 & pw>=.05))), "",  "Dir"))
      ))
  # Cosmetics for the table
  nam <- rownames(tab.rsk.dom)
  nam <- gsub("dom","",nam)
  nam <- gsub(":", " $\\times$ ",nam, fixed=T)
  nam <- sapply(nam, function(x) paste0(expanddomains(substr(x,1,3)), substr(x, 4,nchar(x))))
  nam <- gsub("femFemale","Female",nam,fixed=T)
  tab.rsk.dom$domain <- nam
  tab.rsk.dom$variable <- c("\\multirow{9}{*}{\\rotatebox[origin=c]{90}{Domain}}",rep("",8),"\\rotatebox[origin=c]{90}{Gen.}","\\multirow{9}{*}{\\rotatebox[origin=c]{90}{Domain $\\times$ Gender}}",rep("",8))
  tab.rsk.dom <- tab.rsk.dom[,c(9,8,1:7)]

  # Print table
  sanitize <- function(x) { gsub("<","$<$",x,fixed=T) }

  print(file=paste0("Appendix_table_replication_study", study, ".tex"),
      xtable(tab.rsk.dom,
            label=paste0("tab.rsk.dom",study),
            digits = c(0,0,2,2,2,2,2,2,2,0)),
      include.colnames = FALSE,
      only.contents=TRUE,
      include.rownames=F,
      sanitize.text.function = sanitize,
      hline.after=c(9,10),
      math.style.negative=TRUE)
}