####Script for basic Statistical Analyses of Nominal Variables#####
##########Workshop New Technologies in Linguistic Studies##########
########################Livia Oushiro (2018)#######################

#This file may serve as a basis for your own analyses in R. Copy and paste relevant parts of this script into a new file, save it and adapt it to your data (names of variables, variants, datasets etc.)

#####ANALYSIS OF A NOMINAL VARIABLE#####

#Define working directory
setwd() # !!! 

#Load data
ds <- read.csv("LabovDS.csv", header = T, sep = ",")

#Check data
str(ds)
head(ds)
summary(ds)

####(0) ADJUSTING DATA#####

ds2 <- subset(ds, r %in% c("r0", "r1"))
str(ds2)

ds2$r <- factor(as.character(ds2$r), exclude = "d")
ds2$r

####(1) TABULATING DATA#####

####Frequencies: table()####
r.dist <- with(ds2, table(r))
r.dist

store.dist <- with(ds2, table(store, r))
store.dist

#------> Create a table to visualize the frequencies of r variants in each emphasis context



####Proportions: prop.table()####
r.prop <- with(ds2, prop.table(r.dist))
r.prop

store.prop <- with(ds2, prop.table(store.dist, 1)) * 100 # 1 = per row; 2 = per column; default = whole table
store.prop

#------> Create a proportions table to visualize the proportions of r variants for each word.





####(2) VISUALIZING DATA#####

####Bargraphs: barplot()####

barplot(r.prop, beside = T)

barplot(store.prop, beside = T)

barplot(store.prop, 
        beside = T, 
#        width=c(0.5, 0.5, 0.5, 0.5, 0.5, 0.5),
        horiz = F, 
        main = "Distribution of (-r) in three department stores", 
        xlab = "(-r) variants", 
        ylab = "Proportion %", 
        names.arg = sort(unique(ds2$r)), 
        legend.text = T,
#       xlim = c(0,10), 
        ylim = c(0,100), 
        cex.axis = 1, 
        cex.names = 1, 
        col = c("hotpink", "goldenrod1", "cyan1")
) 

#For more colors, check out http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf

#Customizing the legend:
colors <- c("hotpink", "goldenrod1", "cyan1")  
legend("topright", #Options: "topright", "topleft", "bottomleft", "bottomright"
       c("Klein's", "Macy's", "Saks"), 
       cex = .8, 
       bty = "n", 
       fill = colors) 

#------>What if you'd like the bars to have different colors?

#------>What if you'd like the plot to show only the tokens of r0?

#------>What if you'd like the plot to show the stores in the reverse order (Saks, Macys, Kleins)?

#------>What if you'd like the bars to be named Saks, Macy's and Klein's (with apostrophe)?


####Line charts: plot()####

#plot data
plot(store.prop[3:1, 1], 
     type = "o", 
     pch = 19, 
     lty = 3, 
     axes = F, 
     xlab = "Department stores", 
     ylab = "Proportion of r0", 
     col = "black")

#Define x axis
axis(1, at = 1:3, lab = c("Saks", "Macy's", "S. Klein"))

#Define y axis
axis(2, las = 2)

#Insert box around the figure
box()

#Insert title
title("Proportion of post-vocalic /r/-deletion in three New York City departament stores")

#For more on graphs, check out:
#http://tutorials.iq.harvard.edu/R/Rgraphics/Rgraphics.html
#www.r-graph-gallery.com

####(3) HYPOTHESES####

#The binomial distribution####

#probability of obtaining from 0 to 3 heads in 3 tosses
bb<-dbinom(0:3, 3, .5); bb
cc<-barplot(bb, beside=T, ylim=c(0,.5),names.arg=c("0","1","2","3"))
text(x=cc[,1], y=bb+0.05, labels=bb)

#probability of obtaining from 0 to 6 heads in 6 tosses
bb<-dbinom(0:6, 6, .5); bb
cc<-barplot(bb, beside=T, ylim=c(0,.4),names.arg=c("0","1","2","3", "4", "5", "6"))
text(x=cc[,1], y=bb+0.05, labels=round(bb,3))

