lee_schnalletal2008_analysis.R

## analysis code for Fiona Lee's class project
## Psych 241, Stanford Winter 2013
library(ggplot2)
library(reshape2)
library(psych)
library(bootstrap)
library(lme4)

## for bootstrapping 95% confidence intervals
theta <- function(x,xdata,na.rm=T) {mean(xdata[x],na.rm=na.rm)}
ci.low <- function(x,na.rm=T) {
  mean(x,na.rm=na.rm) - quantile(bootstrap(1:length(x),1000,theta,x,na.rm=na.rm)$thetastar,.025,na.rm=na.rm)}
ci.high <- function(x,na.rm=T) {
  quantile(bootstrap(1:length(x),1000,theta,x,na.rm=na.rm)$thetastar,.975,na.rm=na.rm) - mean(x,na.rm=na.rm)}

####### BEGIN ANALYSIS ######
d <- read.csv("http://www.stanford.edu/~mcfrank/lee_data.csv")

## Exclusions
length(d$subid)
sum(d$GuessedGoalsStrict==1)
sum(d$AttentionCheckFail==1)
d <- subset(d, GuessedGoalsStrict==0 & AttentionCheckFail==0)
length(d$subid)

## Plot emotion ratings
emo.vars <- names(d)[grepl("emo",names(d))]
emo.md <- melt(d, id.vars=c("subid","condition"), 
           measure.vars=emo.vars)
qplot(value, fill=condition, geom="histogram", 
      facets=~variable, position="dodge", binwidth=1,
      data=emo.md)
      
## emotion rating analysis
rs.relaxed = lm(emo.relaxed ~ condition, d); summary(rs.relaxed)
rs.angry = lm(emo.angry ~ condition, d); summary(rs.angry)
rs.happy = lm(emo.happy ~ condition, d); summary(rs.happy)
rs.sad = lm(emo.sad ~ condition, d); summary(rs.sad)
rs.relaxed = lm(emo.relaxed ~ condition, d); summary(rs.relaxed)
rs.afraid = lm(emo.afraid ~ condition, d); summary(rs.afraid)
rs.depressed = lm(emo.depressed ~ condition, d); summary(rs.depressed)
rs.disgusted = lm(emo.disgusted ~ condition, d); summary(rs.disgusted)
rs.upset = lm(emo.upset ~ condition, d); summary(rs.upset)
rs.confused = lm(emo.confused ~ condition, d); summary(rs.confused)

## distribution of morality judgments
moral.vars <- names(d)[grepl("moral",names(d))]
moral.md <- melt(d, id.vars=c("subid","condition","age"), 
                 measure.vars=moral.vars)
qplot(value, fill=condition, geom="histogram", 
      facets=~variable, position="dodge", binwidth=1,
      data=moral.md) 

## Reliability and composite
alpha(d[,moral.vars])
moral.comp <- rowMeans(d[,moral.vars])

## Basic moral judgments plot
ms <- aggregate(moral.comp ~ condition, d, mean)
ms$cih <- aggregate(moral.comp ~ condition, d, ci.high)$moral.comp
ms$cil <- aggregate(moral.comp ~ condition, d, ci.low)$moral.comp
qplot(condition, moral.comp, ymin=moral.comp - cil, ymax = moral.comp + cih, 
      geom=c("bar","linerange"),
      data=ms)

## by scenario
mss <- aggregate(value ~ condition + variable, moral.md, mean)
mss$cih <- aggregate(value ~ condition + variable, moral.md, ci.high)$value
mss$cil <- aggregate(value ~ condition + variable, moral.md, ci.low)$value
qplot(condition, value, ymin=value - cil, ymax = value + cih, 
      geom=c("bar","linerange"), facets=~variable,
      data=mss)

## moral judgments analysis 
rs.comp = lm(moral.comp ~ condition, d); summary(rs.comp)
rs.comp.lmer = lmer(value ~ condition  + (condition | variable) + (1|subid), 
                    data=moral.md)

## individual items
rs.dog = lm(moral.dog ~ condition, d); summary(rs.dog)
rs.wallet = lm(moral.wallet ~ condition, d); summary(rs.wallet)
rs.trolley = lm(moral.trolley ~ condition, d); summary(rs.trolley)
rs.plane = lm(moral.plane ~ condition, d); summary(rs.plane)
rs.resume = lm(moral.resume ~ condition, d); summary(rs.resume)
rs.kitten = lm(moral.kitten ~ condition, d); summary(rs.kitten)

## interaction with age, post-hoc
quartz()
qplot(age,moral.comp, col=condition,
      data=d) + 
  ylab("Moral Judgment Composite") + 
  xlab("Age") + 
  geom_smooth(method="lm")

rs.comp.age = lm(moral.comp ~ age, d); summary(rs.comp.age)
rs.comp = lm(moral.comp ~ condition + age, d); summary(rs.comp)
rs.comp.age = lm(moral.comp ~ condition * age, d); summary(rs.comp.age)

# note: LMER doesn't converge with a condition*age random effect
rs.comp.lmer = lmer(value ~ condition * age  + (condition + age | variable) + (1|subid), 
                    data=moral.md)