####################################################################################
#This is a file with executable R code for chapter 17 of Natalia Levshina's (2015) 
#How to Do Linguistics with R. Amsterdam/Philadelphia: John Benjamins. 
####################################################################################


###Section 17.1

##Main text

install.packages("rworldmap")

library(Rling); library(rworldmap)

data(eWAVE)

str(eWAVE)

rownames(eWAVE)[1:10]

code <- as.numeric(eWAVE$Type)

code

dialmap <- getMap()

plot(dialmap)

points(eWAVE[, 23:24], pch = code + 20, col = code, bg = code)

legend("bottomleft", legend = levels(eWAVE$Type), pch = 21:25, col = 1:5, pt.bg = 1:5)

plot(dialmap, xlim = c(-140, -40), ylim = c(30, 50), asp = 1)

points(eWAVE[c(3, 18), 23:24], col = code[c(3, 18)], cex = 1.2, pch = code[c(3, 18)] + 20, bg = code[c(3, 18)])

text(eWAVE[c(3, 18), 23:24], labels = rownames(eWAVE)[c(3, 18)], cex = 0.8, adj = c(0.5, 1.2))


###Section 17.3
	
##Main text

install.packages(c("fields", "rgl", "MASS"))

library(Rling); library(fields); library(rgl); library(MASS)

data(eWAVE)

geo.dist <- rdist.earth(eWAVE[, 23:24], miles = FALSE)

geo.dist <- as.dist(geo.dist)

geo.mds <- cmdscale(geo.dist, eig = TRUE)

barplot(geo.mds$eig[1:20], xlab = "Number of dimensions", ylab = "Eigenvalues", main = "Scree plot")

plot(geo.mds$points, type = "n", main = "MDS of geographic distances between varieties of English")

text(geo.mds$points, labels = rownames(eWAVE), cex = 0.6)

geo.mds.3d <- cmdscale(geo.dist, k = 3, eig = TRUE)

plot3d(geo.mds.3d$points, type = "n")

text3d(geo.mds.3d$points, texts = rownames(eWAVE), cex = 0.6)

geo.mds.3d$GOF

sqrt(sum((geo.dist - dist(geo.mds.3d$points))^2)/sum(geo.dist^2))

geo.sh <- Shepard(geo.dist, geo.mds.3d$points)

plot(geo.sh, main = "Shepard plot", pch = ".")

lines(geo.sh$x, geo.sh$yf, type = "S")


###Section 17.4

##Main text

install.packages(c("cluster", "MASS", "ggplot2")) # if you haven't installed them yet

library(Rling); library(cluster); library(MASS); library(ggplot2)

data(eWAVE)

eWAVE1 <- eWAVE

eWAVE1[eWAVE1 == "?"] <- NA

eWAVE1[eWAVE1 == "X"] <- NA

eWAVE2 <- lapply(eWAVE1[, 1:20], function(x) ordered(x, levels = c("D", "C", "B", "A")))

eWAVE2 <- data.frame(eWAVE2)

str(eWAVE2)

ling.dist <- daisy(eWAVE2)

summary(ling.dist)

ling.dist1 <- ling.dist

ling.dist1[is.na(ling.dist)] <- mean(ling.dist, na.rm = TRUE)

ling.dist2 <-ling.dist1

ling.dist2[ling.dist2 == 0] <- 0.001

ling.mds.kr <- isoMDS(ling.dist2)

qplot(ling.mds.kr$points[, 1], ling.mds.kr$points[, 2], label = rownames(eWAVE), cex = 0.5, col = eWAVE$Type, geom = "text")


###Section 17.5

##Main text

install.packages("vegan")

library(vegan)

mantel(geo.dist, ling.dist, na.rm = TRUE)