#probability of obtaining from 0 to 12 heads in 12 tosses
bb<-dbinom(0:12, 12, .5); bb
cc<-barplot(bb, beside=T, ylim=c(0,.3))
text(x=cc[,1], y=bb+0.05, labels=round(bb,2))

#probability of obtaining from 0 to 25 heads in 25 tosses
bb<-dbinom(0:25, 25, .5); bb
cc<-barplot(bb, beside=T, ylim=c(0,.15))
#text(x=cc[,1], y=bb+0.05, labels=round(bb,3))

#probability of obtaining from 0 to 50 heads in 50 tosses
bb<-dbinom(0:50, 50, .5); bb
cc<-barplot(bb, beside=T, ylim=c(0,.10))
#text(x=cc[,1], y=bb+0.05, labels=round(bb,3))

#probability of obtaining from 0 to 100 heads in 100 tosses
bb<-dbinom(0:100, 100, .5); bb
cc<-barplot(bb, beside=T, ylim=c(0,.08))
#text(x=cc[,1], y=bb+0.05, labels=round(bb,3))


#Another distribution, with probability = .31 instead of .50
bb<-dbinom(0:100, 100, .31); bb
cc<-barplot(bb, beside=T, ylim=c(0,.08))
#text(x=cc[,1], y=bb+0.05, labels=round(bb,3))

#From your perspective
sum(dbinom(55:100, 100, 0.5))
sum(dbinom(59:100, 100, 0.5))
qbinom(0.05, 100, 0.5, lower.tail=FALSE)

#From an independent observer's perspective
qbinom(0.025, 100, 0.5, lower.tail=TRUE); qbinom(0.025, 100, 0.5, lower.tail=FALSE)

####(4) BASIC INFERENTIAL ANALYSES#####

#Proportions test: prop.test()####

#First, you'd have done some graphic explorations and raised hypotheses...
r.dist <- with(ds2, table(r))
r.dist
r.prop <- with(ds2, prop.table(r.dist))
r.prop
barplot(r.dist, beside = T)

#To test whether the distribution between these variants differ significantly from 50%-50%
prop.test(r.dist)


?prop.test
#-------> What if you want to test the null hypothesis that p = 0.7?

#-------> What if you want to test the null hypothesis that p < 0.7?

#-------> What if you want to test the null hypothesis with a confidence level = 0.99?

#Chi-square test: chisq.test()####

#For more complex distributions envolving more than one variable, one runs a chi-square test
#First, you'd have done some graphic explorations and raised hypotheses...
store.dist <- with(ds2, table(store, r))
store.dist
store.prop <- with(ds2, prop.table(store.dist, 1)) 
store.prop

barplot(store.prop[, 1], beside = T)
abline(h = .68, lty = 2)

#Testing whether distributions among stores are significantly different
x2.store <- chisq.test(store.dist)
x2.store
#If you're practical, this is all you need to run a chi-square test! Make a frequency table and apply chisq.test() function! But if you want to know more about this test...

#R has calculate a lot more than the chi-square output...
str(x2.store)

#The test measures the difference between the observed and an expected distribution (as do all statistical tests) 
O <- x2.store$observed; O
E <- x2.store$expected; E

addmargins(O)
#Compare: (row total * column total)/total N. E.g.:
(499 * 216) / 730 #147.6493
147.6493 / 216

#This is the chi-square formula in R; compare the result with the test output
sum((O - E) ^ 2 / E)

#The residuals measure the difference between observed and expected values. The greater the residuals, the more significant the difference.
x2.store$residuals


#-------> What if you want to test if there's a significant difference only between Macy's and Saks?


#A multivariate analysis with logistic regression: glm()####

model1 <- glm(r ~ store + emphasis + word, data = ds2, family = binomial)
summary(model1)

#(If you have random predictors, such as speaker and lexical item, you should run mixed-effects models with glmer(). See Levshina 2015.)