### Jonathan Howell, Oct 22, 2010 ###
### updated March 10, 2011###

####  Some resources for learning R ####

#Linguists using R

Baayen, R. Harald. 2008. Analyzing Linguistic Data: A Practical Introduction to Statistics using R.   [I have a pre-publication pdf.]

Gries, Stefan Th. 2009.  Quantitative Corpus Linguistics with R: A Practical Introduction. Routledge.

Johnson, Keith. 2008. Quantitative Methods in Linguistics. Wiley-Blackwell.

R-Lang
A mailing list for language researchers who use R
https://mailman.ucsd.edu/mailman/listinfo/ling-r-lang-l


#General resources

Karatzoglou, Alexandros, David Meyer and Kurt Hornik.  2006. Support Vector Machines in R. Journal of Statistical Software 15:9. http://www.jstatsoft.org/v15/i09/paper  # compares libsvm with kernlab, klaR and svmpath
Venables, W.N. and B.D. Ripley. 2002. Modern Applied Statistics with S.  #especially good for graphical stuff


######## Here's some code for using SVMs in R #########

#################Read in datasets#########################
# thanIdid1 (corpus data from Everyzing)
# you'll need to change the path
thanIdid1 <- read.table("/yourdirectory/thanIdid1_with_zeros.dataframe",header=TRUE)

#thanIdid2 (corpus data from play.it)
thanIdid2 <- read.table("your_directory/thanIdid2_with_zeros.dataframe",header=TRUE)


#We can look at the first 5 observations for thanIdid1 with head()
#or the last 5 observations with tail()
#You might need to scroll back up to see the first columns
#Column 1 has the soundname; columns 2-309 have the acoustic features; column 310 has the focus category (subject focus (s) / non-subject focus (ns))
head(thanIdid1)
tail(thanIdid1)


#thanIdid1 has 91 observations with 310 features each
dim(thanIdid1)
#thanIdid2 has 127 observations with 310 features each
dim(thanIdid2)

#You can get a list of colnames with the colnames() function
colnames(thanIdid1)

#Use the $ symbol to call a particular column
thanIdid1$soundname
thanIdid1$duration_V2

#You can also use square brackets to call [row, column]
thanIdid1[,1]  #Column 1
thanIdid1[1,]  #Row 1 
#The colon allows you to specify a range and c() allows you to specify a list
thanIdid1[,c(1,3:5,310)]   #Columns 1, 3, 4, 5, 310

#Notice that we read in .dataframe files.  Here's how I created them
#from the .txt output from Praat and a .txt of category labels I created in Excel
#For example, let's read in some data collected in the lab
thanIdidlab <- read.table("/yourdirectory/thanididlab_withzeros.txt",header=T)
categories <- read.table("/yourdirectory/categories_thanIdid_lab.txt")

#I combine the two files with the column bind function cbind()
thanIdidlab <- cbind(thanIdidlab,categories)

#Finally, I create a dataframe file with write.table
write.table(thanIdidlab,"/yourdirectory/thanIdidlab.dataframe")



#######################Loading the library#################
##The R package e1071 contains libsvm
##(on other platforms the library is simply called libsvm)
# Install package e1071 using your GUI, if you don't already have it

#Load e1071 (once you've installed it)
library(e1071)

#################### Basic commands ###########################

#Specify the training data
train = thanIdid1[,-c(1)]	# the column with the soundname removed

#Create an svm model
#column to be predicted goes before the ~, and the predictor columns follow the ~
#Use . to specify "all the remaining data"

model1 <-svm(focus ~ ., data = train)

#Print a summary of the model
summary(model1)

#To see documentation for a command, put a question mark before it
?svm()


############## Understanding parameters and results ##############


#You can also call details of the model individually
#Here's a list of objects available
names(model1)

# The type is C-classification.  Other possibilities include:
# 	•	nu-classification
# 	•	one-classification (for novelty detection)
# 	•	eps-regression
# 	•	nu-regression


