# # 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) # # &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& # FUNCTION CALLED BY OPTIM(...) Below # # *** Calculate Log-Likelihood -- In this case simple SSE *** # ************************************************************** fr <- function(zz){ xx <- zz # i <- 1 sse <- 0.0 while (i <= np) { ii <- 1 while (ii <= i) { j <- 1 dist_i_ii <- 0.0 while (j <= ndim) { # # Calculate distance between points # dist_i_ii <- dist_i_ii + (xx[ndim*(i-1)+j]-xx[ndim*(ii-1)+j])*(xx[ndim*(i-1)+j]-xx[ndim*(ii-1)+j]) j <- j + 1 } sse <- sse + (sqrt(TT[i,ii]) - sqrt(dist_i_ii))*(sqrt(TT[i,ii]) - sqrt(dist_i_ii)) ii <- ii + 1 } i <- i + 1 } return(sse) } # # &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& # # 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) # # 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 ssesvd should be zero! # # i <- 0 j <- 0 sse <- 0 while (i < np) { i <- i + 1 j <- 0 while (j < np) { j <- j + 1 sse <- sse + (TTT[i,j] - XX[i,j])**2 } } # # 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 } # # # ******************************************* # DO MAXIMUM LIKELIHOOD MAXIMIZATION HERE # ******************************************* # model <- optim(zz,fr,hessian=TRUE,control=list(maxit=50000)) # # Log-Likelihood (inverse -- optim minimizes!!) # logLmax <- model$value # # Parameter Estimates # rhomax <- model$par # # convergence an integer code. # 0 indicates successful convergence. # Error codes are # 1 indicates that the iteration limit maxit had been reached. # 10 indicates degeneracy of the Nelder–Mead simplex. # 51 indicates a warning from the "L-BFGS-B" method; see component message for further details. # 52 indicates an error from the "L-BFGS-B" method; see component message for further details. # nconverge <- model$convergence # # counts # A two-element integer vector giving the number of calls to # fn and gr respectively. ncounts <- model$counts # xhessian <- model$hessian # # Check Rank of Hessian Matrix # evv <- eigen(xhessian) # # evv$values -- shows the eigenvalues # # Put Coordinates into matrix for plot purposes # xmetric <- rep(0,np*ndim) dim(xmetric) <- c(np,ndim) # i <- 1 while (i <= np) { j <- 1 while (j <= ndim) { xmetric[i,j] <- rhomax[ndim*(i-1)+j] j <- j +1 } i <- i + 1 } # 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)