# # metric_mds_nations.r -- Does Metric MDS # # Wish 1971 nations similarities data # Title: Perceived similarity of nations # # Description: Wish (1971) asked 18 students to rate the global # similarity of different pairs of nations such as 'France and China' on # a 9-point rating scale ranging from `1=very different' to `9=very # similar'. The table shows the mean similarity ratings. # # 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 # rm(list=ls(all=TRUE)) library(MASS) library(stats) # &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& # # Read Nations Similarities Data # T <- matrix(scan("C:/uga_course_homework_6_2015/nations.txt",0),ncol=12,byrow=TRUE) # # np <- length(T[,1]) ndim <- 2 nparam <- np*ndim # names <- NULL names[1] <- "Brazil" names[2] <- "Congo" names[3] <- "Cuba" names[4] <- "Egypt" names[5] <- "France" names[6] <- "India" names[7] <- "Israel" names[8] <- "Japan" names[9] <- "China" names[10] <- "USSR" names[11] <- "USA" names[12] <- "Yugo" # # # pos -- a position specifier for the text. Values of 1, 2, 3 and 4, # respectively indicate positions below, to the left of, above and # to the right of the specified coordinates # namepos <- rep(4,np) #namepos[1] <- 4 # Brazil #namepos[2] <- 4 # Congo #namepos[3] <- 4 # Cuba #namepos[4] <- 4 # Egypt #namepos[5] <- 4 # France #namepos[6] <- 4 # India #namepos[7] <- 4 # Israel #namepos[8] <- 4 # Japan #namepos[9] <- 4 # China #namepos[10] <- 4 # USSR #namepos[11] <- 4 # USA #namepos[12] <- 4 # Yugoslavia # # Set Parameter Values # zz <- rep(0,np*ndim) dim(zz) <- c(np*ndim,1) xx <- rep(0,np*ndim) dim(xx) <- c(np*ndim,1) xmetric <- rep(0,np*ndim) dim(xmetric) <- c(np,ndim) # # Set up for Double-Centering # TT <- rep(0,np*np) dim(TT) <- c(np,np) TTT <- rep(0,np*np) dim(TTT) <- c(np,np) # xrow <- NULL xcol <- NULL xcount <- NULL matrixmean <- 0 # # Transform the Matrix # i <- 0 while (i < np) { i <- i + 1 xcount[i] <- i j <- 0 while (j < np) { j <- j + 1 # # This is the Normal Transformation # TT[i,j] <- ((9 - T[i,j]))**2 # # But the R subroutine wants simple Distances! # TTT[i,j] <- (9 - T[i,j]) # } } # # Call Double Center Routine From R Program # cmdscale(....) in stats library # The Input data are DISTANCES!!! Not Squared Distances!! # Note that the R routine does not divide # by -2.0 # # dcenter <- cmdscale(TTT,ndim, eig=FALSE,add=FALSE,x.ret=TRUE) # # returns double-center matrix as dcenter$x if x.ret=TRUE # # Do the Division by -2.0 Here # TTT <- (dcenter$x)/(-2.0) # # Perform Eigenvalue-Eigenvector Decomposition of Double-Centered Matrix # # ev$values are the eigenvectors # ev$vectors are the eigenvectors # ev <- eigen(TTT) # The Two Lines Below Put the Eigenvalues in a # Diagonal Matrix -- The first one creates an # identity matrix and the second command puts # the singular values on the diagonal # Lambda <- diag(np) diag(Lambda) <- ev$val # # Compute U*LAMBDA*U' for check below # # XX <- ev$vec %*% Lambda %*% t(ev$vec) # # # Compute Fit of decomposition -- This is just the sum of squared # error -- Note that xfit should be zero! # xfit <- sum(abs(TTT-XX)) # # Use Eigenvectors for Starts # i <- 1 while (i <= np) { j <- 1 while (j <= ndim) { zz[ndim*(i-1)+j]=ev$vector[i,j] j <- j +1 } i <- i + 1 } # # # SETUP FOR METRIC MDS SCALING USED BY POOLE (1984) PSYCHOMETRIKA # # zhat <- NULL dstar <- rep(0,np*np) dim(dstar) <- c(np,np) dhat <- rep(0,np*np) dim(dhat) <- c(np,np) dstar <- sqrt(TT) SSESSE <- NULL # # Calculate SSE # i <- 1 while (i <= np){ j <- 1 while (j <= np){ k <- 1 while (k <= ndim){ dhat[i,j] <- dhat[i,j] + (zz[ndim*(i-1)+k]-zz[ndim*(j-1)+k])*(zz[ndim*(i-1)+k]-zz[ndim*(j-1)+k]) k <- k + 1 } dhat[i,j] <- sqrt(dhat[i,j]) j <- j + 1 } i <- i + 1 } SSE1 <- sum((dstar-dhat)**2) SSESSE[1] <- SSE1 # # GLOBAL LOOP -- RUNS THE ALGORITHM USING THE BASIC GEOMETRY OF THE FIRST DERIVATIVES # iiii <- 1 while (iiii <= 100){ # i <- 1 while (i <= np){ j <- 1 k <- 1 while (k <= ndim){ zhat[k] <- 0.0 k <- k + 1 } while (j <= np){ dhat[i,j] <- 0 k <- 1 while (k <= ndim){ dhat[i,j] <- dhat[i,j] + (zz[ndim*(i-1)+k]-zz[ndim*(j-1)+k])*(zz[ndim*(i-1)+k]-zz[ndim*(j-1)+k]) k <- k + 1 } if (j != i){ dhat[i,j] <- sqrt(dhat[i,j]) dratio <- dstar[i,j]/dhat[i,j] k <- 1 while (k <= ndim){ zhat[k] <- zhat[k] + zz[ndim*(j-1)+k] + dratio*(zz[ndim*(i-1)+k]-zz[ndim*(j-1)+k]) k <- k + 1 } } j <- j + 1 } k <- 1 while ( k <= ndim){ zz[ndim*(i-1)+k] <- zhat[k]/(np-1) k <- k + 1 } i <- i + 1 } # # Calculate SSE # i <- 1 while (i <= np){ j <- 1 while (j <= np){ k <- 1 dhat[i,j] <- 0.0 while (k <= ndim){ dhat[i,j] <- dhat[i,j] + (zz[ndim*(i-1)+k]-zz[ndim*(j-1)+k])*(zz[ndim*(i-1)+k]-zz[ndim*(j-1)+k]) xmetric[i,k] <- zz[ndim*(i-1)+k] k <- k + 1 } dhat[i,j] <- sqrt(dhat[i,j]) j <- j + 1 } i <- i + 1 } # STORE SSE THIS ITERATION SSE2 <- sum((dstar-dhat)**2) SSESSE[iiii+1] <- SSE2 # iiii <- iiii + 1 } # END OF LOOP # xmax <- max(xmetric) xmin <- min(xmetric) # plot(xmetric[,1],xmetric[,2],type="n",asp=1, main="", xlab="", ylab="", xlim=c(xmin,xmax),ylim=c(xmin,xmax),font=2) points(xmetric[,1],xmetric[,2],pch=16,col="red") axis(1,font=2) axis(2,font=2,cex=1.2) # # Main title mtext("12 Nations Data",side=3,line=1.50,cex=1.2,font=2) # x-axis title #mtext("Liberal-Conservative",side=1,line=2.75,cex=1.2) # y-axis title #mtext("Region (Civil Rights)",side=2,line=2.5,cex=1.2) # text(xmetric[,1],xmetric[,2],names,pos=namepos,offset=00.00,col="blue",font=2) #