# The default kernel function (which maps attributes to another feature space) is RBF (radial basis function)
model1$kernel	#kernel type

#radial basis:  exp(-gamma*|u-v|^2)


# You can specify other kernel types in svm()
# kernels availabls: linear, polynomial, rbf, sigmoid
# See Hsu et al. (2003) "A Practical Guide to SVMs" for discussion
# http://ntur.lib.ntu.edu.tw/bitstream/246246/20060927122852476378/1/freiburg.pdf

#Try an SVM with a linear kernel
model1.linear <-svm(focus ~ ., data = train, kernel = "linear")
summary(model1.linear)

#An SVM with RBF kernel has 2 hyperparameters, cost and gamma
model1$cost; model1$gamma	#cost and gamma
#I believe a linear kernel has only 1 (cost), although both are listed in the summary
model1.linear$cost; model1.linear$gamma


# The default for gamma is 1/(data dimension)
1/309
# Hyperparameters can be tuned chosen using a "grid search" (again see Hsu et al.), although despite the library authors' warnings that "parameters of SVM-models usually must be tuned to yield sensible results!" I've had excellent results without performing this.

# This is the code you'd use to do the tuning.  Note that you must specify the parameter ranges that you want to search within.
obj <- tune(svm, focus ~. , data = thanIdid1, ranges = list(gamma = 2^(-1:1), cost = 2^(2:4)), tunecontrol = tune.control(sampling = "fix"))
# This tuning was not terribly successful (best cross-validation performance 51.6%)
summary(obj)

#You can specify parameters direcetly in the svm() command
model1.tuned <- svm(focus ~ ., data = train, gamma = 0.5, cost = 4)

#Getting back to other objects of the svm() command,
summary(model1)
# Note that this model used 84(!) support vectors (data instances right on the margins), of which 41 belong to class "ns" and 43 belong to class "s".
model1$nSV

#You can call the actual data occurrences which are the support vectors with this command
model1$SV
#Index gives the indices (in the orginal dataframe) for the data instances which are support vectors in the model
model1$index
#In this way, it is possible to know the filenames of the individual support vectors
thanIdid1[model1$index,1]

#In the R implementation of libsvm, data are scaled internally to zero mean and unit variance.  You can get the centre value for each feature with this command
model1$x.scale



#################### Cross-validation ###################
# Initially, we used cross-validation, although in my dissertation\, I just use performance on new data (see below).

# There is a k-fold cross-validation option in svm()
# Here I have selected 10-fold cross-validation
model1.cv <- svm(focus ~ ., data = train, cross=10)
summary(model1.cv)


# Run the command above more than once, and you'll see that the single accuracies and total accuracy will of course be different every time you run it.

#To avoid this, you can do a one-held-out (aka "jack-knifed") cross-validation, although there is no feature for it, so I wrote the code below.

##### One-held-out x-validation
##Run all of the following code at the same time

dataset <- thanIdid1
#Create vector for predicted classification
predicted <- vector(mode="character",length = nrow(dataset))

#Loop through every data instance (i.e. row) in dataset
for (i in 1:nrow(dataset)) {	
	#Get number of data instances (rows)
	index <- 1:nrow(dataset)   
	#Make this data instance the test case
	testindex = i
	#Identify the column that has the name of the soundfile
	soundnameindex <- which(colnames(dataset)=="soundname")
	#Create testset (the single data instance, and exclude filename)
	testset <- dataset[testindex, -soundnameindex]
	#Create trainset (the other data instances, and exclude filename)
	trainset <- dataset[-testindex, -soundnameindex]  
	#Create svm classifier
	svm.model <- svm(focus ~ ., data = trainset)
	#Identify the column containing the human classification
	focusindex <- which(colnames(dataset)=="focus")
	#Get the classifier's prediction for the testset (need to exclude human classification column)
	svm.pred <- predict(svm.model, testset[, -focusindex])
	#Keep track of progress in the terminal
	print(svm.pred[1])
	#Add prediction to the new vector (interpret as character)
	predicted[i] <- as.character(svm.pred)
	}
