# # montecarlo_smacof_MLSMU6.r # # ANALYZING SPATIAL MODELS OF CHOICE AND JUDGMENT WITH R # Dave Armstrong, Ryan Bakker, Royce Carroll, Christopher Hare, Keith T. Poole, and Howard Rosenthal # http://voteview.com/asmcjr.asp # # # Remove all objects just to be safe. rm(list=ls(all=TRUE)) # ngroups <- 12 nlegis <- 100 dims <- 2 ntrials <- 1000 error <- 0.1 # This is the standard deviation on the noise term (drawn from a normal distribution with mean 0) added to distances to compute ratings nmissing <- 60 # Number of missing entries to random insert in the ratings matrix # smacof.legislator.correlations <- matrix(NA, nrow=ntrials, ncol=dims) smacof.interest.group.correlations <- matrix(NA, nrow=ntrials, ncol=dims) MLSMU6.legislator.correlations <- matrix(NA, nrow=ntrials, ncol=dims) MLSMU6.interest.group.correlations <- matrix(NA, nrow=ntrials, ncol=dims) # for(q in 1:ntrials){ # # interest.groups <- matrix(NA, nrow=ngroups, ncol=dims) legislators <- matrix(NA, nrow=nlegis, ncol=dims) # for (k in 1:dims){ interest.groups[,k] <- runif(ngroups, -1, 1) legislators[,k] <- runif(nlegis, -1, 1) } # distances <- matrix(NA, nrow=nlegis, ncol=ngroups) # for (i in 1:nlegis){ for (j in 1:ngroups){ distances[i,j] <- sum(abs(interest.groups[j,1:dims] - legislators[i,1:dims])) }} # ratings <- matrix(NA, nrow=nlegis, ncol=ngroups) # for (i in 1:nlegis){ for (j in 1:ngroups){ ratings[i,j] <- 100 - 50*(distances[i,j]+rnorm(1,0,error)) }} # ratings[ratings > 100] <- 100 ratings[ratings < 0] <- 0 # # # T <- ratings T <- (100-T)/50 # # Insert random missing values: i <- sample(1:nlegis, nmissing, replace=TRUE) j <- sample(1:ngroups, nmissing, replace=TRUE) for (v in 1:nmissing){ T[i[v],j[v]] <- NA } # # SMACOF: # library(smacof) # weightmat <- T weightmat[!is.na(T)] <- 1 weightmat[is.na(T)] <- 0 T[is.na(T)] <- mean(T,na.rm=TRUE) # result <- smacofRect(T, ndim=2, weightmat=weightmat) smacof.legislators <- result$conf.row smacof.interest.groups <- result$conf.col # for (k in 1:dims){ if (abs(cor(legislators,smacof.legislators)[1,1]) > abs(cor(legislators,smacof.legislators)[1,2])) smacof.legislator.correlations[q,] <- c(abs(cor(legislators,smacof.legislators)[1,1]), abs(cor(legislators,smacof.legislators)[2,2])) if (abs(cor(legislators,smacof.legislators)[1,1]) < abs(cor(legislators,smacof.legislators)[1,2])) smacof.legislator.correlations[q,] <- c(abs(cor(legislators,smacof.legislators)[1,2]), abs(cor(legislators,smacof.legislators)[2,1])) # if (abs(cor(interest.groups,smacof.interest.groups)[1,1]) > abs(cor(interest.groups,smacof.interest.groups)[1,2])) smacof.interest.group.correlations[q,] <- c(abs(cor(interest.groups,smacof.interest.groups)[1,1]), abs(cor(interest.groups,smacof.interest.groups)[2,2])) if (abs(cor(interest.groups,smacof.interest.groups)[1,1]) < abs(cor(interest.groups,smacof.interest.groups)[1,2])) smacof.interest.group.correlations[q,] <- c(abs(cor(interest.groups,smacof.interest.groups)[1,2]), abs(cor(interest.groups,smacof.interest.groups)[2,1])) } # # # MLSMU6: # doubleCenterRect <- function(T){ p <- dim(T)[1] n <- dim(T)[2] -(T-matrix(apply(T,1,mean),nrow=p,ncol=n) - t(matrix(apply(T,2,mean),nrow=n,ncol=p)) + mean(T))/2 } # np <- nrow(T) nq <- ncol(T) xrow <- sapply(1:np,function(i) length(rep(1,nq)[!is.na(T[i,])])) xcol <- sapply(1:nq,function(j) length((1:np)[!is.na(T[,j])])) # TT <- T TT[is.na(TT)] <- mean(T,na.rm=TRUE) TTSQ <- T*T TTSQ[is.na(T)] <- (mean(T,na.rm=TRUE))**2 TEIGHT <- as.numeric(TT) dim(TEIGHT) <- c(np,nq) TTSQDC <- doubleCenterRect(TTSQ) # xsvd <- svd(TTSQDC) zz <- xsvd$v[,1:dims] xx <- rep(0,np*dims) dim(xx) <- c(np,dims) for (i in 1:dims){ xx[,i] <- xsvd$u[,i]*sqrt(xsvd$d[i]) } # sumaj <- function(i){ jjj <- 1:nq j <- jjj[!is.na(T[i,])] s <- (xx[i,1]-zz[j,1])^2 + (xx[i,2]-zz[j,2])^2 ### Note this needs to be expanded if estimating more than two dimensions s=sqrt(s) sx=TEIGHT[i,j] sum((s-sx)^2) } # sumvector <- sapply(1:np,sumaj) suma <- sum(sumvector) xxx <- xx zzz <- zz kp=0 ktp=0 sumalast <- 100000.00 SAVEsumalast <- sumalast # for(loop in 1:50){ xxx[,1:dims] <- 0 zzz[,1:dims] <- 0 kp=kp+1 ktp=ktp+1 # for(i in 1:np){ for(j in 1:nq){ if(!is.na(T[i,j])){ s=0 for(k in 1:dims) s <- s+(xx[i,k]-zz[j,k])^2 xc <- ifelse(s==0, 1.0, TEIGHT[i,j]/sqrt(s)) for(k in 1:dims) zzz[j,k]=zzz[j,k]+(xx[i,k]-xc*(xx[i,k]-zz[j,k]))/xcol[j] } } } for(k in 1:dims){ for(i in 1:np){ sw <- 0.0 for(j in 1:nq){ if(!is.na(T[i,j])){ s=0 for(kk in 1:dims) s <- s+(xx[i,kk]-zzz[j,kk])^2 xc <- ifelse(s==0, 1.0, TEIGHT[i,j]/sqrt(s)) xxx[i,k]=xxx[i,k]+(zzz[j,k]-xc*(zzz[j,k]-xx[i,k])) } } xxx[i,k]=xxx[i,k]/xrow[i] } } xx <- xxx zz <- zzz sumvector <- sapply(1:np,sumaj) suma <- sum(sumvector) # print(c(loop,suma)) done=((sumalast-suma)/suma)<0.0005 sumalast=suma if(done) break } # # MLSMU6.legislators <- xx MLSMU6.interest.groups <- zz # for (k in 1:dims){ if (abs(cor(legislators,MLSMU6.legislators)[1,1]) > abs(cor(legislators,MLSMU6.legislators)[1,2])) MLSMU6.legislator.correlations[q,] <- c(abs(cor(legislators,MLSMU6.legislators)[1,1]), abs(cor(legislators,MLSMU6.legislators)[2,2])) if (abs(cor(legislators,MLSMU6.legislators)[1,1]) < abs(cor(legislators,MLSMU6.legislators)[1,2])) MLSMU6.legislator.correlations[q,] <- c(abs(cor(legislators,MLSMU6.legislators)[1,2]), abs(cor(legislators,MLSMU6.legislators)[2,1])) # if (abs(cor(interest.groups,MLSMU6.interest.groups)[1,1]) > abs(cor(interest.groups,MLSMU6.interest.groups)[1,2])) MLSMU6.interest.group.correlations[q,] <- c(abs(cor(interest.groups,MLSMU6.interest.groups)[1,1]), abs(cor(interest.groups,MLSMU6.interest.groups)[2,2])) if (abs(cor(interest.groups,MLSMU6.interest.groups)[1,1]) < abs(cor(interest.groups,MLSMU6.interest.groups)[1,2])) MLSMU6.interest.group.correlations[q,] <- c(abs(cor(interest.groups,MLSMU6.interest.groups)[1,2]), abs(cor(interest.groups,MLSMU6.interest.groups)[2,1])) } # # # # END OF FINAL Q LOOP: } # # windows() par(mfrow=c(1,2)) plot(smacof.legislator.correlations[,1], MLSMU6.legislator.correlations[,1], xlim=c(0,1), ylim=c(0,1), main="Legislator First Dimension Correlations\n5% Missing Data", xlab="SMACOF Correlations", ylab="MLSMU6 Correlations") abline(a=0, b=1, lwd=2) text(0.7, 0.95, table(smacof.legislator.correlations[,1] > MLSMU6.legislator.correlations[,1])[1]) text(0.95, 0.7, table(smacof.legislator.correlations[,1] > MLSMU6.legislator.correlations[,1])[2]) # plot(smacof.legislator.correlations[,2], MLSMU6.legislator.correlations[,2], xlim=c(0,1), ylim=c(0,1), main="Legislator Second Dimension Correlations\n5% Missing Data", xlab="SMACOF Correlations", ylab="MLSMU6 Correlations") abline(a=0, b=1, lwd=2) text(0.7, 0.95, table(smacof.legislator.correlations[,2] > MLSMU6.legislator.correlations[,2])[1]) text(0.95, 0.7, table(smacof.legislator.correlations[,2] > MLSMU6.legislator.correlations[,2])[2]) # # windows() par(mfrow=c(1,2)) plot(smacof.interest.group.correlations[,1], MLSMU6.interest.group.correlations[,1], xlim=c(0,1), ylim=c(0,1), main="Interest Group First Dimension Correlations\n5% Missing Data", xlab="SMACOF Correlations", ylab="MLSMU6 Correlations") abline(a=0, b=1, lwd=2) text(0.7, 0.95, table(smacof.interest.group.correlations[,1] > MLSMU6.interest.group.correlations[,1])[1]) text(0.95, 0.7, table(smacof.interest.group.correlations[,1] > MLSMU6.interest.group.correlations[,1])[2]) # plot(smacof.interest.group.correlations[,2], MLSMU6.interest.group.correlations[,2], xlim=c(0,1), ylim=c(0,1), main="Interest Group Second Dimension Correlations\n5% Missing Data", xlab="SMACOF Correlations", ylab="MLSMU6 Correlations") abline(a=0, b=1, lwd=2) text(0.7, 0.95, table(smacof.interest.group.correlations[,2] > MLSMU6.interest.group.correlations[,2])[1]) text(0.95, 0.7, table(smacof.interest.group.correlations[,2] > MLSMU6.interest.group.correlations[,2])[2]) # mean(smacof.legislator.correlations) mean(MLSMU6.legislator.correlations) mean(smacof.interest.group.correlations) mean(MLSMU6.interest.group.correlations) #