rm(list = ls()) #Load required librarys library(dendextend) library(ggplot2) library(MASS) library(mlogit) library(foreign) library(nnet) library(stargazer) library(reshape2) library(plyr) library(scales) # The palette with black: cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") set.seed(54321) #do you want to calcualte interactions? interactions <- FALSE #Read data dataAll <- read.csv("./interview_data.csv",row.names=1,header=T) dataObjects <- read.csv("./object_coding.csv",row.names=1,header=T) #Load obects data #Hierachical clustering clustPhrases <- as.dendrogram(hclust(dist(dataAll[,1:41]),method = "ward.D2")) clustAtt <- as.dendrogram(hclust(dist(t(dataAll[,42:67])),method = "ward.D2")) #Plot the results of the hierarchical clustering pdf("dendro_phrases_all_horiz.pdf",width = 30, height = 70) par(mar=c(1,1,1,100)) plot(clustPhrases, ylab=NA,xlab=NA, sub=NA, main=NA,axes=FALSE,horiz = TRUE) rect.dendrogram(clustPhrases,k=8,horiz=TRUE,lower_rect = 0, border = "red") dev.off() pdf("dendro_attributes_horiz.pdf",width = 8, height = 11) par(cex = 1.7, mar=c(1,1,1,12)) plot(clustAtt, ylab=NA,xlab=NA, sub=NA, main=NA,axes=FALSE,horiz = TRUE) rect.dendrogram(clustAtt,k=4,horiz=TRUE,lower_rect = 0, border = "red") dev.off() par(mar=c(1,1,1,1)) #reset margins #Cut clustering solution to get processes cutPhrases <- cutree(clustPhrases,k=8) clusterTable <- split(row.names(as.matrix(cutPhrases)),as.matrix(cutPhrases)) # Find which clusters relate to which clips/systems clusterTableClipNo <- split(dataAll$Clip,as.matrix(cutPhrases)) clusterTableSystem <- split(dataAll$System,as.matrix(cutPhrases)) dfPhraseLevel <- data.frame(Phrase = clusterTable[7], Clip = clusterTableClipNo[7], System = clusterTableSystem[7]) dfPhrasePos <- data.frame(Phrase = clusterTable[5], Clip = clusterTableClipNo[5], System = clusterTableSystem[5]) cutLabels <- cutree(clustLabels,k=8) clusterLabelsTable <- split(row.names(as.matrix(cutLabels)),as.matrix(cutLabels)) #Make main data frame for regressions df <- data.frame(Process = as.numeric(cutPhrases), System = dataAll$System, Clip = dataAll$Clip, attEnv = dataAll$Enveloping, attHozWid = dataAll$Horizontal_width, attSpecBal = dataAll$Overall_spectral_balance, attSpatNat = dataAll$Spatial_naturalness, attDist = dataAll$Amount_of_distortion, attSpatClar = dataAll$Spatial_clarity, attSenSpac = dataAll$Sense_of_space, attSpatOpen = dataAll$Spatial_openness, attDepth = dataAll$Depth_of_field, attBand = dataAll$Bandwidth, attReal = dataAll$Realism, attSpecClar = dataAll$Spectral_clarity, attPhas = dataAll$Phasiness, attQualRev = dataAll$Subjective_quality_of_reverb, attRevLev = dataAll$Level_of_reverb, attSpecRes = dataAll$Spectral_resonances, attEnsBal = dataAll$Ensemble_balance, attSpatMov = dataAll$Spatial_movement, attSurr = dataAll$Surrounding, attClar = dataAll$Clarity, attHarsh = dataAll$Harshness, attPos = dataAll$Position_of_sound, attRich = dataAll$Richness_of_sound, attBass = dataAll$Bass, attDetail = dataAll$Detail, attEase = dataAll$Ease_of_listening ) #Merge the 2 position processes df$Process[df$Process == 5] <- 4 df$Process[df$Process == 6] <- 5 df$Process[df$Process == 7] <- 6 df$Process[df$Process == 8] <- 7 #Make copy of main data frame an turn attributes to factors df2 <- df df2[,names(df2)] <- lapply(df[,names(df2)],factor) #Set reference level for multinomial logit regression ref = 7 #7 sets No Change as the baseline #Reshape data frame for mlogit df3 <- mlogit.data(df2,shape="wide",choice="Process") #Check system-process relationships fit.mlog.system <-mlogit(Process ~ 1 | Clip + System, data = df3, reflevel = ref, method = "nr", print.level = 0) #Check individual mlogit regressions fit.mlog.attEnv <- mlogit(Process ~ 1 | attEnv, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attHozWid <- mlogit(Process ~ 1 | attHozWid, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attSpecBal <- mlogit(Process ~ 1 | attSpecBal, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attSpatNat <- mlogit(Process ~ 1 | attSpatNat, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attDist <- mlogit(Process ~ 1 | attDist, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attSpatClar <- mlogit(Process ~ 1 | attSpatClar, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attSenSpac <- mlogit(Process ~ 1 | attSenSpac, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attSpatOpen <- mlogit(Process ~ 1 | attSpatOpen, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attDepth <- mlogit(Process ~ 1 | attDepth, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attBand <- mlogit(Process ~ 1 | attBand, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attReal <- mlogit(Process ~ 1 | attReal, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attSpecClar <- mlogit(Process ~ 1 | attSpecClar, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attPhas <- mlogit(Process ~ 1 | attPhas, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attQualRev <- mlogit(Process ~ 1 | attQualRev, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attRevLev <- mlogit(Process ~ 1 | attRevLev, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attSpecRes <- mlogit(Process ~ 1 | attSpecRes, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attEnsBal <- mlogit(Process ~ 1 | attEnsBal, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attSpatMov <- mlogit(Process ~ 1 | attSpatMov, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attSurr <- mlogit(Process ~ 1 | attSurr, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attClar <- mlogit(Process ~ 1 | attClar, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attHarsh <- mlogit(Process ~ 1 | attHarsh, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attPos <- mlogit(Process ~ 1 | attPos, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attRich <- mlogit(Process ~ 1 | attRich, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attBass <- mlogit(Process ~ 1 | attBass, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attDetail <- mlogit(Process ~ 1 | attDetail, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.attEase <- mlogit(Process ~ 1 | attEase, data = df3, reflevel = ref, method = "nr", print.level = 0) #Try blockwise entry of remaining attributes fit.mlog.baseline <- mlogit(Process ~ 1 | Clip, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.block1 <- mlogit(Process ~ 1 | Clip + attEnv + attDepth + attEnsBal + attSpatMov, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.block2 <- mlogit(Process ~ 1 | Clip + attEnv + attDepth + attEnsBal + attSpatMov + attSpatClar + attSpatNat + attSurr + attSenSpac + attPos, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.block3 <- mlogit(Process ~ 1 | Clip + attEnv + attDepth + attEnsBal + attSpatMov + attSpatClar + attSpatNat + attSurr + attSenSpac + attPos+ attDetail + attSpecClar + attClar + attReal + attBass + attRich, data = df3, reflevel = ref, method = "nr", print.level = 0) fit.mlog.block4 <- mlogit(Process ~ 1 | Clip + attEnv + attDepth + attEnsBal + attSpatMov + attSpatClar + attSpatNat + attSurr + attSenSpac + attPos+ attDetail + attSpecClar + attClar + attReal + attBass + attRich + attEase + attDist + attQualRev + attHarsh, data = df3, reflevel = ref, method = "nr", print.level = 0) #Fit attribte-system models fit.system.attEnv = glm(attEnv ~ System, data = df2, family=binomial(link=logit)) fit.system.attDepth = glm(attDepth ~ System, data = df2, family=binomial(link=logit)) fit.system.attEnsBal = glm(attEnsBal ~ System, data = df2, family=binomial(link=logit)) fit.system.attSpatMov = glm(attSpatMov ~ System, data = df2, family=binomial(link=logit)) fit.system.attSpatClar = glm(attSpatClar ~ System, data = df2, family=binomial(link=logit)) fit.system.attSpatNat = glm(attSpatNat ~ System, data = df2, family=binomial(link=logit)) fit.system.attSurr = glm(attSurr ~ System, data = df2, family=binomial(link=logit)) fit.system.attSenSpac = glm(attSenSpac ~ System, data = df2, family=binomial(link=logit)) fit.system.attPos = glm(attPos ~ System, data = df2, family=binomial(link=logit)) fit.system.attDetail = glm(attDetail ~ System, data = df2, family=binomial(link=logit)) fit.system.attSpecClar = glm(attSpecClar ~ System, data = df2, family=binomial(link=logit)) fit.system.attClar = glm(attClar ~ System, data = df2, family=binomial(link=logit)) fit.system.attReal = glm(attReal ~ System, data = df2, family=binomial(link=logit)) fit.system.attBass = glm(attBass ~ System, data = df2, family=binomial(link=logit)) fit.system.attRich = glm(attRich ~ System, data = df2, family=binomial(link=logit)) fit.system.attEase = glm(attEase ~ System, data = df2, family=binomial(link=logit)) fit.system.attDist = glm(attDist ~ System, data = df2, family=binomial(link=logit)) fit.system.attQualRev = glm(attQualRev ~ System, data = df2, family=binomial(link=logit)) fit.system.attHarsh = glm(attHarsh ~ System, data = df2, family=binomial(link=logit)) glm_pval <- function(x){ #Function to calcualte model significane from a glm object an <- anova(x) pval <- pchisq(an$Deviance[2], an$Df[2], lower.tail = FALSE) return(pval) } #Plots of processes per system dataAll$Process <- df$Process dataAll$System <- df$System processLabel <- c("Spread", "EQ and processing", "Reverb", "Position", "Bass", "Level", "Nothing") systemLabel <- c("VBAP 4+5+0","VBAP 0+5+0","VBAP 0+2+0","Matrix 0+5+0","Matrix 0+2+0") plotTextSize <- 20 dfProcessSys <- data.frame(process = as.numeric(table(dataAll$Process)), label = processLabel) dfProcessSys$label = factor(dfProcessSys$label,levels=dfProcessSys[order(dfProcessSys$process),"label"]) #Set order as a factor, so that we can plot in ascending order p <- ggplot(data=dfProcessSys,aes(x = label,y = process)) + geom_bar(stat="identity") p <- p + theme_bw() p <- p + theme(plot.title = element_text(size = plotTextSize), axis.title = element_text(size = plotTextSize), axis.text.x=element_text(angle=45, hjust=1, size = plotTextSize),axis.text.y=element_text(size = plotTextSize)) + ylab("Count") + xlab("") pdf(paste("process_all_data.pdf",sep=""),width = 12, height = 10) print(p) dev.off() #individual plots for (n in 2:6){ dfProcessSys <- data.frame(process = as.numeric(table(dataAll$Process[which(dataAll$System == n)])), label = processLabel) dfProcessSys$label = factor(dfProcessSys$label,levels=dfProcessSys[order(dfProcessSys$process),"label"]) #Set order as a factor, so that we can plot in ascending order p <- ggplot(data=dfProcessSys,aes(x = label,y = process)) + geom_bar(stat="identity") p <- p + theme_bw() p <- p + theme(plot.title = element_text(size = plotTextSize), axis.title = element_text(size = plotTextSize), axis.text.x=element_text(angle=45, hjust=1, size = plotTextSize),axis.text.y=element_text(size = plotTextSize)) + ggtitle(paste("Frequency of process use for: ",systemLabel[n-1])) + ylab("Count") + xlab("") pdf(paste("process_",systemLabel[n-1],".pdf",sep=""),width = 12, height = 8) print(p) dev.off() } #grouped plot dfProcessSys <- df2 levels(dfProcessSys$System) = c("VBAP 4+5+0","VBAP 0+5+0","VBAP 0+2+0","Matrix 0+5+0","Matrix 0+2+0") levels(dfProcessSys$Process) = c("Spread", "EQ and processing", "Reverb", "Position", "Bass", "Level", "Nothing") p <- ggplot(dfProcessSys, aes(System, fill = Process)) + geom_bar(position="dodge") + theme_bw() p <- p + theme(legend.title=element_blank(), plot.title = element_text(size = plotTextSize), axis.title = element_text(size = plotTextSize), axis.text.x=element_text(angle=45, hjust=1, size = plotTextSize),axis.text.y=element_text(size = plotTextSize)) + ylab("Count") + xlab("") p <- ggplot(dfProcessSys, aes(System, fill = Process)) + geom_bar(position="fill",color="black") + theme_bw() p <- p + scale_y_continuous(labels = percent_format()) + theme(legend.text=element_text(size=15),legend.title=element_blank(),plot.title = element_text(size = 15), axis.title = element_text(size = 15), axis.text.x=element_text(angle=45, hjust=1, size = 15),axis.text.y=element_text(size = 15)) + theme(panel.border = element_blank()) p <- p + ylab("Cumulative percentage") + xlab("") p <- p + scale_fill_brewer(palette="OrRd") + coord_flip() pdf(paste("system_all_grouped.pdf",sep=""),width = 12, height = 10) print(p) dev.off() #Object type plots dataObjects$System <- df$System dataObjects$Clip <- df$Clip dataObjects$Process <- df$Process dfObjects <- data.frame(dataObjects) row.names(dfObjects) <- NULL dfObjects <- dfObjects[!dfObjects$Process==7,] #Remove 'Nothing' process data processLabel <- processLabel[1:6] objectAll <- colSums(dfObjects[,1:8]) objectProcess <- rowsum(dfObjects[,1:8],group = dfObjects$Process) objectSystem <- rowsum(dfObjects[,1:8],group = dfObjects$System) dfObjectProcess <- data.frame(prop.table(as.matrix(objectProcess),margin=2)*100,processLabel) dfObjectProcess.m <- melt(dfObjectProcess,id.vars = 'processLabel') dfObjectProcess.m$processLabel <- factor(dfObjectProcess.m$processLabel,levels(dfObjectProcess.m$processLabel)[c(6,2,5,4,1,3)]) cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") dfObjectProcess.m$variable <- revalue(dfObjectProcess.m$variable, c("Action_and_movement"="Action and movement", "Continuous_background" = "Continuous background", "Transient_background" = "Transient background", "Clear_speech" = "Clear speech", "Non_diegetic_music" = "Non diegetic music", "Prescence_of_people" = "Prescence of people", "Prominent_transients" = "Prominent transients", "No_specific_object_mentioned"="No specific object mentioned")) p <- ggplot(dfObjectProcess.m,aes(x = variable, y = value,fill = processLabel)) + geom_bar(position = "fill",stat = "identity",color="black") p <- p + scale_y_continuous(labels = percent_format()) + theme_bw() + theme(legend.title=element_blank()) p <- p + theme(legend.text=element_text(size=15),plot.title = element_text(size = 15), axis.title = element_text(size = 15), axis.text.x=element_text(angle=45, hjust=1, size = 15),axis.text.y=element_text(size = 15)) + theme(panel.border = element_blank()) p <- p + ylab("Cumulative percentage") + xlab("") p <- p + scale_fill_brewer(palette="OrRd") + coord_flip() pdf(paste("objects_processes_grouped.pdf",sep=""),width = 15, height = 10) print(p) dev.off()