#Add prediction vector to dataset matrix	
dataset$predicted <- predicted

#Print a contingency table of predicted against true
table(pred=dataset$predicted,true=dataset$focus)

#Print classification accuracy and balanced error rate
accuracy <- (length(dataset$focus)[1]-(table[2]+table[3]))/length(dataset$focus)[1]*100; round(accuracy,1)
BER <- 0.5*((table[3]/(table[1]+table[3]))+(table[2]/(table[2]+table[4])))*100; round(BER,1)

#############

############## Testing on novel data #########

#Specify training data and vector of correct responses
test <- thanIdid2[,c(-1,-310)]
correct <- thanIdid2$focus

#In order to test a model on new data, use the predict() function
pred <- predict(model1,test)

#Print contingency table
table <- table(pred,correct); table

#Print classification accuracy
accuracy <- (dim(test)[1]-(table[2]+table[3]))/dim(test)[1]*100; round(accuracy,1)

#Print balanced error rate
BER <- 0.5*((table[3]/(table[1]+table[3]))+(table[2]/(table[2]+table[4])))*100; round(BER,1)


############## Making a custom function #########

#Since I needed to test a lot of training/test pairs, I wrote a little function for it.
# Read in the code below and then you can use the function get.values.svm()

#### start of function ###

	#arguments should have 310 columns
get.values.svm <- function(train,test,features){

	model.rbf <- svm(focus ~ ., data = train[,c(features,310)])
	pred.rbf <- predict(model.rbf,test[,c(-1,-310)])
	table.rbf <- table(pred.rbf,test$focus)
	accuracy.rbf <- round((dim(test)[1]-(table.rbf[2]+table.rbf[3]))/dim(test)[1]*100,1)
	BER.rbf <- round(0.5*((table.rbf[3]/(table.rbf[1]+table.rbf[3]))+(table.rbf[2]/(table.rbf[2]+table.rbf[4])))*100,1)

	model.lin <-svm(focus ~ ., data = train[,c(features,310)], kernel = "linear")
	pred.lin <- predict(model.lin,test[,c(-1,-310)])
	table.lin <- table(pred.lin,test$focus)
	accuracy.lin <- round((dim(test)[1]-(table.lin[2]+table.lin[3]))/dim(test)[1]*100,1)
	BER.lin <- round(0.5*((table.lin[3]/(table.lin[1]+table.lin[3]))+(table.lin[2]/(table.lin[2]+table.lin[4])))*100,1)
	result = cat("SVM-RBF","\n",accuracy.rbf," (",BER.rbf,")","\n","SVM-lin","\n",accuracy.lin," (",BER.lin,")", sep="")

	print(result)
	}
### end of function ####
#############################################

#The function get.values.svm() expects three arguments in order: training set, test set, and vector with feature indices.

# For example:
f.all <- c(2:309)    # features 2 through 209

get.values.svm(thanIdid1, thanIdid2, f.all)

#Not sure what the "NULL" is doing there, but the output is the classification accuracy and balanced error rate for an SVM with RBF kernel and SVM with linear kernel.

# For a list of features, use this command
colnames(thanIdid1)

# So, if I want to test a model with duration_V2 and duration_V3, I can put in the corresponding indices

f.voweldurations <- c(2,3)

# and plug them into the function

get.values.svm(thanIdid1, thanIdid2, f.voweldurations)

# Most likely, the best model for your data will not use all variables, but instead some subset.  There are many, many ways to do pick a subset of variables (aka feature selection, feature elimination, feature filtering, etc.)

# My most successful method of feature selection was human trial-and-error, although I also found the package/library varSelRF useful.



######################################
############## Plotting Data #########
######################################

#### Plotting within libsvm  ####
# There is a plot function within libsvm, although it is limited to 2 dimensions, does not plot the decision function, and is difficult/impossible to customize.  See the examples in 
?plot.svm()

#See also other svm packages (e.g. Karatzoglou et al. above)

### Scatterplot matrices ###
# Plotting in more than two dimensions is complicated.  One way to visualize the distribution of multivariate data is to plot variables 2x2 in a matrix of scatterplots.
#I use the function splom() "S(catter)PL(ot)M(atrix)", which is found in the lattice package.

#Load package lattice (once you've installed it)
library(lattice)
#Documentation for splom() here:
?splom()

# Let's start with the basic command, and then slowly customize the graphics.
# Plot the variables "duration_V2", "duration_V3", "f0_ratio" and "f1f2Time50_V2" (indices 2,5,17,179) against each other
splom(~thanIdid1[c(2,5,17,179)])

#Group the data according to the focus label
splom(~thanIdid1[c(2,5,17,179)], data = thanIdid1, groups = focus)

#Specify the colours of the two groups
splom(~thanIdid1[c(2,5,17,179)], data = thanIdid1, groups = focus, col=c("red","blue"))

#A list of 657 is available with this command
colors()

#Specify the variable names (e.g. I wanted to shorten mine), 
splom(~thanIdid1[c(2,5,17,179)], data = thanIdid1, groups = focus, col=c("red","blue"), varnames = list("dur_V2","dur_C3", "f0_ratio","f1f2_50_V2"))

#Specify the font size of the variable names
splom(~thanIdid1[c(2,5,17,179)], data = thanIdid1, groups = focus, col=c("red","blue"), varnames = list("dur_V2","dur_C3", "f0_ratio","f1f2_50_V2"), varname.cex = .8)

#Specify the font size of the axes
splom(~thanIdid1[c(2,5,17,179)], data = thanIdid1, groups = focus, col=c("red","blue"), varnames = list("dur_V2","dur_C3", "f0_ratio","f1f2_50_V2"), varname.cex = .8, axis.text.cex = .5)

#Add a title and legend
splom(~thanIdid1[c(2,5,17,179)],col=c("red","blue"), aspect = 9/12, groups = focus, data = thanIdid1,
		varnames = list("dur_V2", "dur_C3", "f0_ratio", "f1f2_50_V2"), varname.cex = .8,
		axis.text.cex = .5,
      key = list(title = "Model 1",
                 columns = 2, 
                 points = list(pch = c(1), col=c("red","blue")),
                 text = list(c("s", "ns")
                 ))
      )


#Write the graphic to a file (e.g. pdf, postscript, jpg)
pdf("/yourdirectory/image_name.pdf")

splom(~thanIdid1[c(2,5,17,179)],col=c("red","blue"), aspect = 9/12, groups = focus, data = thanIdid1,
		varnames = list("dur_V2", "dur_C3", "f0_ratio", "f1f2_50_V2"), varname.cex = .8,
		axis.text.cex = .5,
      key = list(title = "Model 1",
                 columns = 2, 
                 points = list(pch = c(1), col=c("red","blue")),
                 text = list(c("s", "ns")
                 ))
      )
dev.off()


#####################
### RGGobi
#####################

#It is also possible to model the svm model in 3D, although it requires a fair bit of setup work. See the following resources for GGobi and its R interface RGGobi here:

Book available (free) online:
Interactive and Dynamic Graphics for Data Analysis: With Examples Using R and GGobi. Dianne Cook and Deborah F. Swayne
http://www.ggobi.org/
http://www.ggobi.org/rggobi/


#Load the necessary packages
library(rggobi)
library(classifly)
library(lattice)

#Here is the code I use to invoke GGobi
#There are three predictor variables and the kernel is radial
classifly(thanIdid1,focus~duration_V2+duration_C3+f1f2Time50_V2,best.svm,probability=TRUE,kernel="radial")

#In the terminal, View > 2D Tour allows one to rotate the data in 3D
#Interaction > Identify allows one to know the value of an observation by hovering the cursor over it

#Have fun playing, and remember how to call the help for this function.

?classifly()






