# # smacof_metric_nations.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)) # library(smacof) # load(url("ftp://voteview.com/wf1/nations.Rda")) # # Data Must Be Transformed to Squared Distances Below # # BRAZIL 9.00 4.83 5.28 3.44 4.72 4.50 3.83 3.50 2.39 3.06 5.39 3.17 # CONGO 4.83 9.00 4.56 5.00 4.00 4.83 3.33 3.39 4.00 3.39 2.39 3.50 # CUBA 5.28 4.56 9.00 5.17 4.11 4.00 3.61 2.94 5.50 5.44 3.17 5.11 # EGYPT 3.44 5.00 5.17 9.00 4.78 5.83 4.67 3.83 4.39 4.39 3.33 4.28 # FRANCE 4.72 4.00 4.11 4.78 9.00 3.44 4.00 4.22 3.67 5.06 5.94 4.72 # INDIA 4.50 4.83 4.00 5.83 3.44 9.00 4.11 4.50 4.11 4.50 4.28 4.00 # ISRAEL 3.83 3.33 3.61 4.67 4.00 4.11 9.00 4.83 3.00 4.17 5.94 4.44 # JAPAN 3.50 3.39 2.94 3.83 4.22 4.50 4.83 9.00 4.17 4.61 6.06 4.28 # CHINA 2.39 4.00 5.50 4.39 3.67 4.11 3.00 4.17 9.00 5.72 2.56 5.06 # USSR 3.06 3.39 5.44 4.39 5.06 4.50 4.17 4.61 5.72 9.00 5.00 6.67 # USA 5.39 2.39 3.17 3.33 5.94 4.28 5.94 6.06 2.56 5.00 9.00 3.56 # YUGOSLAVIA 3.17 3.50 5.11 4.28 4.72 4.00 4.44 4.28 5.06 6.67 3.56 9.00 # d <- (9-nations)^2 # smacof_metric_result <- smacofSym(delta=d, ndim=2, weightmat=NULL, type="ratio", itmax = 1000, eps=0.000001) # conf <- smacof_metric_result$conf print(smacof_metric_result$niter) print(smacof_metric_result$stress) # plot(conf[,1],conf[,2], main="SMACOF (Metric) Solution", xlab="", ylab="", xlim=c(-1,1), ylim=c(-1,1), asp=1, type="n") points(conf[,1], conf[,2], pch=16, col="gray") text(-0.9, 0.9, paste("Stress = ", round(smacof_metric_result$stress,3)), font=2) text(conf[,1], conf[,2], rownames(nations), pos=rep(4,nrow(nations)), offset=00.25, col="black", font=2) abline(a=0, b=-3/5, lwd=2, lty=2) abline(a=0, b=5/3, lwd=2, lty=2) # ndim <- 5 result <- vector("list", ndim) for (i in 1:ndim){ result[[i]] <- smacofSym(delta=d, ndim=i, type="ratio") } stress <- sapply(result, function(x)x$stress) # windows() plot(1:ndim,stress[1:ndim], main="Stress Values, 1-5 Dimensions", xlab="Dimensions", ylab="Stress", xlim=c(1,5), ylim=c(0,0.5), pch=16, lwd=2, font=2, type="o") #