Create a gist now

Instantly share code, notes, and snippets.

## 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment