################### ### Ferni 1 ### ################### library(nlme) source('summarySE.r') library(lme4) long1b = read.csv(file = 'ferni1.csv', header = T) ctrl <- lmeControl(opt='optim'); ### Analysis of AT scores atMod1.3 = lme(AT ~ 1 + abs_dev + log_dev, random = ~ 1 + abs_dev + log_dev|VP, data = long1b, method = 'ML') summary(atMod1.3) ### Analysis of change relativized by intial estimates as a function of distance relativized by initial estimates shiftMod1.1 = lme(shift ~ 1 + abs_dev, random = ~1+abs_dev|VP, data=long1b, method="ML", na.action = na.omit, control = ctrl) shiftMod1.2 = lme(shift ~ 1 + abs_dev + log_dev, random = ~1+abs_dev+log_dev|VP, data=long1b, method="ML", na.action = na.omit, control = ctrl) shiftMod1.3 = lme(shift ~ 1 + abs_dev + abs_dev_sq, random = ~1+abs_dev+abs_dev_sq|VP, data=long1b, method="ML", na.action = na.omit, control = ctrl) anova(shiftMod1.1, shiftMod1.2) anova(shiftMod1.1, shiftMod1.3) summary(shiftMod1.3) # best-fitting model ### proportion of AT = 0 vs advice distance long1b$ignore = 1 long1b$ignore[long1b$AT > 0] = 0 ignoreMod2 = suppressWarnings(glmer(ignore ~ abs_dev + log_dev + (1 + abs_dev + log_dev | VP), data = long1b, family = binomial, nAGQ = 1)) summary(ignoreMod2) ### confidence shifts dcMod1.3 = lme(DC ~ abs_dev+log_dev, random = ~ 1 + abs_dev + log_dev|VP, data = long1b, method = 'ML', control = ctrl, na.action = na.omit) summary(dcMod1.3) ### preparation of Plots long1b$fitted = fitted(ignoreMod2) ignore = rep(0,15) meanAT1 = rep(0,15) meanDC1 = rep(0,15) meanDEV1 = rep(0,15) meanShift1 = rep(0,15) j=1 for (i in unique(long1b$abs_dev)) { a = subset(long1b, abs_dev == i) meanAT1[j] = mean(a$AT, na.rm = T) meanDC1[j] = mean(a$DC, na.rm = T) meanDEV1[j] = i meanShift1[j] = mean(a$shift, na.rm = T) ignore[j] = mean(a$fitted, na.rm = T) j = j+1 } shiftInt = fixef(shiftMod1.3)[1] + ranef(shiftMod1.3)[1] shiftSlope = fixef(shiftMod1.3)[2] + ranef(shiftMod1.3)[2] shiftSlopeSq = fixef(shiftMod1.3)[3] + ranef(shiftMod1.3)[3] atInt = fixef(atMod1.3)[1] + ranef(atMod1.3)[1] atSlope = fixef(atMod1.3)[2] + ranef(atMod1.3)[2] atSlopeLog = fixef(atMod1.3)[3] + ranef(atMod1.3)[3] dcInt = fixef(dcMod1.3)[1] + ranef(dcMod1.3)[1] dcSlope = fixef(dcMod1.3)[2] + ranef(dcMod1.3)[2] dcSlopeLog = fixef(dcMod1.3)[3] + ranef(dcMod1.3)[3] plotShift = function(dev) { shift = fixef(shiftMod1.3)[1] + fixef(shiftMod1.3)[2] * dev + fixef(shiftMod1.3)[3] * (dev^2) plot(dev, shift, type = 'l', ylim = c(0,0.5), xlim = c(0, 80), axes = F, ylab = NA, xlab = NA, pin = c(6,4)) par(new = F) } plotAT = function(dev) { at = fixef(atMod1.3)[1] + fixef(atMod1.3)[2] * dev + fixef(atMod1.3)[3] * log(dev) plot(dev, at, type = 'l', ylim = c(0,0.5), xlim = c(0, 80), axes = F, ylab = NA, xlab = NA, pin = c(6,4)) par(new = F) } plotDC = function(dev) { dc = fixef(dcMod1.3)[1] + fixef(dcMod1.3)[2] * dev + fixef(dcMod1.3)[3] * log(dev) plot(dev, dc, type = 'l', ylim = c(-1,2), xlim = c(0,80), axes = F, ylab = NA, xlab = NA, pin = c(6,4)) par(new = F) } shiftEst = function(Int, Slope, SlopeSq, Dev) { shift = Int + Slope*Dev + SlopeSq*(Dev^2) plot(Dev, shift, type = 'l', ylim = c(0,0.5), xlim = c(0,80), axes = F, ylab = NA, xlab = NA, pin = c(6,4), col = '#BBBBBB') par(new = T) } ATest = function(Int, Slope, SlopeLog, Dev) { at = Int + Slope*Dev + SlopeLog*log(Dev) plot(Dev, at, type = 'l', ylim = c(0,0.5), xlim = c(0,80), axes = F, ylab = NA, xlab = NA, pin = c(6,4), col = '#BBBBBB') par(new = T) } DCest = function(Int, Slope, SlopeLog, Dev) { dc = Int + Slope*Dev + SlopeLog*log(Dev) plot(Dev, dc, type = 'l', ylim = c(-1,2), xlim = c(0,80), axes = F, ylab = NA, xlab = NA, pin = c(6,4), col = '#BBBBBB') par(new = T) } iPlotDataVP = suppressWarnings(summarySE(long1b, measurevar = 'fitted', groupvars = c('abs_dev', 'VP'))) iPlotDataMean = summarySE(iPlotDataVP, measurevar = 'fitted', groupvars = 'abs_dev') ### final plot - all 4 in 1 plot1 = function(){ layout(t(matrix(c(1:4), nrow = 2))) par(mar = c(3,3,.5,.5)) for (j in 1:length(atInt[,1])) { ATest(atInt[j,1], atSlope[j,1], atSlopeLog[j,1], seq(5,75,5)) } plot(meanDEV1, meanAT1, ylim = c(0,0.5), xlim = c(0,80), axes = F, ylab = NA, xlab = NA, pin = c(6,4)) par(new=T) plotAT(seq(5,75,5)) box() axis(1, tck = -.020, labels = NA) axis(2, tck = -.020, labels = NA) axis(1, lwd = 0, line = -.7, cex.axis = .8) axis(2, lwd = 0, line = -.7, cex.axis = .8) mtext(side = 1, 'advice distance', line = 1.5, cex = 1.2) mtext(side = 2, 'advice weighting', line = 1.5, cex = 1.2) par(mar = c(3,3,.5,.5)) for (j in 1:length(shiftInt[,1])) { shiftEst(shiftInt[j,1], shiftSlope[j,1], shiftSlopeSq[j,1], seq(5,75,5)) } plot(meanDEV1, meanShift1, ylim = c(0,0.5), xlim = c(0,80), axes = F, ylab = NA, xlab = NA, pin = c(6,4)) par(new=T) plotShift(seq(5,75,5)) box() axis(1, tck = -.020, labels = NA) axis(2, tck = -.020, labels = NA) axis(1, lwd = 0, line = -.7, cex.axis = .8) axis(2, lwd = 0, line = -.7, cex.axis = .8) mtext(side = 1, 'advice distance', line = 1.5, cex = 1.2) mtext(side = 2, 'absolute opinion shift', line = 1.5, cex = 1.2) par(mar = c(3,3,.5,.5)) for(i in unique(long1b$VP)){ plot(iPlotDataVP$abs_dev[iPlotDataVP$VP == i], iPlotDataVP$fitted[iPlotDataVP$VP == i]*100, ylim = c(0,100), xlim = c(0,80), axes = F, ylab = NA, xlab = NA, pin = c(6,4), type = 'l', col = '#BBBBBB') par(new = T) } plot(meanDEV1, ignore*100, ylim = c(0,100), xlim = c(0,80), axes = F, ylab = NA, xlab = NA, pin = c(6,4)) par(new=T) plot(iPlotDataMean$abs_dev, iPlotDataMean$fitted*100, ylim = c(0,100), xlim = c(0,80), type = 'l', , axes = F, ylab = NA, xlab = NA, pin = c(6,4)) box() axis(1, tck = -.020, labels = NA) axis(2, tck = -.020, labels = NA) axis(1, lwd = 0, line = -.7, cex.axis = .8) axis(2, lwd = 0, line = -.7, cex.axis = .8) mtext(side = 1, 'advice distance', line = 1.5, cex = 1.2) mtext(side = 2, 'advice ignored (in %)', line = 1.5, cex = 1.2) par(mar = c(3,3,.5,.5)) for (j in 1:length(dcInt[,1])) { DCest(dcInt[j,1], dcSlope[j,1], dcSlopeLog[j,1], seq(5,75,5)) } plot(meanDEV1, meanDC1, ylim = c(-1,2), xlim = c(0,80), axes = F, ylab = NA, xlab = NA, pin = c(6,4)) par(new=T) plotDC(seq(5,75,5)) abline (0,0) box() axis(1, tck = -.020, labels = NA) axis(2, tck = -.020, labels = NA) axis(1, lwd = 0, line = -.7, cex.axis = .8) axis(2, lwd = 0, line = -.7, cex.axis = .8) mtext(side = 1, 'advice distance', line = 1.5, cex = 1.2) mtext(side = 2, 'confidence shift', line = 1.5, cex = 1.2) } postscript('Fig3_6x6.eps', width = 6, height = 6) plot1() dev.off() png('Fig3_6x6.png', width = 6, height = 6, units = 'in', res = 300) plot1() dev.off()