# http://gastonsanchez.com/blog/how-to/2013/01/23/MDS-in-R.html rm(list=ls()) # installing and loading packages #install.packages(c("vegan", "ecodist", "labdsv", "ape", "ade4", "smacof")) # Data eurodist # We will use the dataset eurodist that gives the road distances (in km) # between 21 cities in Europe. Notice that eurodist is already # an object of class "dist"(matrix distance). euromat = as.matrix(eurodist) euromat[1:5, 1:5] # 1) MDS with cmdscale() # cmdscale(d, k = 2, eig = FALSE, add = FALSE, x.ret = FALSE) mds1 = cmdscale(eurodist, k = 2) plot(mds1[,1],mds1[,2],type="n",xlab="",ylab="",axes=FALSE,main="cmdscale (stats)") text(mds1[,1],mds1[,2],labels(eurodist),cex=0.9,xpd=TRUE) # 2) MDS with wcmdscale() # The package "vegan" provides the function wcmdscale() # (Weighted Classical Multidimensional Scaling). # Its general usage has the following form: # wcmdscale(d, k, eig = FALSE, add = FALSE, x.ret = FALSE, w) # If we specify the vector of the weights w as a vector of # ones, wcmdscale() will give ordinary multidimensional scaling. install.packages("vegan") library(vegan) mds2 = wcmdscale(eurodist,k=2,w=rep(1,21)) plot(mds2[,1],mds2[,2],type="n",xlab="",ylab="",axes=FALSE,main="wcmdscale (vegan)") text(mds2[,1],mds2[,2],labels(eurodist),cex=0.9,xpd=TRUE) # 3) MDS with pco() (package ecodist) # The package "ecodist" provides the function pco() # (Principal Coordinates Analysis). Its general usage has the following form: # pco(x, negvals = "zero", dround = 0) install.packages("ecodist") library(ecodist) mds3 = pco(eurodist) plot(mds3$vectors[,1], mds3$vectors[,2], type = "n", xlab = "", ylab = "",axes = FALSE, main = "pco (ecodist)") text(mds3$vectors[,1], mds3$vectors[,2], labels(eurodist),cex = 0.9, xpd = TRUE) #4) MDS with pco() (package labdsv) #The package "labdsv" also provides a function pco() (Principal Coordinates Analysis). #Its general usage has the following form: #pco(dis, k = 2) install.packages("labdsv") library(labdsv) mds4 = pco(eurodist, k = 2) plot(mds4$points[,1], mds4$points[,2], type = "n", xlab = "", ylab = "",axes = FALSE, main = "pco (labdsv)") text(mds4$points[,1], mds4$points[,2], labels(eurodist), cex = 0.9, xpd = TRUE) #5) MDS with pcoa() #The package "ape" provides the function pcoa() (Principal Coordinates Analysis). #Its general usage has the following form: #pcoa(D, correction="none", rn = NULL) install.packages("ape") library(ape) mds5 = pcoa(eurodist) plot(mds5$vectors[,1], mds5$vectors[,2], type = "n", xlab = "", ylab = "",axes = FALSE, main = "pcoa (ape)") text(mds5$vectors[,1], mds5$vectors[,2], labels(eurodist), cex = 0.9, xpd = TRUE) # 6) MDS with dudi.pco() #The package "ade4" provides the function dudi.pco() (Principal Coordinates Analysis). #Its general usage has the following form: #dudi.pco(d, row.w = "uniform", scannf = TRUE, nf = 2, full = FALSE, tol = 1e-07) install.packages("ade4") library(ade4) mds6 = dudi.pco(eurodist, scannf = FALSE, nf = 2) plot(mds6$li[,1], mds6$li[,2], type = "n", xlab = "", ylab = "",axes = FALSE, main = "dudi.pco (ade4)") text(mds6$li[,1], mds6$li[,2], labels(eurodist), cex = 0.9) # 7) MDS with smacofSym() # The package "smacof" provides the function smacofSym() (Multidimensional scaling # (stress minimization: SMACOF) on #symmetric dissimilarity matrix.). This function # uses a majorization approach to get the solution (more info in this ade4. # Its general usage has the following form: # smacofSym(delta, ndim = 2, weightmat = NULL, init = NULL, metric = TRUE, # ties = "primary", verbose = FALSE, #relax = FALSE, modulus = 1, itmax = 1000, eps = 1e-06) install.packages("smacof") library(smacof) mds7 = smacofSym(eurodist, ndim = 2) plot(mds7$conf[,1], mds7$conf[,2], type = "n", xlab = "", ylab = "",axes = FALSE, main = "smacofSym (smacof)") text(mds7$conf[,1], mds7$conf[,2], labels(eurodist),cex = 0.9, xpd = TRUE) install.packages("dlm") library(dlm) dados = read.csv("UKdriversKSI.csv",header=TRUE) y = as.numeric(dados[,1]) ts.plot(y) dlm3 = dlmModPoly() + dlmModSeas(12) buildFun = function(x) { diag(W(dlm3))[2:3] = exp(x[1:2]) V(dlm3) = exp(x[3]) return(dlm3) } fit3 = dlmMLE(y,parm=c(0.1,0.1,0.1),build=buildFun) dlm3 = buildFun(fit3$par) ySmooth = dlmSmooth(y, mod = dlm3) x = cbind(y, dropFirst(ySmooth$s[, c(1, 3)])) colnames(x) = c("KSI", "Trend", "Seasonal") plot(x, type = "o", main = "UK - Kill Seriously Injured") par(mfrow=c(1,2)) ts.plot(x[,1],ylim=range(x[,1:2])) lines(x[,2],col=2,lwd=3) par(mfrow=c(1,1)) plot(y[1:48],type="o") abline(v=12,lty=2) abline(v=24,lty=2) abline(v=36,lty=2) abline(v=48,lty=2) erro=y-x[,2]-x[,3] par(mfrow=c(1,3)) ts.plot(erro) acf(erro) pacf(erro) beta=0.25 t=1:100 e = rnorm(100) y = t*e plot(t,beta*t) plot(t,log(y))