# This script loads a data spreadsheet and creates the R data file that is loaded by the analysis script # # Marc Badger # "Avoiding topsy-turvy: how Anna's Hummingbirds (Calypte anna) fly through upward gusts" # Marc Badger, Hao Wang, Robert Dudley # June 2018 # SET YOUR OWN WORKING DIRECTORY! setwd("C:/Users/USERNAME/Google Drive/Gust Manuscript/For Dryad") rm(list=ls()) library(R.matlab) library(xlsx) # Import data directly from Matlab struct: matlabFile <- readMat('gustDataForR_2018-06-25_21-11-15.mat') colNames <- names(matlabFile$dataStructForR[,,1]) datList <- matlabFile$dataStructForR datList <- lapply(datList, unlist, use.names=FALSE) AIRKNIFE <- as.data.frame(datList) names(AIRKNIFE) <- colNames colsToFactor <- c('noMBAORAirOff','MBABeTcCm','MBABeTcLg','MBAAfTrEnd','airKnifeOn') AIRKNIFE[colsToFactor] <- lapply(AIRKNIFE[colsToFactor] , factor) #str(AIRKNIFE, list.len=149) AIRKNIFE['WTreatTNSC'] <- (AIRKNIFE$trialNumber-mean(AIRKNIFE$trialNumber))/(max(AIRKNIFE$trialNumber)-mean(AIRKNIFE$trialNumber)) AIRKNIFE['isMBAinAF'] <- factor(1*!(AIRKNIFE$tminBodyAngleSegment == 'AfTcLg' | AIRKNIFE$tminBodyAngleSegment == 'AfTrEnd')) levels(AIRKNIFE$birdID)[levels(AIRKNIFE$birdID)=="b4"] <- "b3" # rename bird levels so they are 1-4 instead of 1,2,4,5 levels(AIRKNIFE$birdID)[levels(AIRKNIFE$birdID)=="b5"] <- "b4" # Variable for minBodyAngularAccelerationAfterBBCrossMidline AIRKNIFE["minPitchAccel"] <- AIRKNIFE["minBodyAngularAcceleration"] # A variable for minimum pitch angular acceleration after entering the air gust AIRKNIFE[AIRKNIFE$tminBodyAngularAccelerationSegment == "BeBbCm","minPitchAccel"] <- AIRKNIFE[AIRKNIFE$tminBodyAngularAccelerationSegment == "BeBbCm","bodyAngularAccelerationInBodyAtBBCrossMidline"] str(AIRKNIFE, list.len=200) # Data frames used: # AIRKNIFE - all observations in the dataset # AIRKNIFEGU - all transits in which the air gust was on # AIRKNIFEGUG - all transits whith air knife onn where birds recovered minimum body angle before going off screen # AIRKNIFEGUGOR - like AIRKNIFEGUG, but with trial XXX (an outlier) removed # AIRKNIFEOR - All trials (AIRKNIFE), but with trial XXX (an outlier) removed # AKDAT - AIRKNIFE, but with a subset of variables with abbreviated names for ease of use # AIRKNIFEG is all trials in which minimum body angle found before the end of recording (MBAAfTrEnd == 0) AIRKNIFEG <- subset(AIRKNIFE, MBAAfTrEnd == 0) # AIRKNIFEGU is all trials in which the air was on AIRKNIFEGU <- subset(AIRKNIFE, airKnifeOn == 1) #airKnifeOn:1=on,0=off # AIRKNIFEGUG is all trials in which the air was on and minimum body angle found before the end of recording (MBAAfTrEnd == 0) AIRKNIFEGUG <- subset(AIRKNIFE, noMBAORAirOff == 0 & airKnifeOn == 1) # Can remove outliers if needed e.g.: # AIRKNIFEGUGOR <- subset(AIRKNIFEGUG, filename != 'b2_jUP_p18psi_t010_dB_xyzpts.csv') # AIRKNIFEGUGOR["DURSC"]<-scale(AIRKNIFEGUGOR[,"durationtLg"]) # rescale the durationtLg using only the selected trials AIRKNIFEGU <- droplevels(AIRKNIFEGU) AIRKNIFEGUG <- droplevels(AIRKNIFEGUG) # Variables were averaged over several different intervals, which can be selected using the last part of the name. # E.g. MAfTcCmBeTcLg is mean from after tail center crossed gust midline but before tail center left gust # MAfBbCmBeTcLg is mean from after bill base crossed gust midline but before tail center left gust # Create a dataset with selected variables and abbreviated names for ease of use responseVars <- c("WTreatTNSC", "zBBAtBBCrossMidline", "vxBBAtBBCrossMidline", "vzBBAtBBCrossMidline", "bodyAngleAtBBCrossMidline", paste("tailExtensionAngle", "MAfTcCmBeTcLg", sep=''), paste("tailElevation", "InBody", "MAfTcCmBeTcLg", sep=''), paste("tailElevation", "InBirdZ", "MAfTcCmBeTcLg", sep=''), paste("rawWingElevationAngle", "InBody", "MAfBbCmBeTcLg", sep=''), paste("rawWingSweepAngle", "InBody", "MAfBbCmBeTcLg", sep=''), paste("proportionSweptBackAndElevated", "InBody", "MAfBbCmBeTcLg", sep=''), paste("instantaneousVX","MAfBbCmBeTcLg", sep=''), "durationtLg", "minBodyAngle", "minPitchAccel", "deltaHeightLg", "wingPhaseAtBBCrossMidline", "birdID", "flightDirection", "airKnifeOn", "filename") AKDAT <- AIRKNIFE[, responseVars] colnames(AKDAT) <- c("TN","ZB","VX","VZ","IBA","TF","TAB","TAH","EL","SW","WH","ATS", "DU","MBA","MPA","DH","IWP","BID","FD","AKO","FN") levels(AKDAT$AKO) <- c("AirKnifeOff", "AirKnifeOn") # Rescale TN so it corresponds to within bird x treatment trial number # TN will range from -0.5 to 0.5 within each bird x treatment combination. # We need to do this in order to have meaningful interpretation of other coefficients # when TN or TN interactions are in the model. thisAKO <- levels(AKDAT[,"AKO"]) thisBID <- levels(AKDAT[,"BID"]) for (aNum in 1:2) { for (bNum in 1:4) { thisDat <- AKDAT[AKDAT$BID == thisBID[[bNum]] & AKDAT$AKO == thisAKO[[aNum]], "TN"] if (min(thisDat) == max(thisDat)) { AKDAT[AKDAT$BID == thisBID[[bNum]] & AKDAT$AKO == thisAKO[[aNum]], "TN"] <- rep(0, length(thisDat)) } else { AKDAT[AKDAT$BID == thisBID[[bNum]] & AKDAT$AKO == thisAKO[[aNum]], "TN"] <- (thisDat - min(thisDat))/(max(thisDat)-min(thisDat)) - 0.5 } } } # TN = TrialNumber (within treatment) # ZB = Vertical position of beak base above air knife opening at beak base # crossing gust midline ("gust entry") (in cm) # VX = Horizontal velocity of beak base at gust entry (in m/s) # VZ = Vertical velocity of beak base at gust entry (in m/s) # IBA = Initial body angle at gust entry (in deg above horizontal) # TF = Mean tail fan angle from tail center cross gust midline to tail center leave gust # (time region of greatest variation in tail angles) (in deg from centerline) # TAB = Mean tail elevation angle from tail center cross gust midline to tail center leave gust (deg from body axis) # TAH = Mean tail elevation angle from tail center cross gust midline to tail center leave gust (deg from horizontal) # EL = Mean wing elevation angle in body frame (+ = dorsal, - = ventral) during # gust transit (from beak base cross gust midline to tail center leave gust) (in deg) # SW = Mean wing sweep angle in body frame (+ = forward toward head, - = backward toward tail) # during gust transit (in deg) # WH = Percent of transit durationtLg that wings were both elevated and swept backward # ATS = Average transit speed - mean instantaneous horizontal speed from # when bill base crossed midline to tail center left gust (in m/s) # DU = Transit duration (from beak base cross gust midline to tail center leave gust) (in seconds) # MBA = Minimum body angle (proxy for intensity of perturbation) (in degrees) # MPA = Minimum pitch angular acceleration in body frame (proxy for intensity of perturbation) (in deg/s^2) # DH = Height gained during time from bill base cross midline to tail center leave gust (in cm) # IWP = Initial wing phase at gust entry (in degrees) # BID = Bird ID # AKO = Air knife on (AirKnifeOn) or off (AirKnifeOff) # Preliminary clustering to look at distributions of values library(mclust) # Cluster using variables that could be important for entry conditions, technique, and performance: standardizedAKDAT <- as.data.frame(scale(AKDAT[c("ZB","IBA","TF","TAH","EL","SW","ATS","MBA","MPA","IWP")])) fit <- Mclust(standardizedAKDAT) # or "VVI", "VVE", "VII" summary(fit) # plot(fit) # Three clusters selected via BIC # We see that clusters are most separated in tail fan and elevation angles. Wing sweep, minimum body angle, and minimum pitch angular acceleration also appear to differ among clusters. CLUSTDAT <- data.frame(AKDAT, Technique=factor(fit$classification)) levels(CLUSTDAT$Technique) <- c("Control", "WingsDominated", "TailDominated") createPlots <- FALSE if (createPlots) { # Create some descriptive plots (warning: can take a minute to run) library(GGally) # Run one of these: subsetToPlot <- CLUSTDAT[,c("BID","AKO","FD","ZB","VX","IBA","Technique")] subsetToPlot <- CLUSTDAT[,c("TF","TAH","EL","SW","WH","Technique")] subsetToPlot <- CLUSTDAT[,c("IWP","MBA","MPA","ATS","Technique")] subsetToPlot <- CLUSTDAT[,c("BID","ZB","IBA","TF","TAH","WH","MBA","Technique")] subsetToPlot <- CLUSTDAT[,c("BID","ZB","VX","IBA","TF","TAH","EL","SW","WH","MBA","Technique")] levels(subsetToPlot$Technique) = c("C", "WD","TD") # Then one of these: ggpairs(subsetToPlot, aes(colour = Technique, alpha = 0.4)) ggpairs(subsetToPlot, aes(colour = BID, alpha = 0.4)) } ############################ ##### HELPER FUNCTIONS ##### ############################ theme_simple <- function (base_size = 12, base_family = "Helvetica") { theme_gray(base_size = base_size, base_family = base_family) %+replace% theme( axis.text = element_text(colour = "black"), axis.title.x = element_text(colour = "black", size=rel(1)), axis.title.y = element_text(colour = "black", angle=90), axis.line.x = element_line(colour = "black", size=rel(1)), axis.line.y = element_line(colour = "black", size=rel(1)), axis.ticks=element_line(colour = "black"), panel.background = element_blank(), panel.grid= element_blank(), plot.background = element_blank(), legend.key=element_rect(colour="white") ) } findBestVar <- function(my.respVar, my.fitForm, my.varList, my.d) { glsEqualVar <- do.call('gls', list(model=as.formula(paste(my.respVar, "~", my.fitForm, collapse=" ")), data = substitute(my.d), na.action = na.omit)) modelList <- list() modelAICs <- list() modelLRTs <- list() returnList <- list() for (i in 1:length(my.varList)) { modelList[[i]] <- do.call('gls', list(model=as.formula(paste(my.respVar, "~", my.fitForm, collapse=" ")), weights = substitute(eval(parse(text=eval(my.varList[i])))), data = substitute(my.d), na.action = na.omit)) modelAICs[[i]] <- AICc(modelList[[i]]) theLRT <- anova(glsEqualVar, modelList[[i]]) modelLRTs[[i]] <- theLRT$'p-value'[2] } if (min(p.adjust(unlist(modelLRTs), method = "BY")) > 0.05) { returnList[[1]] <- glsEqualVar returnList[[2]] <- "NA" return(returnList) } else { minPos = which.min(modelAICs) returnList[[1]] <- modelList[[minPos]] returnList[[2]] <- my.varList[[minPos]] cat(my.varList[[minPos]], ' selected for ', my.respVar, '\n') return(returnList) } } doLRTs <- function(my.coefsToTest, my.respVar, my.fullModel, my.varForm, my.d) { my.pVals <- data.frame(response = my.respVar) my.pVals[my.coefsToTest] <- NA if (my.varForm == "NA") { my.fullModelFit <- do.call('gls', list(model=as.formula(paste(my.respVar, "~", my.fullModel, collapse=" ")), data = substitute(my.d), na.action = na.omit, method="ML")) } else { my.fullModelFit <- do.call('gls', list(model=as.formula(paste(my.respVar, "~", my.fullModel, collapse=" ")), weights = substitute(eval(parse(text=eval(my.varForm)))), data = substitute(my.d), na.action = na.omit, method="ML")) } for (intVNum in 1:length(my.coefsToTest)) { my.redModelFit <- update(my.fullModelFit, as.formula(paste(".~.-",my.coefsToTest[[intVNum]],sep=""))) my.LRT <- anova(my.fullModelFit, my.redModelFit) my.pVals[my.coefsToTest[[intVNum]]] <- my.LRT$'p-value'[2] } return(my.pVals) } to.upper<-function(X) t(X)[lower.tri(X,diag=FALSE)] save.image("AirKnife_Data_R_Workspace_11.30.18.RData") write.xlsx(AKDAT, "AirKnife_Data_Small_11.30.18.xlsx") # filenameID <- as.vector(CLUSTDAT$FN) # technique <- as.vector(CLUSTDAT$Technique) # writeMat("TechniqueClustersForMatlab.mat", filenameID = filenameID, technique = technique)