Example 3 for Wright, Strubler & Vallano
(2011) Legal and Criminological Psychology
This page goes through the R code and the R output in more detail than the paper. As time goes on, if people want more done we can add more here. Following the practice in the paper, dark will be R code and light will be output. Comments are written in purple. Here is a brief outline which allows you to go where you want in this page.
Wrongly running a non-multilevel model on multilevel data
Doing the correct multilevel model
Bootstrapping dispersion parameter
Graphing means for simulations
Graphing standard errors from simulations
This reads the data in from the web.
> example3 <-
read.table("http://www.fiu.edu/~dwright/jurstat/ex3jury.dat", header=TRUE)
> attach(example3)
>
> #sit 1 - binary multilevel logistic regression
> library(lme4)
Loading required package: Matrix
Loading required package: lattice
Attaching package: 'Matrix'
> lev1 <- glm(guilty ~ type+author,binomial)
> summary(lev1)
This shows statistically significant effects for both type of crime and authoritarianism. But, it assumes that the residuals are independent. They aren't, so the p values are wrong.
Call:
glm(formula = guilty ~ type + author, family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.6382 -1.1398 0.8087 1.1026 1.5728
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.902898 0.332360
-2.717 0.00660 **
type
0.792817 0.292764 2.708 0.00677 **
author
0.011533 0.005116 2.254 0.02417 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 276.94 on 199 degrees of freedom
Residual deviance: 263.68 on 197 degrees of freedom
AIC: 269.68
Number of Fisher Scoring iterations: 4
> anova(update(lev1, .~. -type),lev1,test="Chisq")
The update function is really useful. It saves typing. For some objects types you have to specify the test for ANOVA.
Analysis of Deviance Table
Model 1: guilty ~ author
Model 2: guilty ~ type + author
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 198 271.149
2 197 263.678 1
7.471 0.006
> anova(update(lev1, .~. -author),lev1,test="Chisq")
Analysis of Deviance Table
Model 1: guilty ~ type
Model 2: guilty ~ type + author
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 198 268.871
2 197 263.678 1
5.193 0.023
>
mod1 <- lmer(guilty ~ type + author + (1|juryno),family=binomial)
> summary(mod1)
This is the more appropriate multilevel analysis. Laplace approximation is used. You used to have to request this because it was slightly slower than other methods, but the speed differences was slight and it was a lot more accurate, so Doug Bates made it the only method (though they are working on further advances).
Generalized linear mixed
model fit by the Laplace approximation
Formula: guilty ~ type + author + (1 | juryno)
AIC BIC logLik deviance
241.9 255.1 -116.9 233.9
Random effects:
Groups Name Variance Std.Dev.
juryno (Intercept) 2.1005 1.4493
Number of obs: 200, groups: juryno, 20
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.060110 0.602051 -1.761
0.0783 .
type 1.134075
0.737541 1.538 0.1241
author 0.012161 0.005822
2.089 0.0367 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) type
type -0.616
author -0.491 -0.008
> anova(update(mod1, .~. -type),mod1)
Carefully using the update function when changing the random variables because +(author|jury) would add the random intercept also.
Data:
Models:
c("update", "mod1", ". ~ . - type"): guilty ~ author + (1 | juryno)
mod1: guilty ~ type + author + (1 | juryno)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
c("update", "mod1", ". ~ . - type") 3 242.13 252.03 -118.06
mod1 4 241.89 255.08 -116.94 2.2398 1 0.1345
> anova(update(mod1, .~. -author),mod1)
Data:
Models:
c("update", "mod1", ". ~ . - author"): guilty ~ type + (1 | juryno)
mod1: guilty ~ type + author + (1 | juryno)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
c("update", "mod1", ". ~ . - author") 3 244.07 253.96 -119.03
mod1 4 241.89 255.08 -116.94 4.1765 1 0.04099 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
> par(mfrow=c(1,3))
> # Figure 7a
> jauthor <- split(author,juryno)
> jtype <- split(type, juryno)
> jguilty <- split(guilty,juryno)
> plot(c(133,993),c(0,1),xlim=c(0,100),ylim=c(0,1),xlab="Authoritarianism",
+ ylab="Probability guilty",cex.axis=1.6,cex.lab=1.6,
+ main="Dashed misdemeanor. Solid felony.",cex.main=1.3)
> for (i in 1:20){
+ jaut <- jauthor[[i]]; jguil <- jguilty[[i]]
+ modjury <- glm(jguil ~ jaut, binomial)
+ lines(0:100,predict(modjury,data.frame(jaut=0:100),type="response"),lty=2-mean(jtype[[i]]))
+ points(mean(jaut),predict(modjury,data.frame(jaut=mean(jaut)),type="response"),
+ cex=1.3,pch=21,bg="grey")
+ }
This error was due to probabilities very near zero.
Warning message:
In glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart,
:
fitted probabilities numerically 0 or 1 occurred
>
>
> # Figure 7b
> mod2 <- lmer(guilty~author + (1|juryno),family=binomial)
> plot(c(133,993),c(0,1),xlim=c(0,100),ylim=c(0,1),xlab="Authoritarianism",
+ ylab="Probability guilty",cex.axis=1.6,cex.lab=1.6,main="Individual curves,
multilevel")
> datamod2 <- data.frame(cbind(rep(1,101),0:100))
> for (i in 1:20){
+ coefmodi <- fixef(mod2) + c(ranef(mod2)$juryno[i,],0)
+ preds <- coefmodi%*%t(datamod2)
+ predresponse <- exp(preds)/(1+exp(preds))
+ lines(0:100,predresponse)
+ }
>
>
> #Figure 7c
> plot(c(133,993),c(0,1),xlim=c(0,100),ylim=c(0,1),xlab="Authoritarianism",
+ ylab="Probability guilty",cex.axis=1.6,cex.lab=1.6,main="Group curve with +-
sd")
> datamod2 <- data.frame(cbind(rep(1,101),0:100))
> coefmodup <- fixef(mod2) + c(1.55,0)
> coefmoddown <- fixef(mod2) - c(1.55,0)
There is no predict function for the object from lmer, so this is me doing it.
> preds <- fixef(mod2)%*%t(datamod2)
> predsup <- coefmodup%*%t(datamod2)
> predsdown <- coefmoddown%*%t(datamod2)
> predresponse <- exp(preds)/(1+exp(preds))
> predresponseup <- exp(predsup)/(1+exp(predsup))
> predresponsedown <- exp(predsdown)/(1+exp(predsdown))
> polygon(c(0:100,100:0),c(predresponsedown,
+ sort(predresponseup,decreasing=TRUE)),col="grey80",border=NA)
> lines(0:100,predresponse,lwd=2)
>
> par(mfrow=c(1,1)) >
>
>
> # sit 2
>
> # or use aggregate
The aggregate function was written to avoid doing multiple tapply functions, but given these all have different statistics it is easier to keep track of them with tapply.
> medincome <- tapply(household,juryno,median)
> nguilty <- tapply(guilty,juryno,sum)
> ntype <- tapply(type,juryno,mean)
> njurors <- tapply(guilty,juryno,length)
>
>
ninnocent <- njurors - nguilty
> verdict <- cbind(nguilty,ninnocent)
> mod1 <- glm(verdict ~ ntype + medincome,binomial)
> summary(mod1)
Call:
glm(formula = verdict ~ ntype + medincome, family = binomial)
Deviance Residuals:
Min 1Q Median
3Q Max
-3.0378 -1.4432 -0.4835 1.8862 3.1581
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.961999 0.324195 -2.967
0.00300 **
ntype 0.754222
0.294212 2.564 0.01036 *
medincome 0.013843 0.005397
2.565 0.01032 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 92.501 on 19 degrees of freedom
Residual deviance: 77.297 on 17 degrees of freedom
AIC: 124.03
Number of Fisher Scoring iterations: 5
> anova(mod1,test="Chisq")
Analysis of Deviance
Table
Model: binomial, link: logit
Response: verdict
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid.
Dev P(>|Chi|)
NULL 19 92.501
ntype 1 8.068
18 84.433 0.005
medincome 1 7.137 17
77.297 0.008
>
>
> mod2 <- glm(verdict~ ntype +
medincome,quasibinomial); summary(mod2)
The statistics behind quasi-likelihood methods are complex. At present they are not often used but this may change (hence their inclusion in the paper!).
Call:
glm(formula = verdict ~ ntype + medincome, family = quasibinomial)
Deviance Residuals:
Min 1Q Median
3Q Max
-3.0378 -1.4432 -0.4835 1.8862 3.1581
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.96200 0.64037 -1.502
0.151
ntype
0.75422 0.58115 1.298 0.212
medincome 0.01384 0.01066
1.298 0.211
(Dispersion parameter for quasibinomial family taken to be 3.901685)
Null deviance: 92.501 on 19 degrees of freedom
Residual deviance: 77.297 on 17 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 5
> anova(glm(verdict~ ntype,quasibinomial),mod2,test="Chisq")
Analysis of Deviance Table
Model 1: verdict ~ ntype
Model 2: verdict ~ ntype + medincome
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 18 84.433
2 17 77.297 1
7.137 0.176
> anova(glm(verdict~ medincome,quasibinomial),mod2,test="Chisq")
Analysis of Deviance Table
Model 1: verdict ~ medincome
Model 2: verdict ~ ntype + medincome
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 18 83.970
2 17 77.297
1 6.673 0.191
> anova(glm(verdict~
1,quasibinomial),glm(verdict~ ntype,quasibinomial),test="Chisq")
Analysis of Deviance
Table
Model 1: verdict ~ 1
Model 2: verdict ~ ntype
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 19 92.501
2 18 84.433 1
8.068 0.146
> anova(glm(verdict~ 1,quasibinomial),glm(verdict~
medincome,quasibinomial),test="Chisq")
Analysis of Deviance
Table
Model 1: verdict ~ 1
Model 2: verdict ~ medincome
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 19 92.501
2 18 83.970 1
8.531 0.138
> summary(glm(verdict~1,quasibinomial))
Call:
glm(formula = verdict ~ 1, family = quasibinomial)
Deviance Residuals:
Min 1Q Median
3Q Max
-3.8314 -1.4078 -0.1265 1.3290 3.6164
Coefficients:
Estimate
Std. Error t value Pr(>|t|)
(Intercept) 0.08004 0.27969 0.286
0.778
(Dispersion parameter for quasibinomial family taken to be 3.905202)
Null deviance: 92.501 on 19 degrees of freedom
Residual deviance: 92.501 on 19 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 3
>
> library(boot)
Attaching package: 'boot'
> vars <-
as.data.frame(cbind(nguilty,ninnocent))
> dispfunc <- function(x,i)
+ summary(glm(cbind(vars$nguilty[i],vars$ninnocent[i])~ 1,
+ quasibinomial))$dispersion
> disboot <- boot(vars,dispfunc,R=2000)
> boot.ci(disboot,type="bca")
BOOTSTRAP CONFIDENCE INTERVAL
CALCULATIONS
Based on 2000 bootstrap replicates
CALL :
boot.ci(boot.out = disboot, type = "bca")
Intervals :
Level BCa
95% ( 2.581, 6.074 )
Calculations and Intervals on Original Scale
>
> # Figure 8
> binmod <- glm(verdict ~ medincome,binomial)
> qbinmod <- glm(verdict ~ medincome,quasibinomial)
> predbin <- predict(binmod,data.frame(medincome=0:150),type="response",se.fit=TRUE)
> predqbin <- predict(qbinmod,data.frame(medincome=0:150),type="response",se.fit=TRUE)
> plot(0:150,predqbin$fit,xlab="Median income", ylab="Proportion guilty
verdicts",
+ cex.axis=1.3,cex.lab=1.3)
> polygon(c(0:150,150:0),c(predqbin$fit[1:151]-predqbin$se[1:151],
+ predqbin$fit[151:1]+predqbin$se[151:1]),col="grey85",border=NA)
> polygon(c(0:150,150:0),c(predbin$fit[1:151]-predbin$se[1:151],
+ predbin$fit[151:1]+predbin$se[151:1]),col="grey60",border=NA)
> lines(0:150,predqbin$fit,lwd=2)
> plot(0:150,predqbin$fit,xlab="Median income", ylab="Proportion guilty
verdicts",
+ cex.axis=1.3,cex.lab=1.3,ylim=c(0,1),col="white",las=1)
> polygon(c(0:150,150:0),c(predqbin$fit[1:151]-2*predqbin$se[1:151],
+ predqbin$fit[151:1]+2*predqbin$se[151:1]),col="grey85",border=NA)
> polygon(c(0:150,150:0),c(predbin$fit[1:151]-2*predbin$se[1:151],
+ predbin$fit[151:1]+2*predbin$se[151:1]),col="grey60",border=NA)
> lines(0:150,predqbin$fit,lwd=2)
> points(medincome,nguilty/(nguilty+ninnocent),pch=21,bg="grey90",cex=1.5)
> legend(86,.14,c("binomial","quasi-binomial"),fill=c("grey60","grey85"),cex=1.3)
> box()
>
>
>
The numbers for these simulations is set at 10 rather than 1000 so that it runs more quickly. The means and standard deviations and that stuff though come from the 1000 replications per condition simulations. The reps variable controls this.
>
> # Simulation 1
> # Sit JURY
> # JURY, with jury variable
> totjurors <- 200
> antilogit <- function(x) 1/(1+exp(-x))
> jurorno <- 1:totjurors; counter <- 0; reps <- 10
> sim1 <- numeric(reps*5*4*8)
> dim(sim1) <- c(reps*5*8,4)
> for (i in 1:reps){
+ if (i%%25 == 0) print(i)
+ for (jef in 0:4){
+ for (nj in c(5,10,15,20,30,50,100,200)){
+ counter <- counter + 1
+ juryno <- as.numeric(cut(jurorno,nj))
+ juryef1 <- rnorm(length(unique(juryno)),0,jef/2)
+ juryef <- rep(juryef1,as.numeric(table(juryno)))
+ x1a <- rnorm(length(unique(juryno)),0,1)
+ x1 <- rep(x1a,as.numeric(table(juryno)))
+ guilty <- rbinom(totjurors,1,antilogit(juryef+x1))
+ nog <- tapply(guilty,juryno,sum)
+ nong <- tapply(1-guilty,juryno,sum)
+ x1s <- glm(cbind(nog,nong)~x1a,binomial)$coefficients[2]
+ sim1[counter,] <- c(i,jef,nj,x1s)
+ }}}
> sim1 <- as.data.frame(sim1)
> colnames(sim1) <- c("reps","jef","nj","x1s")
> diff1 <- function(x) sum(abs(1-x))/length(x)
> s1sd <- tapply(sim1$x1s,list(sim1$jef,sim1$nj),sd)
> s1diff <- tapply(sim1$x1s,list(sim1$jef,sim1$nj),diff1)
> s1mean <- tapply(sim1$x1s,list(sim1$jef,sim1$nj),mean)
>
> # Simulation 2
> # jury, with juror variable
> totjurors <- 200
> antilogit <- function(x) 1/(1+exp(-x))
> jurorno <- 1:totjurors; counter <- 0; reps <- 10
> sim2 <- numeric(reps*5*4*8)
> dim(sim2) <- c(reps*5*8,4)
> for (i in 1:reps){
+ if (i%%25 == 0) print(i)
+ for (jef in 0:4){
+ for (nj in c(5,10,15,20,30,50,100,200)){
+ counter <- counter + 1
+ juryno <- as.numeric(cut(jurorno,nj))
+ juryef1 <- rnorm(length(unique(juryno)),0,jef/2)
+ juryef <- rep(juryef1,as.numeric(table(juryno)))
+ x2 <- rnorm(length(juryno))
+ guilty <- rbinom(totjurors,1,antilogit(juryef+x2))
+ nog <- tapply(guilty,juryno,sum)
+ nong <- tapply(1-guilty,juryno,sum)
+ x2a <- tapply(x2,juryno,mean)
+ x2s <- glm(cbind(nog,nong)~x2a,binomial)$coefficients[2]
+ sim2[counter,] <- c(i,jef,nj,x2s)
+ }}}
> sim2 <- as.data.frame(sim2)
> colnames(sim2) <- c("reps","jef","nj","x2s")
> diff1 <- function(x) sum(abs(1-x))/length(x)
> s2sd <- tapply(sim2$x2s,list(sim2$jef,sim2$nj),sd)
> s2diff <- tapply(sim2$x2s,list(sim2$jef,sim2$nj),diff1)
> s2mean <- tapply(sim2$x2s,list(sim2$jef,sim2$nj),mean)
>
>
> library(lme4)
You will get some warnings with the next two simulations. I looked through about 50 of times when there are warnings (there is not a way that I could find to flag warnings and save those in a particular way). These were ones with values really near 0 or 1 and mostly with 50 in each jury (which usually would not be done).
> # Simulation 3
> # Response juror
> # with jury variable
> # needs lme4 if that is not already active
> options(warn=1)
> totjurors <- 200
> antilogit <- function(x) 1/(1+exp(-x))
> jurorno <- 1:totjurors; counter <- 0; reps <- 10
> sim3 <- numeric(reps*5*4*8)
> dim(sim3) <- c(reps*5*8,4)
> for (i in 1:reps){
+ if (i%%25 == 0) print(i)
+ for (jef in 0:4){
+ for (nj in c(5,10,15,20,30,50,100,200)){
+ counter <- counter + 1
+ juryno <- as.numeric(cut(jurorno,nj))
+ juryef1 <- rnorm(length(unique(juryno)),0,jef/2)
+ juryef <- rep(juryef1,as.numeric(table(juryno)))
+ x3a <- rnorm(length(unique(juryno)),0,1)
+ x3 <- rep(x3a,as.numeric(table(juryno)))
+ guilty <- rbinom(totjurors,1,antilogit(juryef+x3))
+ x3s <- summary(lmer(guilty~x3+(1|juryno),family=binomial))@coefs[2,1]
+ sim3[counter,] <- c(i,jef,nj,x3s)
+ }}}
> write.table(sim3,"c:\\temp\\sim3.dat")
> sim3 <- as.data.frame(sim3)
> colnames(sim3) <- c("reps","jef","nj","x3s")
> options(warn=0)
> sim3 <- as.data.frame(sim3)
> diff1 <- function(x) sum(abs(1-x))/length(x)
> s3sd <- tapply(sim3$x3s,list(sim3$jef,sim3$nj),sd)
> s3diff <- tapply(sim3$x3s,list(sim3$jef,sim3$nj),diff1)
> s3mean <- tapply(sim3$x3s,list(sim3$jef,sim3$nj),mean)
>
>
>
> ############
>
>
> library(lme4)
> # Simulation 4
> # Response juror
> # with juror variable
> options(warn=1)
> totjurors <- 200
> antilogit <- function(x) 1/(1+exp(-x))
> jurorno <- 1:totjurors; counter <- 0; reps <- 10
> sim4 <- numeric(reps*5*4*8)
> dim(sim4) <- c(reps*5*8,4)
> for (i in 1:reps){
+ if (i%%25 == 0) print(i)
+ for (jef in 0:4){
+ for (nj in c(5,10,15,20,30,50,100,200)){
+ counter <- counter + 1
+ juryno <- as.numeric(cut(jurorno,nj))
+ juryef1 <- rnorm(length(unique(juryno)),0,jef/2)
+ juryef <- rep(juryef1,as.numeric(table(juryno)))
+ x4 <- rnorm(200,0,1)
+ guilty <- rbinom(totjurors,1,antilogit(juryef+x4))
+ x4s <- summary(lmer(guilty~x4+(1|juryno),family=binomial))@coefs[2,1]
+ sim4[counter,] <- c(i,jef,nj,x4s)
+ }}}
> write.table(sim4,"c:\\temp\\sim4.dat")
> sim4 <- as.data.frame(sim4)
> colnames(sim4) <- c("reps","jef","nj","x4s")
>
> options(warn=0)
> sim4 <- as.data.frame(sim4)
> s4sd <- tapply(sim4$x4s,list(sim4$jef,sim4$nj),sd)
> s4mean <- tapply(sim4$x4s,list(sim4$jef,sim4$nj),mean)
>
> ####################
>
>
>
>
> library(splines)
> s1mean <- c(
+ 1.0233492, 1.0180094, 1.0125180, 1.0233319, 1.0215106, 1.0198464, 1.0152214,
1.0180310,
+ 1.0150230, 0.9796248, 0.9835153, 0.9892937, 0.9734357, 0.9701524, 0.9773750,
0.9715888,
+ 0.9492091, 0.9047893, 0.8689214, 0.8969002, 0.8604160, 0.8650683, 0.8614112,
0.8602811,
+ 0.8938063, 0.8178848, 0.7873018, 0.7653331, 0.7517878, 0.7453435, 0.7279942,
0.7390122,
+ 0.9541119, 0.7347108, 0.6803388, 0.6805217, 0.6631648, 0.6455990, 0.6391322,
0.6154472)
> s1sd <- c(
+ 0.2558871, 0.2074931, 0.2019210, 0.1993456, 0.1956842, 0.1952699, 0.1954283,
0.1968623,
+ 0.4881609, 0.2834966, 0.2590079, 0.2311456, 0.2177313, 0.2055821, 0.1994074,
0.1880976,
+ 0.7429666, 0.4092906, 0.3167434, 0.2927429, 0.2527559, 0.2170875, 0.1969479,
0.1785792,
+ 1.0521131, 0.5281757, 0.4126196, 0.3583045, 0.2852777, 0.2379186, 0.1926662,
0.1721626,
+ 1.6983403, 0.6600759, 0.4703960, 0.3886030, 0.3205154, 0.2471658, 0.1904367,
0.1615976)
> s2mean <- c(
+ 0.8101896, 0.8270381, 0.8462382, 0.8353982, 0.8612579, 0.8739147, 0.9292625,
1.0165103,
+ 0.7661301, 0.8545842, 0.8382510, 0.7877157, 0.8323192, 0.8456999, 0.8914131,
0.9718037,
+ 0.9756622, 0.7809839, 0.7808141, 0.7915066, 0.7729451, 0.7805628, 0.7766208,
0.8476952,
+ 0.8473234, 0.6116651, 0.7240962, 0.6724986, 0.6674575, 0.6782188, 0.6906285,
0.7426497,
+ 0.9298414, 0.6441807, 0.6478278, 0.5738834, 0.5987168, 0.6156720, 0.6106082,
0.6276753)
> s2sd <- c(
+ 1.289713, 0.7079645, 0.5443964, 0.4557863, 0.3600675, 0.2905618, 0.2314770,
0.1919032,
+ 2.204718, 1.0352492, 0.7228366, 0.5769365, 0.4352699, 0.3227683, 0.2311042,
0.1866620,
+ 3.782651, 1.5552894, 0.9753109, 0.7328684, 0.5338099, 0.3525301, 0.2363640,
0.1756816,
+ 5.263356, 1.9738421, 1.2211028, 0.9583618, 0.6279809, 0.4022240, 0.2617814,
0.1741623,
+ 7.705801, 2.3052927, 1.3783467, 1.0642116, 0.6870980, 0.4257074, 0.2536335,
0.1686875)
> s3sd <- c(
+ 0.2672647, 0.2199326, 0.2075911, 0.1971651, 0.2061662, 0.2072476, 0.2130259,
0.1928112,
+ 0.4504687, 0.2883622, 0.2668242, 0.2480416, 0.2339451, 0.2200699, 0.2127237,
0.1898366,
+ 0.9104342, 0.4604099, 0.3874367, 0.3374590, 0.2976313, 0.2559326, 0.2473418,
0.1824161,
+ 1.2790115, 0.6850313, 0.5413326, 0.4623169, 0.4106435, 0.3290730, 0.2816922,
0.1615023,
+ 1.6172746, 0.9515100, 0.7337054, 0.5861010, 0.5215870, 0.4204276, 0.3748176,
0.1675852)
> s3mean <- c(
+ 1.019361, 1.032583, 1.038166, 1.026884, 1.034368, 1.037783, 1.0536807,
1.0260184,
+ 1.030941, 1.030202, 1.015011, 1.011430, 1.030089, 1.011929, 1.0315768,
0.9715992,
+ 1.094175, 1.020490, 1.017812, 1.006094, 1.024823, 1.013596, 0.9818869,
0.8581432,
+ 1.093622, 0.994017, 1.035233, 1.026690, 1.025811, 1.002776, 0.9660568,
0.7331220,
+ 1.070821, 1.107004, 1.020095, 1.019185, 1.067298, 1.017065, 0.9825400,
0.6296542)
> s4sd <- c(
+ 0.1977376, 0.1956897, 0.1990413, 0.1976083, 0.1988865, 0.2066219, 0.2133723,
0.1948147,
+ 0.1958121, 0.2031994, 0.1938832, 0.1984780, 0.2028359, 0.2143121, 0.2170887,
0.1887306,
+ 0.2021009, 0.2173261, 0.2122282, 0.2218343, 0.2180687, 0.2219821, 0.2329598,
0.1773513,
+ 0.2311124, 0.2233433, 0.2351271, 0.2221499, 0.2371908, 0.2401055, 0.2706758,
0.1754606,
+ 0.2394915, 0.2377767, 0.2393044, 0.2485314, 0.2487319, 0.2594596, 0.2858890,
0.1693721)
> s4mean <- c(
+ 1.025399, 1.035918, 1.039732, 1.039419, 1.034599, 1.043545, 1.0543309,
1.0208207,
+ 1.017997, 1.018147, 1.022816, 1.017625, 1.003781, 1.033450, 1.0349845,
0.9752180,
+ 1.022549, 1.011347, 1.019080, 1.030066, 1.036627, 1.021812, 0.9740666,
0.8501249,
+ 1.020098, 1.035326, 1.027258, 1.018737, 1.014026, 1.014551, 0.9599101,
0.7443075,
+ 1.033918, 1.022175, 1.019426, 1.022649, 1.026183, 1.008189, 0.9481097,
0.6284726)
>
>
>
> xvals <- c(5,10,15,20,30,50,100,200)
> dim(s1sd) <- c(8,5); dim(s2sd) <- c(8,5)
> dim(s3sd) <- c(8,5); dim(s4sd) <- c(8,5)
> dim(s1mean) <- c(8,5); dim(s2mean) <- c(8,5)
> dim(s3mean) <- c(8,5); dim(s4mean) <- c(8,5)
> s1sd <- t(s1sd); s2sd <- t(s2sd)
> s3sd <- t(s3sd); s4sd <- t(s4sd)
> s1mean <- t(s1mean); s2mean <- t(s2mean)
> s3mean <- t(s3mean); s4mean <- t(s4mean)
> colnames(s1sd) <- c(5,10,15,20,30,50,100,200)
> colnames(s2sd) <- c(5,10,15,20,30,50,100,200)
> colnames(s3sd) <- c(5,10,15,20,30,50,100,200)
> colnames(s4sd) <- c(5,10,15,20,30,50,100,200)
> colnames(s1mean) <- c(5,10,15,20,30,50,100,200)
> colnames(s2mean) <- c(5,10,15,20,30,50,100,200)
> colnames(s3mean) <- c(5,10,15,20,30,50,100,200)
> colnames(s4mean) <- c(5,10,15,20,30,50,100,200)
> rownames(s1sd) <- c(0,1,2,3,4); rownames(s2sd) <- c(0,1,2,3,4)
> rownames(s3sd) <- c(0,1,2,3,4); rownames(s4sd) <- c(0,1,2,3,4)
> rownames(s1mean) <- c(0,1,2,3,4); rownames(s2mean) <- c(0,1,2,3,4)
> rownames(s3mean) <- c(0,1,2,3,4); rownames(s4mean) <- c(0,1,2,3,4)
>
> # The means, not in paper
> par(mfrow=c(2,2))
> plot(log(xvals),s1mean[1,],ylim=c(.5,1.5),col=0,xaxt="n",
+ xlab="Number of juries (logged)",ylab="Mean",
+ cex.axis=1.3,cex.lab=1.3)
> axis(1,log(c(xvals)),colnames(s1mean),cex=1.3)
> axis(3,log(200/c(40,20,10,5,3,2,1)),c(40,20,10,5,3,2,1),cex=1.3)
> mtext("Number in each jury",3,line=2,cex=1)
> mtext("A", 3, line=3, cex=1.3)
> lines(spline(log(xvals),s1mean[1,]),lwd=.5)
> lines(spline(log(xvals),s1mean[2,]),lwd=1)
> lines(spline(log(xvals),s1mean[3,]),lwd=1.5)
> lines(spline(log(xvals),s1mean[4,]),lwd=2)
> lines(spline(log(xvals),s1mean[5,]),lwd=3)
> abline(h=1,lty=3)
> legend(3,1.57,c("Jury effect small","Jury effect large"),lwd=c(.5,3),cex=1)
>
>
> plot(log(xvals),s2mean[1,],ylim=c(-3,2),col=0,xaxt="n",
+ xlab="Number of juries (logged)",ylab="Mean",
+ cex.axis=1.3,cex.lab=1.3)
> axis(1,log(c(xvals)),colnames(s2mean),cex=1.3)
> axis(3,log(200/c(40,20,10,5,3,2,1)),c(40,20,10,5,3,2,1),cex=1.3)
> mtext("Number in each jury",3,line=2,cex=1)
> mtext("B", 3, line=3, cex=1.3)
> lines(spline(log(xvals),s2mean[1,]),lwd=.5)
> lines(spline(log(xvals),s2mean[2,]),lwd=1)
> lines(spline(log(xvals),s2mean[3,]),lwd=1.5)
> lines(spline(log(xvals),s2mean[4,]),lwd=2)
> lines(spline(log(xvals),s2mean[5,]),lwd=3)
> abline(h=1,lty=3)
> legend(3,-1.84,c("Jury effect small","Jury effect large"),lwd=c(.5,3),cex=1)
>
>
> plot(log(xvals),s3mean[1,],ylim=c(.5,1.5),col=0,xaxt="n",
+ xlab="Number of juries (logged)",ylab="Mean",
+ cex.axis=1.3,cex.lab=1.3)
> axis(1,log(c(xvals)),colnames(s3mean),cex=1.3)
> axis(3,log(200/c(40,20,10,5,3,2,1)),c(40,20,10,5,3,2,1),cex=1.3)
> mtext("Number in each jury",3,line=2,cex=1)
> mtext("C", 3, line=3, cex=1.3)
> lines(spline(log(xvals),s3mean[1,]),lwd=.5)
> lines(spline(log(xvals),s3mean[2,]),lwd=1)
> lines(spline(log(xvals),s3mean[3,]),lwd=1.5)
> lines(spline(log(xvals),s3mean[4,]),lwd=2)
> lines(spline(log(xvals),s3mean[5,]),lwd=3)
> abline(h=1,lty=3)
> legend(3,1.57,c("Jury effect small","Jury effect large"),lwd=c(.5,3),cex=1)
>
> plot(log(xvals),s4mean[1,],ylim=c(.5,1.5),col=0,xaxt="n",
+ xlab="Number of juries (logged)",ylab="Mean",
+ cex.axis=1.3,cex.lab=1.3)
> axis(1,log(c(xvals)),colnames(s4mean),cex=1.3)
> axis(3,log(200/c(40,20,10,5,3,2,1)),c(40,20,10,5,3,2,1),cex=1.3)
> mtext("Number in each jury",3,line=2,cex=1)
> mtext("D", 3, line=3, cex=1.3)
> lines(spline(log(xvals),s4mean[1,]),lwd=.5)
> lines(spline(log(xvals),s4mean[2,]),lwd=1)
> lines(spline(log(xvals),s4mean[3,]),lwd=1.5)
> lines(spline(log(xvals),s4mean[4,]),lwd=2)
> lines(spline(log(xvals),s4mean[5,]),lwd=3)
> abline(h=1,lty=3)
> legend(3,1.57,c("Jury effect small","Jury effect large"),lwd=c(.5,3),cex=1)
>
> >
> # The standard error, not in paper
(but curves basically the same)
> par(mfrow=c(2,2))
> plot(log(xvals),s1sd[1,],ylim=c(0,2),col=0,xaxt="n",
+ xlab="Number of juries (logged)",ylab="Standard error",
+ cex.axis=1.3,cex.lab=1.3,las=2)
> axis(1,log(c(xvals)),colnames(s1sd),cex=1.3)
> axis(3,log(200/c(40,20,10,5,3,2,1)),c(40,20,10,5,3,2,1),cex=1.3)
> mtext("Number in each jury",3,line=2,cex=1)
> mtext("A", 3, line=3, cex=1.3)
> lines(spline(log(xvals)[8:1],s1sd[1,8:1]),lwd=.5)
> lines(spline(log(xvals),s1sd[2,]),lwd=1)
> lines(spline(log(xvals),s1sd[3,]),lwd=1.5)
> lines(spline(log(xvals),s1sd[4,]),lwd=2)
> lines(spline(log(xvals),s1sd[5,]),lwd=3)
> legend(3,2.1,c("Jury similarity small","Jury similarity large"),lwd=c(.5,3),cex=1)
>
> plot(log(xvals),s2sd[1,],ylim=c(0,10),col=0,xaxt="n",
+ xlab="Number of juries (logged)",ylab="Standard error",
+ cex.axis=1.3,cex.lab=1.3)
> axis(1,log(c(xvals)),colnames(s2sd),cex=1.3)
> axis(3,log(200/c(40,20,10,5,3,2,1)),c(40,20,10,5,3,2,1),cex=1.3)
> mtext("Number in each jury",3,line=2,cex=1)
> mtext("B", 3, line=3, cex=1.3)
> lines(spline(log(xvals),s2sd[1,]),lwd=.5)
> lines(spline(log(xvals),s2sd[2,]),lwd=1)
> lines(spline(log(xvals),s2sd[3,]),lwd=1.5)
> lines(spline(log(xvals),s2sd[4,]),lwd=2)
> lines(spline(log(xvals),s2sd[5,]),lwd=3)
> legend(3,10.5,c("Jury similarity small","Jury similarity large"),lwd=c(.5,3),cex=1)
>
>
> plot(log(xvals),s3sd[1,],ylim=c(0,2),col=0,xaxt="n",
+ xlab="Number of juries (logged)",ylab="Standard error",
+ cex.axis=1.3,cex.lab=1.3)
> axis(1,log(c(xvals)),colnames(s3sd),cex=1.3)
> axis(3,log(200/c(40,20,10,5,3,2,1)),c(40,20,10,5,3,2,1),cex=1.3)
> mtext("Number in each jury",3,line=2,cex=1)
> mtext("C", 3, line=3, cex=1.3)
> lines(spline(log(xvals),s3sd[1,]),lwd=.5)
> lines(spline(log(xvals),s3sd[2,]),lwd=1)
> lines(spline(log(xvals),s3sd[3,]),lwd=1.5)
> lines(spline(log(xvals),s3sd[4,]),lwd=2)
> lines(spline(log(xvals),s3sd[5,]),lwd=3)
> legend(3,2.1,c("Jury similarity small","Jury similarity large"),lwd=c(.5,3),cex=1)
>
> plot(log(xvals),s4sd[1,],ylim=c(0,.5),col=0,xaxt="n",
+ xlab="Number of juries (logged)",ylab="Standard error",
+ cex.axis=1.3,cex.lab=1.3)
> axis(1,log(c(xvals)),colnames(s4sd),cex=1.3)
> axis(3,log(200/c(40,20,10,5,3,2,1)),c(40,20,10,5,3,2,1),cex=1.3)
> mtext("Number in each jury",3,line=2,cex=1)
> mtext("D", 3, line=3, cex=1.3)
> lines(spline(log(xvals),s4sd[1,]),lwd=.5)
> lines(spline(log(xvals),s4sd[2,]),lwd=1)
> lines(spline(log(xvals),s4sd[3,]),lwd=1.5)
> lines(spline(log(xvals),s4sd[4,]),lwd=2)
> lines(spline(log(xvals),s4sd[5,]),lwd=3)
> legend(3,.115,c("Jury similarity small","Jury similarity large"),lwd=c(.5,3),cex=1)
>
>
>
> s1nsd <- s1sd/s1sd[,8]
> s2nsd <- s2sd/s2sd[,8]
> s3nsd <- s3sd/s3sd[,8]
> s4nsd <- s4sd/s4sd[,8]
>
> # Figure 9
> par(mfrow=c(2,2))
> plot(200/xvals[8:1],s1nsd[1,8:1],col=0,
+ xlab="Number in each jury",ylab="Standard error multiplier",las=1,
+ cex.axis=1.3,cex.lab=1.3,ylim=c(1,10))
> mtext("A", 3, line=1, cex=1.3)
> abline(h=1)
> lines(200/xvals[8:1],s1nsd[1,8:1],lwd=.5)
> lines(200/xvals[8:1],s1nsd[2,8:1],lwd=1)
> lines(200/xvals[8:1],s1nsd[3,8:1],lwd=1.5)
> lines(200/xvals[8:1],s1nsd[4,8:1],lwd=2)
> lines(200/xvals[8:1],s1nsd[5,8:1],lwd=3)
> legend(-1,10.5,c("Jury similarity large","Jury similarity small"),lwd=c(3,.5),cex=1)
>
>
> plot(200/xvals[8:1],s2nsd[1,8:1],col=0,
+ xlab="Number in each jury",ylab="Standard error multiplier",las=1,
+ cex.axis=1.3,cex.lab=1.3,ylim=c(1,50))
> mtext("B", 3, line=1, cex=1.3)
> abline(h=1)
> lines(200/xvals[8:1],s2nsd[1,8:1],lwd=.5)
> lines(200/xvals[8:1],s2nsd[2,8:1],lwd=1)
> lines(200/xvals[8:1],s2nsd[3,8:1],lwd=1.5)
> lines(200/xvals[8:1],s2nsd[4,8:1],lwd=2)
> lines(200/xvals[8:1],s2nsd[5,8:1],lwd=3)
> legend(-1,52,c("Jury similarity large","Jury similarity small"),lwd=c(3,.5),cex=1)
>
> plot(200/xvals[8:1],s3nsd[1,8:1],col=0,
+ xlab="Number in each jury",ylab="Standard error multiplier",las=1,
+ cex.axis=1.3,cex.lab=1.3,ylim=c(1,10))
> mtext("C", 3, line=1, cex=1.3)
> abline(h=1)
> lines(200/xvals[8:1],s3nsd[1,8:1],lwd=.5)
> lines(200/xvals[8:1],s3nsd[2,8:1],lwd=1)
> lines(200/xvals[8:1],s3nsd[3,8:1],lwd=1.5)
> lines(200/xvals[8:1],s3nsd[4,8:1],lwd=2)
> lines(200/xvals[8:1],s3nsd[5,8:1],lwd=3)
> legend(-1,10.5,c("Jury similarity large","Jury similarity small"),lwd=c(3,.5),cex=1)
>
>
> plot(200/xvals[8:1],s4nsd[1,8:1],col=0,
+ xlab="Number in each jury",ylab="Standard error multiplier",las=1,
+ cex.axis=1.3,cex.lab=1.3,ylim=c(1,2))
> mtext("D", 3, line=1, cex=1.3)
> abline(h=1)
> lines(200/xvals[8:1],s4nsd[1,8:1],lwd=.5)
> lines(200/xvals[8:1],s4nsd[2,8:1],lwd=1)
> lines(200/xvals[8:1],s4nsd[3,8:1],lwd=1.5)
> lines(200/xvals[8:1],s4nsd[4,8:1],lwd=2)
> lines(200/xvals[8:1],s4nsd[5,8:1],lwd=3)
> legend(12,2.05,c("Jury similarity large","Jury similarity small"),lwd=c(3,.5),cex=1)
>
>
>
> # Figure 10
> par(mfrow=c(2,2))
> plot(200/xvals[8:1],s1nsd[1,8:1],col=0,
+ xlab="Number in each jury",ylab="Standard error multiplier",las=1,
+ cex.axis=1.3,cex.lab=1.3,ylim=c(1,3),xlim=c(0,10))
> mtext("A", 3, line=1, cex=1.3)
> abline(h=1)
> lines(200/xvals[8:1],s1nsd[1,8:1],lwd=.5)
> lines(200/xvals[8:1],s1nsd[2,8:1],lwd=1)
> lines(200/xvals[8:1],s1nsd[3,8:1],lwd=1.5)
> lines(200/xvals[8:1],s1nsd[4,8:1],lwd=2)
> lines(200/xvals[8:1],s1nsd[5,8:1],lwd=3)
> legend(-1,3.1,c("Jury similarity large","Jury similarity small"),lwd=c(3,.5),cex=1)
>
>
> plot(200/xvals[8:1],s2nsd[1,8:1],col=0,
+ xlab="Number in each jury",ylab="Standard error multiplier",las=1,
+ cex.axis=1.3,cex.lab=1.3,ylim=c(1,8),xlim=c(0,10))
> mtext("B", 3, line=1, cex=1.3)
> abline(h=1)
> lines(200/xvals[8:1],s2nsd[1,8:1],lwd=.5)
> lines(200/xvals[8:1],s2nsd[2,8:1],lwd=1)
> lines(200/xvals[8:1],s2nsd[3,8:1],lwd=1.5)
> lines(200/xvals[8:1],s2nsd[4,8:1],lwd=2)
> lines(200/xvals[8:1],s2nsd[5,8:1],lwd=3)
> legend(-1,8.26,c("Jury similarity large","Jury similarity small"),lwd=c(3,.5),cex=1)
>
> plot(200/xvals[8:1],s3nsd[1,8:1],col=0,
+ xlab="Number in each jury",ylab="Standard error multiplier",las=1,
+ cex.axis=1.3,cex.lab=1.3,ylim=c(1,4),xlim=c(0,10))
> mtext("C", 3, line=1, cex=1.3)
> abline(h=1)
> lines(200/xvals[8:1],s3nsd[1,8:1],lwd=.5)
> lines(200/xvals[8:1],s3nsd[2,8:1],lwd=1)
> lines(200/xvals[8:1],s3nsd[3,8:1],lwd=1.5)
> lines(200/xvals[8:1],s3nsd[4,8:1],lwd=2)
> lines(200/xvals[8:1],s3nsd[5,8:1],lwd=3)
> legend(-1,4.14,c("Jury similarity large","Jury similarity small"),lwd=c(3,.5),cex=1)
>
> plot(200/xvals[8:1],s4nsd[1,8:1],ylim=c(1,2),col=0,xlim=c(0,10),
+ xlab="Number in each jury",ylab="Standard error multiplier",las=1,
+ cex.axis=1.3,cex.lab=1.3)
> mtext("D", 3, line=1, cex=1.3)
> abline(h=1)
> lines(200/xvals[8:1],s4nsd[1,8:1],lwd=.5)
> lines(200/xvals[8:1],s4nsd[2,8:1],lwd=1)
> lines(200/xvals[8:1],s4nsd[3,8:1],lwd=1.5)
> lines(200/xvals[8:1],s4nsd[4,8:1],lwd=2)
> lines(200/xvals[8:1],s4nsd[5,8:1],lwd=3)
> legend(2.7,2.05,c("Jury similarity large","Jury similarity small"),lwd=c(3,.5),cex=1)
>
>
>
>
