### R code from vignette source 'lecture6_2015.Rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: readpolpar ################################################### library(mokken) library(foreign) dat <- read.dta("http://www.quantoid.net/files/essex/anes_2012_polpar.dta") pp <- dat[,grep("polpar", colnames(dat))] H <- coefH(pp) H$Hi H$H ################################################### ### code chunk number 2: checkni ################################################### summary(check.iio(pp)) ################################################### ### code chunk number 3: checkni2 ################################################### pp2 <- pp[,-3] summary(check.iio(pp2)) ################################################### ### code chunk number 4: testscores ################################################### scores <- rowSums(pp2) table(scores) ################################################### ### code chunk number 5: readwvs ################################################### u <- url("http://www.quantoid.net/files/essex/wvs6.rda") load(u) close(u) sub.wvs <- wvs6[which(wvs6$Country == "Algeria"),] ################################################### ### code chunk number 6: ocsc ################################################### u <- url("http://www.quantoid.net/files/essex/supreme_court_2010-2013.rda") load(u) close(u) library(pscl) library(oc) sc.rc <- rollcall(t(scv[,-10]), desc = "Supreme Court 2010-2013", legis.names =colnames(scv)[-10]) sc.oc <- oc(sc.rc, dims=2, polarity=c(1,2)) ################################################### ### code chunk number 7: sumscoc ################################################### sc.oc$fits ################################################### ### code chunk number 8: judges ################################################### sc.oc$legislators ################################################### ### code chunk number 9: plotoc ################################################### plot(sc.oc) ################################################### ### code chunk number 10: readehrc ################################################### u <- url("http://www.quantoid.net/files/essex/echr.rda") load(u) close(u) ################################################### ### code chunk number 11: pairdifffun ################################################### pair.diffs <- function(X){ if(!is.matrix(X))X <- as.matrix(X) d <- combn(ncol(X), 2) D <- matrix(0, ncol=ncol(d), nrow=ncol(X)) D[cbind(d[1,], 1:ncol(D))] <- 1 D[cbind(d[2,], 1:ncol(D))] <- -1 Xtmp <- X Xtmp[which(is.na(Xtmp), arr.ind=T)] <- 10000 diffs <- Xtmp %*% D diffs[which(abs(diffs) > 100, arr.ind = T)] <- NA diffs[which(diffs == 0, arr.ind=T)] <- NA diffs <- diffs > 0 colnames(diffs) <- paste(colnames(X)[d[1,]], colnames(X)[d[2,]], sep="-") diffs } ################################################### ### code chunk number 12: octherm ################################################### u <- url("http://www.quantoid.net/files/essex/candidatetherms2008.rda") load(u) close(u) therms <- as.matrix(candidatetherms2008) pdt <- pair.diffs(therms) rownames(pdt) <- 1:nrow(pdt) pairs.rc <- rollcall(pdt, vote.names=colnames(pdt)) pair.oc <- oc(pairs.rc, polarity=c(3,3), dims=2, minvotes=20) ################################################### ### code chunk number 13: sumpairoc ################################################### summary(pair.oc) pair.oc$fits ################################################### ### code chunk number 14: plotpoc ################################################### plot(pair.oc) ################################################### ### code chunk number 15: findanglesfun ################################################### findOCangles <- function(obj, data){ oc1 <- obj$legislators[,7] oc2 <- obj$legislators[,8] PRE <- obj$rollcalls[,5] N1 <- obj$rollcalls[,6] N2 <- obj$rollcalls[,7] ws <- obj$rollcalls[,8] xws <- ws * N1 yws <- ws * N2 C1 <- N2 C2 <- -N1 for (i in 1:nrow(obj$rollcalls)){ if (C1[i] < 0 & !is.na(C2[i])) C2[i] <- -C2[i] if (C1[i] < 0 & !is.na(C1[i])) C1[i] <- -C1[i] } theta <- atan2(C2,C1) theta4 <- theta * (180/pi) data.frame(colnames(data$votes), theta4) } ################################################### ### code chunk number 16: showangles ################################################### pairs.angles <- findOCangles(pair.oc, pairs.rc) pairs.angles[grep("mccain", pairs.angles[,1]), ] ################################################### ### code chunk number 17: findstims ################################################### plotOCstims <- function(obj, data, who="",...){ oc1 <- obj$legislators[,7] oc2 <- obj$legislators[,8] PRE <- obj$rollcalls[,5] N1 <- obj$rollcalls[,6] N2 <- obj$rollcalls[,7] ws <- obj$rollcalls[,8] xws <- ws * N1 yws <- ws * N2 C1 <- N2 C2 <- -N1 for (i in 1:nrow(obj$rollcalls)){ if (C1[i] < 0 & !is.na(C2[i])) C2[i] <- -C2[i] if (C1[i] < 0 & !is.na(C1[i])) C1[i] <- -C1[i] } plot(oc1, oc2, ..., xlab="First Dimension (Liberal - Conservative)", ylab="Second Dimension", xlim=c(-1,1),ylim=c(-1,1), asp=1, type="n") # inds <- grep(who, colnames(data$votes)) for (i in inds){ vote <- as.integer(data$votes[,i]) polarity <- oc1*N1[i] + oc2*N2[i] - ws[i] errors1 <- vote==1 & polarity >= 0 errors2 <- vote==0 & polarity <= 0 errors3 <- vote==1 & polarity <= 0 errors4 <- vote==0 & polarity >= 0 kerrors12 <- sum(errors1==1,na.rm=T) + sum(errors2==1,na.rm=T) kerrors34 <- sum(errors3==1,na.rm=T) + sum(errors4==1,na.rm=T) if(kerrors12 < kerrors34){ xwslow <- (ws[i] - 0.1) * N1[i] ywslow <- (ws[i] - 0.1) * N2[i] } if(kerrors12 >= kerrors34){ xwslow <- (ws[i] + 0.1) * N1[i] ywslow <- (ws[i] + 0.1) * N2[i] } segments(xws[i]+N2[i], yws[i]-N1[i], xws[i]-N2[i], yws[i]+N1[i], lwd=2, col="black") arrows(xws[i]+N2[i], yws[i]-N1[i], xwslow+N2[i], ywslow-N1[i], length=0.05, lwd=2, col="gray33") arrows(xws[i]-N2[i], yws[i]+N1[i], xwslow-N2[i], ywslow+N1[i], length=0.05, lwd=2, col="gray33") } } ################################################### ### code chunk number 18: plotstims ################################################### plotOCstims(pair.oc, pairs.rc, "mccain") ################################################### ### code chunk number 19: readtherm ################################################### u <- url("http://www.quantoid.net/files/essex/thermometers2004.rda") load(u) close(u)