---
title: "Code_RHP"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Load libraries
```{r, message = FALSE}
library(tidyverse)
library(lme4)
library(broom)
library(broom.mixed)
library(emmeans)
library(lmerTest)
library(infer)
library(dabestr)
library(emmeans)
library(MASS)
library(factoextra)
library(ggalt)
library(randomForest)
library(ggpubr)
library(cowplot)
library(boot)
```

# Data Preparation and Functions
Load relevant datasets
```{r, message = FALSE}
#Full database on all captured hummingbirds
captures <- read_csv("capswithnames.csv")

#Genetic sex information on all birds
sexes <- read_csv("geneticsexingcompiled.csv")

#Provides all individual ID info -- color tags, leg bands, RFID, mount name, etc
nameswide <- read_csv("names.csv")

#This is all the lifting data for all years
weightswide <- read_csv("weightdata.csv")

#This has the weights of each level on the chain
weightlevels <- read_csv("weightchain.csv")

#This has all the detailed wing measures
area <- read_csv("WingMeasurementsJJF.csv")
```

Gather the names database to make it long format, better for joins
```{r}
names <- nameswide %>%
  gather(idtype, ids, c(band, 
                        rfid, 
                        shortrfid, 
                        secondrfid, 
                        secondshortrfid, 
                        secondband), 
         na.rm = TRUE)
```

Prepare the capture dataset to join plumage and sex information
```{r}
caps <- captures %>%
  filter(species == "fmel") %>%
  dplyr::select(name, 
         date_captured, 
         plutype, 
         mass, 
         corrugations, 
         unusualplu, 
         uplu_orangeface,
         wingchord,
         wingdepth) %>%
  mutate(date_captured = lubridate::dmy(date_captured)) %>%
  arrange(name, date_captured) %>%
  left_join(sexes, by = "name") %>%
  dplyr::select(-extractionID, band, rfid, uniquerfid, extractionyear) %>%
  mutate(type = if_else(
    gensex == "m" & plutype == "m", "mlm", if_else(
      gensex == "f" & plutype == "f", "flf", if_else(
        gensex == "f" & plutype == "m", "mlf", NULL))))
```

Makes a confusion matrix for future analyses
```{r}
confusion <- function(actual, predicted, names = NULL, printit = TRUE, prior = NULL) {
  # names and priors are null unless otherwise specified; R should print our
  # output
  if (is.null(names)) 
    names <- levels(actual)
  tab <- table(actual, predicted)
  acctab <- t(apply(tab, 1, function(x) x/sum(x)))
  dimnames(acctab) <- list(Actual = names, Predicted = names)
  if (is.null(prior)) {
    relnum <- table(actual)
    prior <- relnum/sum(relnum)
    acc <- sum(tab[row(tab) == col(tab)])/sum(tab)
  } else {
    acc <- sum(prior * diag(acctab))
    names(prior) <- names
  }
  if (printit) 
    print(round(c(`Overall accuracy` = acc, `Prior frequency` = prior), 
                4))
  if (printit) {
    cat("\nConfusion matrix", "\n")
    print(round(acctab, 4))
  }
  invisible(acctab)
}
```

Use this to find mean and bootstrapped 95% CI 
```{r}
# get mean and 95% CI of values x via bootstrapping
boot_ci <- function(x, perms=5000, bca=F) {
  get_mean <- function(x, d) {
    return(mean(x[d]))
  }
  x <- as.vector(na.omit(x))
  mean <- mean(x)
  if(bca){
    boot <- boot.ci(boot(data=x, 
                         statistic=get_mean, 
                         R=perms, 
                         parallel = "multicore", 
                         ncpus = 4), 
                    type="bca")
    low <- boot$bca[1,4]
    high <- boot$bca[1,5] 
  }else{
    boot <- boot.ci(boot(data=x, 
                         statistic=get_mean, 
                         R=perms, 
                         parallel = "multicore", 
                         ncpus = 4), 
                    type="perc")
    low <- boot$perc[1,4]
    high <- boot$perc[1,5] 
  }
  c(low=low,mean=mean,high=high, N=round(length(x)))
}

#This function takes data from a permutation test and puts it into a histograph with p values
# plot permutation test results
hist_perm <- function(exp=exp, obs=obs, perms=perms, label=''){
  exp.range <- round(quantile(exp, probs= c(0.025, 0.975)),3)
  ggplot()+
    geom_histogram(aes(x=exp), color="black",fill="light blue")+
    geom_vline(aes(xintercept=obs), color="red", size=1)+
    xlab("expected values from null model")+
    ggtitle(label, subtitle = paste('obs = ',round(obs,3), ', exp = ', exp.range[1], ' to ', exp.range[2], ", Prob exp >= obs: p", ifelse(mean(exp>=obs)==0,paste("<",1/perms), paste("=",signif(mean(exp>=obs),digits=2))),", permutations=",perms, sep=""))
}
```

# Muscle Capacity Analysis
## Data Prep
Select the relevant data from the weights dataset and join it with names
```{r}
weightnames <- weightswide %>%
  filter(!is.na(date)) %>%
  filter(species == "flo") %>%
  dplyr::select(date, newdate, species, ID, '1', '2', '3', '4', '5', '6', '7', '8', '9', '10') %>%
  left_join(names, by = c("ID" = "ids")) %>%
  dplyr::select(-Mountname, -idtype) %>%
  mutate(date = lubridate::dmy(date)) %>%
  filter(!is.na(name))
```

Some experiment dates were unaccounted for in the capture data set. I went through and fixed the dates by using the dates on the videos and put these into the "newdate" column. Need to make another date column with all the correct dates
```{r}
#Creates a dataframe that gives unique name to each capture "namedate" in the capture dataset
capsx <- caps %>%
  dplyr::select(date_captured, name) %>%
  unite(namedate, c(name, date_captured))

#Creates a dataframe that gives unique name to each weight experiment "namedate" in the weight dataset.
#Then makes a column that checks whether the namedate in weights is in the capture dataset
weightx <- weightnames %>%
  dplyr::select(date, name) %>%
  unite(namedate, c(name, date)) %>%
  mutate(exist = namedate %in% capsx$namedate)

#List all the experiment dates that are not in the capture dataset
notincaps <- weightx %>% filter(!exist)

#Corrected dates are now in the column "newdate" 

#If there is something in newdate, use that date, if not, use the original date
weightdate <- weightnames %>% 
  mutate(date_correct = if_else(!is.na(newdate), newdate, date)) %>% 
  dplyr::select(-date, -newdate)

#Now check if there are any more individuals that don't have a date in the capture set
weightx2 <- weightdate %>%
  dplyr::select(date_correct, name) %>%
  unite(namedate, c(name, date_correct)) %>%
  mutate(exist = namedate %in% capsx$namedate)

#The only one that is not accounted for is Nicole on 2018-11-30. This date is not in the capture data for Nicole so it should be removed
weightx2 %>% filter(!exist)

```

Now we can match the mass on the exact date of the experiment
```{r}

#First prepare the sex datasex. Just take name and genetic sex
sex <- sexes %>% dplyr::select(name, gensex)

#Join it with the capture data. Bring together name and date to make an individual ID for each capture
capsex <- caps %>%
  unite(namedate, c(name, date_captured))

#Join with the capture data with the weight lifting data using namedate
weightcaps <- weightdate %>%
  unite(namedate, c(name, date_correct)) %>%
  inner_join(capsex, by = "namedate")
```

Juvenile cutoff 
```{r}
cutoff <- 40
```

Now prep the data and put it in the right format if we want to use the mean lift to do the analysis
```{r}
masslevels <- weightlevels %>% rename(beadmass = mass)

weightdata_means <- weightcaps %>%
  #First make all the lift data into long format
  pivot_longer(3:12, names_to = "liftnum", values_to = "beads") %>% 
  #Join the actual mass of the beads
  left_join(masslevels, by = c("beads" = "level")) %>%
  #Make the number of beads a factor
  mutate(liftnum = as.numeric(liftnum)) %>%
  #Number of lifts goes all the way to 10, but if the bird didn't lift that many it leaves an NA, so take those NAs out
  filter(!is.na(beadmass)) %>%
  #Group everything by the experiment namedate and summarize the mean and number of lifts, plus keep all the rest of the data
  group_by(namedate) %>%
  summarize(meanlift = mean(beadmass),
            lifts = n(),
            gensex = gensex[1], 
            plutype = plutype[1], 
            mass = mass[1], 
            corrugations = corrugations[1], 
            unusualplu = unusualplu[1], 
            uplu_orangeface = uplu_orangeface[1],
            wingchord = wingchord[1]) %>%
  #Take out any individuals that aren't sexed
  filter(!is.na(gensex)) %>%
  #Take out the mix plumage types
  filter(plutype != "mix") %>%
  #Divide up the namedate column back into two seperate name and date collumns
  mutate(name = sapply(strsplit(namedate, split = "_", fixed = TRUE), function(x) (x[1]))) %>%
  mutate(date = sapply(strsplit(namedate, split = "_", fixed = TRUE), function(x) (x[2]))) %>%
  #Change the date back into a date type
  mutate(date = lubridate::ymd(date)) %>%
  #Only keep experiments where there are at least 3 lifts
  filter(lifts > 3) %>%
  #Add a column for lift total (including bird mass), and the amount lifted per gram of body mass
  mutate(totlift = mass + meanlift,
         pergram = meanlift/mass) %>%
  arrange(name, date) %>%
  distinct(name, .keep_all = TRUE) %>%
  mutate(ageclass = if_else(corrugations <= cutoff, "adult", "juv")) %>%
  mutate(type = if_else(
    gensex == "m" & plutype == "m", "mlm", if_else(
      gensex == "f" & plutype == "f", "flf", if_else(
        gensex == "f" & plutype == "m", "mlf", NULL))))

#Just adults
weightmeans.adult <- weightdata_means %>% filter(corrugations <= cutoff)
```

## Analysis

We can do a quick regression to see if adults and juveniles are different in their muscle capacity
```{r}
agemodel <- lm(totlift ~ ageclass*gensex, data = weightdata_means)
emmeans(agemodel, pairwise ~ ageclass|gensex, type = "response")
```

For most of the analyses, we'll only pay attention to the adults

Get the mean and 95% CI for each type
```{r}
#Separate the heterochrome females
weightmeans.adultflf <- weightmeans.adult %>% filter(type == "flf")
#Separate the androchrome females
weightmeans.adultmlf <- weightmeans.adult %>% filter(type == "mlf")
#Separate the androchrome males
weightmeans.adultmlm <- weightmeans.adult %>% filter(type == "mlm")

#Get mean and 95% CI for each group
flfad <- boot_ci(weightmeans.adultflf$totlift)
#      low     mean     high        N 
# 17.76265 18.43948 19.12429 52.00000 
mlfad <- boot_ci(weightmeans.adultmlf$totlift)
#      low     mean     high        N 
# 17.87107 19.07117 20.20099 10.00000 
mlmad <- boot_ci(weightmeans.adultmlm$totlift)
 #      low      mean      high         N 
 # 21.52853  22.00946  22.48610 105.00000 
```

Get p-values with permutation tests
```{r}
#How many of each type?
table(weightmeans.adult$type)

#What's the observed difference between each type?
diffsobs_muscle <- weightmeans.adult %>%
  group_by(type) %>%
  summarize(mean = mean(totlift)) %>%
  pivot_wider(names_from = type, values_from = mean) %>%
  mutate(flfmlf = flf-mlf,
         flfmlm = flf-mlm,
         mlfmlm = mlf-mlm) %>% data.frame()

#Now permute the bird types and find the differences between each type in each permutation
perms <- 5000
exp <- matrix(NA, ncol=3, nrow=perms)
set.seed(123)

for(i in 1:perms){
  diffperm <- weightmeans.adult %>% mutate(perm = sample(type))
  diffpermcalc <- diffperm %>%
    group_by(perm) %>%
    summarize(perm_totlift = mean(totlift)) %>%
    pivot_wider(names_from = perm, values_from = perm_totlift) %>%
    mutate(flfmlf = flf-mlf,
           flfmlm = flf-mlm,
           mlfmlm = mlf-mlm) %>%
    data.frame()
  exp[i,] <- c(diffpermcalc[1,4], diffpermcalc[1,5], diffpermcalc[1,6])
}
```

Use the permuted datasets to calculate p values 
```{r}
#diff between flf-mlf
hist_perm(exp= exp[,1], obs= diffsobs_muscle[,4], perms=perms)
pflfmlf <- ifelse(mean(exp[,1]>=diffsobs_muscle[,4])==0, 
                  1/perms, signif(mean(exp[,1]>=diffsobs_muscle[,4]),digits=2))

#diff between flf-mlm
hist_perm(exp= exp[,2], obs= diffsobs_muscle[,5], perms=perms)
pflfmlm <- ifelse(mean(exp[,1]>=diffsobs_muscle[,5])==0, 
                  1/perms, signif(mean(exp[,1]>=diffsobs_muscle[,5]),digits=2))

#diff between mlf-mlm
hist_perm(exp= exp[,3], obs= diffsobs_muscle[,6], perms=perms)
pmlfmlm <- ifelse(mean(exp>=diffsobs_muscle[,6])==0, 
                  1/perms, signif(mean(exp[,1]>=diffsobs_muscle[,6]),digits=2))

musclead_p <- c(pflfmlf, pflfmlm, pmlfmlm)
```

Make a graph
```{r}
meansCI.ad <- tibble(type = c("flf", "mlf", "mlm"),
                     mean = c(flfad[[2]], mlfad[[2]], mlmad[[2]]),
                     LCL = c(flfad[[1]], mlfad[[1]], mlmad[[1]]),
                     UCL = c(flfad[[3]], mlfad[[3]], mlmad[[3]]))

musclead_contrast_p <- tibble(contrast = c("mlf-flf", "flf-mlm", "mlf-mlm"),
                              p = c(musclead_p[1],musclead_p[2],musclead_p[3]),
                              group1 = c("flf", "flf", "mlf"),
                              group2 = c("mlf", "mlm", "mlm"), 
                              p.adj.signif = if_else(p > 0.9998, "****",
                                if_else(p > 0.999, "***",
                                        if_else(p > 0.99, "**",
                                                if_else(p > 0.975, "*", "N.S.")))))

muscleplotad <- ggplot(data = weightmeans.adult) +
  theme_classic() +
  geom_violin(aes(x = type, y = totlift, fill = type), trim = TRUE, scale = "count") +
  scale_fill_manual(values = c("darkseagreen1", "lightsteelblue1", "gray52")) +
  theme(legend.position = "none",
        axis.ticks.x = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 8),
        axis.text.y = element_text(size = 6)) +
  geom_errorbar(data = meansCI.ad, aes(x = type, ymin = LCL, ymax = UCL), color = c("seagreen", "steelblue3", "gray3"), width = 0) +
  geom_point(data = meansCI.ad, aes(x = type, y = mean), size = 0.7, stroke = 1, color = c("seagreen", "steelblue3",  "gray3")) +
  ylab("Total lift capacity (g)") + 
  scale_x_discrete(labels = NULL) +
  stat_pvalue_manual(musclead_contrast_p,
                     tip.length = 0,
                     label = "p.adj.signif", 
                     y.position = 26,
                     step.increase = 0.07,
                     size = 2.5)

muscleplotad
```

# Wing Load Analysis
## Data Prep
```{r}
wingload <- capsex %>%
  dplyr::select(namedate, plutype, gensex, mass, corrugations, wingchord, wingdepth) %>%
  mutate(name = sapply(strsplit(namedate, split = "_", fixed = TRUE), function(x) (x[1]))) %>%
  mutate(date = sapply(strsplit(namedate, split = "_", fixed = TRUE), function(x) (x[2]))) %>%
  distinct(name, .keep_all = TRUE) %>%
  mutate(wingchordm = wingchord / 1000,
         wingdepthm = wingdepth / 1000,
         wingloading = (mass/1000) / (wingchordm * wingdepthm),
         type = if_else(gensex == "m" & plutype == "m", "mlm", if_else(
           gensex == "f" & plutype == "f", "flf", if_else(
             gensex == "f" & plutype == "m", "mlf", NULL)))) %>%
  filter(!is.na(type)) %>%
  filter(!is.na(wingloading))
```

Remove juveniles
```{r}
cutoff <- 40
load.adult <- wingload %>% filter(corrugations <= cutoff)
```

## Analysis

Mean and bootstrapped 95%CI
```{r}
#Separate the heterochrome females
loadad_flf <- load.adult %>% filter(type == "flf")
#Separate the androchrome females
loadad_mlf <- load.adult %>% filter(type == "mlf")
#Separate the androchrome males
loadad_mlm <- load.adult %>% filter(type == "mlm")

#Calculate the mean and 95% CI for each group
loadflf <- boot_ci(loadad_flf$wingloading)
 #      low      mean      high         N 
 # 4.620477  4.715493  4.815115 78.000000 
loadmlf <- boot_ci(loadad_mlf$wingloading)
 #      low      mean      high         N 
 # 4.504915  4.663932  4.833966 20.000000 
loadmlm <- boot_ci(loadad_mlm$wingloading)
  #     low       mean       high          N 
  # 4.504150   4.558456   4.616044 184.000000 

meansCI.load <- tibble(type = c("flf", "mlf", "mlm"),
                     mean = c(loadflf[[2]], loadmlf[[2]], loadmlm[[2]]),
                     LCL = c(loadflf[[1]], loadmlf[[1]], loadmlm[[1]]),
                     UCL = c(loadflf[[3]], loadmlf[[3]], loadmlm[[3]]))
```

Get p-values with a permutation test
```{r}
#How many adults of each type?
table(load.adult$type)

#What's the observed difference in wing loading between each type? 
diffsobs_load <- load.adult %>%
  group_by(type) %>%
  summarize(mean = mean(wingloading)) %>%
  pivot_wider(names_from = type, values_from = mean) %>%
  mutate(flfmlf = flf-mlf,
         flfmlm = flf-mlm,
         mlfmlm = mlf-mlm) %>% data.frame()

#Now permute the bird types and recalculate the difference 5000 times to calculate p values
perms <- 5000
exp <- matrix(NA, ncol=3, nrow=perms)
set.seed(123)

for(i in 1:perms){
  diffperm <- load.adult %>% mutate(perm = sample(type))
  diffpermcalc <- diffperm %>%
    group_by(perm) %>%
    summarize(perm_wingload = mean(wingloading)) %>%
    pivot_wider(names_from = perm, values_from = perm_wingload) %>%
    mutate(flfmlf = flf-mlf,
           flfmlm = flf-mlm,
           mlfmlm = mlf-mlm) %>%
    data.frame()
  exp[i,] <- c(diffpermcalc[1,4], diffpermcalc[1,5], diffpermcalc[1,6])
}
```

Compare the observed with permuted differences
```{r}
#diff between flf-mlf
hist_perm(exp= exp[,1], obs= diffsobs_load[,4], perms=perms)
pflfmlf <- ifelse(mean(exp[,1]>=diffsobs_load[,4])==0, 
                  1/perms, signif(mean(exp[,1]>=diffsobs_load[,4]),digits=2))

#diff between flf-mlm
hist_perm(exp= exp[,2], obs= diffsobs_load[,5], perms=perms)
pflfmlm <- ifelse(mean(exp[,2]>=diffsobs_load[,5])==0, 
                  1/perms, signif(mean(exp[,2]>=diffsobs_load[,5]),digits=2))

#diff between mlf-mlm
hist_perm(exp= exp[,3], obs= diffsobs_load[,6], perms=perms)
pmlfmlm <- ifelse(mean(exp[,3]>=diffsobs_load[,6])==0, 
                  1/perms, signif(mean(exp[,3]>=diffsobs_load[,6]),digits=2))

load_p <- c(pflfmlf, pflfmlm, pmlfmlm)
```

Make a graph
```{r}

load_contrast_p <- tibble(contrast = c("flf-mlf", "flf-mlm", "mlf-mlm"),
                              p = c(load_p[1],load_p[2],load_p[3]),
                              group1 = c("flf", "flf", "mlf"),
                              group2 = c("mlf", "mlm", "mlm"), 
                              p.adj.signif = if_else(p < 0.0002, "****",
                                if_else(p < 0.001, "***",
                                        if_else(p < 0.01, "**",
                                                if_else(p < 0.025, "*", "N.S.")))))

loadplotad <- ggplot(data = load.adult) +
  theme_classic() +
  geom_violin(aes(x = type, y = wingloading, fill = type), trim = TRUE, scale = "count") +
  scale_fill_manual(values = c("darkseagreen1", "lightsteelblue1", "gray52")) +
  geom_errorbar(data = meansCI.load, aes(x = type, ymin = LCL, ymax = UCL), color = c("seagreen", "steelblue3", "gray3"), width = 0) +
  geom_point(data = meansCI.load, aes(x = type, y = mean), size = 0.7, stroke = 1, color = c("seagreen", "steelblue3",  "gray3")) +
  ylab(expression(Wing~Load~(kg/m^{"2"}))) + 
  scale_x_discrete(labels = c("Speckle-throated\n(Heterochrome)\nFemale", "Blue-headed\n(Androchrome)\nFemale", "Blue-headed\n(Androchrome)\nMale")) +
  ylim(2.8, 7.5) +
  stat_pvalue_manual(load_contrast_p,
                     tip.length = 0,
                     label = "p.adj.signif", 
                     y.position = 6.5,
                     step.increase = 0.1, 
                     size = 2.5) +
  theme(legend.position = "none",
        axis.ticks.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_text(size = 8),
        axis.title.y = element_text(size = 8),
        axis.text.y = element_text(size = 6)) 

loadplotad
```
Put the muscle-capacity and wing-load graphs together
```{r}
muscleandload <- plot_grid(muscleplotad, loadplotad, labels = c("A","B"), ncol = 1, align = "v")
muscleandload 
#ggsave("muscleandload.png", muscleandload, width = 4, height = 5)
```

# Wing Shape Measures
## Data Preparation
```{r}
sexesprep <- sexes %>% dplyr::select(name, gensex)

#join the genetic sex data with the capture data to get plumage types
capsprep <- captures %>%
  dplyr::select(name, species,date_captured, plutype, mass, culmen, nare, tarsus, wingchord, wingdepth, tail, corrugations, uplu_orangeface) %>%
  inner_join(sexesprep) %>%
  mutate(day = lubridate::dmy(date_captured)) %>%
  mutate(namedate = str_c("wing_",day,"_",name)) %>%
  mutate(type = if_else(
    gensex == "m" & plutype == "m", "mlm", if_else(
      gensex == "f" & plutype == "f", "flf", if_else(
        gensex == "f" & plutype == "m", "mlf", "mix"))))

#I'm lazy and the names of the columns are too long! I made them shorter. 
areaprep <- area %>%
  rename(wingarea = `Wing_area_70X_mm^2`, 
         wingarea_detail = Wing_area_Every_detail_mm2,
         imageID = `Image ID`,
         length = `R_ Wing_Length_mm`, 
         width = Wing_width_mm,
         shaperatio = `Rs_Shape ratio`,
         aspectratio = Ra_Aspect_ratio,
         taper = `Tw_Wing taper`) %>%
  dplyr::select(imageID, 
         Wing_spread_properly, 
         Molting, 
         finger_hold,
         length,
         width,
         shaperatio,
         wingarea,
         wingarea_detail,
         aspectratio,
         taper)

#bring together all the capture data with the wing data and make some calculations
areaready <- areaprep %>%
  inner_join(capsprep, by = c("imageID" = "namedate")) %>%
  filter(species == "fmel") %>%
  filter(Wing_spread_properly == "Yes") %>%
  filter(Molting == "No") %>%
  filter(!is.na(gensex)) %>%
  mutate(wingload = mass/(wingarea*2)) %>%
  mutate(wingloadcalc = mass/(wingchord*wingdepth)) %>%
  mutate(areacalc = wingchord*wingdepth/2) %>%
  replace_na(list(uplu_orangeface = 0)) %>%
  mutate(uplu_orangeface = as.factor(uplu_orangeface)) %>%
  mutate(cuberootmass = mass^(1/3)) %>%
  mutate(lengthpermass = length/cuberootmass) %>%
  mutate(widthpermass = width/cuberootmass) %>%
  mutate(masscm = mass/10) %>%
  filter(finger_hold == "pull")

#Gather the names database to make it long format, better for joins
names <- nameswide %>%
  gather(idtype, ids, c(band, 
                        rfid, 
                        shortrfid, 
                        secondrfid, 
                        secondshortrfid, 
                        secondband), 
         na.rm = TRUE)

ad_areaready <- areaready %>% filter(corrugations <= 40)
```

## Analysis

How well do the estimated wing areas correlate with measured ones?
```{r}
ggplot(data = areaready, aes(x = wingarea, y = areacalc)) +
  geom_point(aes(col = gensex)) +
  geom_smooth(method = lm)

areacalc.mod <- lm(areacalc ~ wingarea, data = areaready)
summary(areacalc.mod)

cor(areaready$areacalc, areaready$wingarea)
```

To keep things consistent, let's just look at the adults.  
```{r}
mlmwing <- areaready %>% filter(type == "mlm") %>% filter(corrugations <= 40)
flfwing <- areaready %>% filter(type == "flf") %>% filter(corrugations <= 40)
mlfwing <- areaready %>% filter(type == "mlf") %>% filter(corrugations <= 40)
```

### Wing Area
Calculate mean and 95% CI
```{r}
boot_ci(flfwing$wingarea)
#      low     mean     high        N 
# 1133.971 1160.964 1185.804   31.000  
boot_ci(mlfwing$wingarea)
#      low     mean     high        N 
# 1152.005 1194.294 1244.547    8.000 
boot_ci(mlmwing$wingarea)
#      low     mean     high        N 
# 1251.354 1270.184 1289.625   77.000 
```

Get the observed and permuted differences
```{r}
diffsobs_wingarea <- ad_areaready %>%
  group_by(type) %>%
  summarize(mean = mean(wingarea)) %>%
  pivot_wider(names_from = type, values_from = mean) %>%
  mutate(flfmlf = flf-mlf,
         flfmlm = flf-mlm,
         mlfmlm = mlf-mlm) %>% data.frame()

perms <- 5000
exp <- matrix(NA, ncol=3, nrow=perms)
set.seed(123)

for(i in 1:perms){
  diffperm <- ad_areaready %>% mutate(perm = sample(type))
  diffpermcalc <- diffperm %>%
    group_by(perm) %>%
    summarize(perm_wingarea = mean(wingarea)) %>%
    pivot_wider(names_from = perm, values_from = perm_wingarea) %>%
    mutate(flfmlf = flf-mlf,
           flfmlm = flf-mlm,
           mlfmlm = mlf-mlm) %>%
    data.frame()
  exp[i,] <- c(diffpermcalc[1,4], diffpermcalc[1,5], diffpermcalc[1,6])
}
```

Extract the p values
```{r}
#diff between flf-mlf
hist_perm(exp= exp[,1], obs= diffsobs_wingarea[,4], perms=perms)

#diff between flf-mlm
hist_perm(exp= exp[,2], obs= diffsobs_wingarea[,5], perms=perms)

#diff between mlf-mlm
hist_perm(exp= exp[,3], obs= diffsobs_wingarea[,6], perms=perms)
```

### Wing Length
Calculate mean and 95% CI
```{r}
boot_ci(flfwing$length)
#      low     mean     high        N 
# 68.87484 69.72610 70.57230 31.00000 
boot_ci(mlfwing$length)
#      low     mean     high        N 
# 69.06354 70.35438 71.69270  8.00000 
boot_ci(mlmwing$length)
#      low     mean     high        N 
# 72.70533 73.29773 73.90591 77.00000 
```

Get the observed and permuted differences
```{r}
diffsobs_length <- ad_areaready %>%
  group_by(type) %>%
  summarize(mean = mean(length)) %>%
  pivot_wider(names_from = type, values_from = mean) %>%
  mutate(flfmlf = flf-mlf,
         flfmlm = flf-mlm,
         mlfmlm = mlf-mlm) %>% data.frame()

perms <- 5000
exp <- matrix(NA, ncol=3, nrow=perms)
set.seed(123)

for(i in 1:perms){
  diffperm <- ad_areaready %>% mutate(perm = sample(type))
  diffpermcalc <- diffperm %>%
    group_by(perm) %>%
    summarize(perm_length = mean(length)) %>%
    pivot_wider(names_from = perm, values_from = perm_length) %>%
    mutate(flfmlf = flf-mlf,
           flfmlm = flf-mlm,
           mlfmlm = mlf-mlm) %>%
    data.frame()
  exp[i,] <- c(diffpermcalc[1,4], diffpermcalc[1,5], diffpermcalc[1,6])
}
```

Extract the p values
```{r}
#diff between flf-mlf
hist_perm(exp= exp[,1], obs= diffsobs_length[,4], perms=perms)

#diff between flf-mlm
hist_perm(exp= exp[,2], obs= diffsobs_length[,5], perms=perms)

#diff between mlf-mlm
hist_perm(exp= exp[,3], obs= diffsobs_length[,6], perms=perms)
```

### Wing Width
Calculate mean and 95% CI
```{r}
boot_ci(flfwing$width)
#      low     mean     high        N 
# 22.18010 22.46635 22.75550 31.00000 
boot_ci(mlfwing$width)
#      low     mean     high        N 
# 22.65800 22.97525 23.21336  8.00000 
boot_ci(mlmwing$width)
#      low     mean     high        N 
# 24.43319 24.67655 24.90649 77.00000 
```

Get the observed and permuted differences
```{r}
diffsobs_width <- ad_areaready %>%
  group_by(type) %>%
  summarize(mean = mean(width)) %>%
  pivot_wider(names_from = type, values_from = mean) %>%
  mutate(flfmlf = flf-mlf,
         flfmlm = flf-mlm,
         mlfmlm = mlf-mlm) %>% data.frame()

perms <- 5000
exp <- matrix(NA, ncol=3, nrow=perms)
set.seed(123)

for(i in 1:perms){
  diffperm <- ad_areaready %>% mutate(perm = sample(type))
  diffpermcalc <- diffperm %>%
    group_by(perm) %>%
    summarize(perm_width = mean(width)) %>%
    pivot_wider(names_from = perm, values_from = perm_width) %>%
    mutate(flfmlf = flf-mlf,
           flfmlm = flf-mlm,
           mlfmlm = mlf-mlm) %>%
    data.frame()
  exp[i,] <- c(diffpermcalc[1,4], diffpermcalc[1,5], diffpermcalc[1,6])
}
```

Extract the p values
```{r}
#diff between flf-mlf
hist_perm(exp= exp[,1], obs= diffsobs_width[,4], perms=perms)

#diff between flf-mlm
hist_perm(exp= exp[,2], obs= diffsobs_width[,5], perms=perms)

#diff between mlf-mlm
hist_perm(exp= exp[,3], obs= diffsobs_width[,6], perms=perms)
```

### Wing Load
Calculate mean and 95% CI
```{r}
boot_ci(flfwing$wingload)
 #         low         mean         high            N 
 # 0.002851424  0.002955644  0.003072278 31.000000000 
boot_ci(mlfwing$wingload)
#          low        mean        high           N 
# 0.002680717 0.002783887 0.002886628 8.000000000 
boot_ci(mlmwing$wingload)
 #         low         mean         high            N 
 # 0.002890085  0.002961759  0.003036603 77.000000000
```

Get the observed and permuted differences
```{r}
diffsobs_wingload <- ad_areaready %>%
  group_by(type) %>%
  summarize(mean = mean(wingload)) %>%
  pivot_wider(names_from = type, values_from = mean) %>%
  mutate(flfmlf = flf-mlf,
         flfmlm = flf-mlm,
         mlfmlm = mlf-mlm) %>% data.frame()

perms <- 5000
exp <- matrix(NA, ncol=3, nrow=perms)
set.seed(123)

for(i in 1:perms){
  diffperm <- ad_areaready %>% mutate(perm = sample(type))
  diffpermcalc <- diffperm %>%
    group_by(perm) %>%
    summarize(perm_wingload = mean(wingload)) %>%
    pivot_wider(names_from = perm, values_from = perm_wingload) %>%
    mutate(flfmlf = flf-mlf,
           flfmlm = flf-mlm,
           mlfmlm = mlf-mlm) %>%
    data.frame()
  exp[i,] <- c(diffpermcalc[1,4], diffpermcalc[1,5], diffpermcalc[1,6])
}
```

Extract the p values
```{r}
#diff between flf-mlf
hist_perm(exp= exp[,1], obs= diffsobs_wingload[,4], perms=perms)

#diff between flf-mlm
hist_perm(exp= exp[,2], obs= diffsobs_wingload[,5], perms=perms)

#diff between mlf-mlm
hist_perm(exp= exp[,3], obs= diffsobs_wingload[,6], perms=perms)
```

### Aspect Ratio
Calculate mean and 95% CI
```{r}
boot_ci(flfwing$aspectratio)
 #      low      mean      high         N 
 # 8.448839  8.531161  8.624290 31.000000 
boot_ci(mlfwing$aspectratio)
#      low     mean     high        N 
# 8.295138 8.466500 8.618250 8.000000 
boot_ci(mlmwing$aspectratio)
 #       low      mean      high         N 
 # 8.542298  8.623766  8.706388 77.000000 
```

Get the observed and permuted differences
```{r}
diffsobs_aspectratio <- ad_areaready %>%
  group_by(type) %>%
  summarize(mean = mean(aspectratio)) %>%
  pivot_wider(names_from = type, values_from = mean) %>%
  mutate(flfmlf = flf-mlf,
         flfmlm = flf-mlm,
         mlfmlm = mlf-mlm) %>% data.frame()

perms <- 5000
exp <- matrix(NA, ncol=3, nrow=perms)
set.seed(123)

for(i in 1:perms){
  diffperm <- ad_areaready %>% mutate(perm = sample(type))
  diffpermcalc <- diffperm %>%
    group_by(perm) %>%
    summarize(perm_aspectratio = mean(aspectratio)) %>%
    pivot_wider(names_from = perm, values_from = perm_aspectratio) %>%
    mutate(flfmlf = flf-mlf,
           flfmlm = flf-mlm,
           mlfmlm = mlf-mlm) %>%
    data.frame()
  exp[i,] <- c(diffpermcalc[1,4], diffpermcalc[1,5], diffpermcalc[1,6])
}
```

Extract the p values
```{r}
#diff between flf-mlf
hist_perm(exp= exp[,1], obs= diffsobs_aspectratio[,4], perms=perms)

#diff between flf-mlm
hist_perm(exp= exp[,2], obs= diffsobs_aspectratio[,5], perms=perms)

#diff between mlf-mlm
hist_perm(exp= exp[,3], obs= diffsobs_aspectratio[,6], perms=perms)
```

### Shape Ratio
Calculate mean and 95% CI
```{r}
boot_ci(flfwing$shaperatio)
 #      low      mean      high         N 
 # 3.088871  3.133000  3.178000 31.000000   
boot_ci(mlfwing$shaperatio)
#       low     mean     high        N 
# 3.014753 3.094250 3.169250 8.000000 
boot_ci(mlmwing$shaperatio)
 #       low      mean      high         N 
 # 2.965610  3.000974  3.035843 77.000000 
```

Get the observed and permuted differences
```{r}
diffsobs_shaperatio <- ad_areaready %>%
  group_by(type) %>%
  summarize(mean = mean(shaperatio)) %>%
  pivot_wider(names_from = type, values_from = mean) %>%
  mutate(flfmlf = flf-mlf,
         flfmlm = flf-mlm,
         mlfmlm = mlf-mlm) %>% data.frame()

perms <- 5000
exp <- matrix(NA, ncol=3, nrow=perms)
set.seed(123)

for(i in 1:perms){
  diffperm <- ad_areaready %>% mutate(perm = sample(type))
  diffpermcalc <- diffperm %>%
    group_by(perm) %>%
    summarize(perm_shaperatio = mean(shaperatio)) %>%
    pivot_wider(names_from = perm, values_from = perm_shaperatio) %>%
    mutate(flfmlf = flf-mlf,
           flfmlm = flf-mlm,
           mlfmlm = mlf-mlm) %>%
    data.frame()
  exp[i,] <- c(diffpermcalc[1,4], diffpermcalc[1,5], diffpermcalc[1,6])
}
```

Extract the p values
```{r}
#diff between flf-mlf
hist_perm(exp= exp[,1], obs= diffsobs_shaperatio[,4], perms=perms)

#diff between flf-mlm
hist_perm(exp= exp[,2], obs= diffsobs_shaperatio[,5], perms=perms)

#diff between mlf-mlm
hist_perm(exp= exp[,3], obs= diffsobs_shaperatio[,6], perms=perms)
```

### Taper
Calculate mean and 95% CI
```{r}
boot_ci(flfwing$taper)
 #        low       mean       high          N 
 # 0.3443895  0.3629677  0.3830000 31.0000000   
boot_ci(mlfwing$taper)
#        low      mean      high         N 
# 0.3347500 0.3692500 0.4013719 8.0000000
boot_ci(mlmwing$taper)
 #       low       mean       high          N 
 # 0.4254688  0.4387532  0.4522727 77.0000000  
```

Get the observed and permuted differences
```{r}
diffsobs_taper <- ad_areaready %>%
  group_by(type) %>%
  summarize(mean = mean(taper)) %>%
  pivot_wider(names_from = type, values_from = mean) %>%
  mutate(flfmlf = flf-mlf,
         flfmlm = flf-mlm,
         mlfmlm = mlf-mlm) %>% data.frame()

perms <- 5000
exp <- matrix(NA, ncol=3, nrow=perms)
set.seed(123)

for(i in 1:perms){
  diffperm <- ad_areaready %>% mutate(perm = sample(type))
  diffpermcalc <- diffperm %>%
    group_by(perm) %>%
    summarize(perm_taper = mean(taper)) %>%
    pivot_wider(names_from = perm, values_from = perm_taper) %>%
    mutate(flfmlf = flf-mlf,
           flfmlm = flf-mlm,
           mlfmlm = mlf-mlm) %>%
    data.frame()
  exp[i,] <- c(diffpermcalc[1,4], diffpermcalc[1,5], diffpermcalc[1,6])
}
```

Extract the p values
```{r}
#diff between flf-mlf
hist_perm(exp= exp[,1], obs= diffsobs_taper[,4], perms=perms)

#diff between flf-mlm
hist_perm(exp= exp[,2], obs= diffsobs_taper[,5], perms=perms)

#diff between mlf-mlm
hist_perm(exp= exp[,3], obs= diffsobs_taper[,6], perms=perms)
```

# Body Size Analysis
```{r}
caps <- captures %>%
  dplyr::select(name, primary_molt, retrix_molt, date_captured, plutype, mass, corrugations, wingchord, wingdepth, tail, nare, culmen, tarsus) %>%
  arrange(name) %>%
  distinct(name, .keep_all = TRUE) %>%
  inner_join(sexes, by = "name") %>%
  mutate(type = if_else(gensex == "m" & plutype == "m", "mlm", 
                        if_else(gensex == "f" & plutype == "f", "flf", 
                                if_else(gensex == "f" & plutype == "m", "mlf", "mix")))) %>%
  dplyr::select(name, corrugations, primary_molt, retrix_molt, date_captured, type, plutype, mass, corrugations, wingchord, wingdepth, gensex, type, tail, nare, culmen, tarsus) %>%
  filter(!is.na(type)) %>%
  filter(type != "mix") %>%
  filter(!is.na(mass)) %>%
  filter(!is.na(wingchord)) %>%
  filter(!is.na(wingdepth)) %>%
  filter(!is.na(tail)) %>%
  filter(!is.na(nare)) %>%
  filter(!is.na(culmen)) %>%
  filter(!is.na(tarsus)) %>%
  filter(name != "Irving") #extremely short tail length which is probably from a molt
```

## PCA
```{r}
justmorph <- caps %>% dplyr::select(name, 
                             gensex,
                             type, 
                             corrugations,
                             mass,
                             wingchord,
                             wingdepth, 
                             tail,
                             nare,
                             culmen,
                             tarsus) %>%
  filter(!is.na(corrugations)) %>%
  mutate(age = if_else(corrugations <= 40, "ad", "im")) %>%
  mutate(typeage = str_c(type, age))

justmorphad <- justmorph %>% filter(age == "ad")

PC.morph <- prcomp(justmorphad[,c(-1,-2,-3,-4,-12,-13)], scale = TRUE)
```

Pull some of the info from the PC Analysis
```{r}
#Shows the loadings
PC.morph$rotation
#Shows the amount of variation explained by each PC axis
summary(PC.morph)
#Scree plot. Elbow might occur at 3 or 5 dimensions
fviz_eig(PC.morph)
#Only the first 2 dimensions have eigen value of 1 or greater. Covers 63.9% of variation together.
get_eigenvalue(PC.morph)
#Create a data frame with PC values and name, gensex, and date
pcas.morph <- data.frame(justmorphad[,c(1,2,3,4,12,13)],PC.morph$x) 

fviz_pca_var(PC.morph, col.var = "black")
```
 
Make a graph
```{r}
morphspace <- ggplot(data = pcas.morph, aes(x = PC1, y = PC2)) +
  labs(x = "Morphology PC1", y = "Morphology PC2") +
  geom_encircle(aes(fill = as.factor(type)), alpha = 0.3, expand = 0.04, s_shape = 0.8, show.legend = FALSE) +
  theme_classic() +
  theme(axis.text = element_text(size = 9),
        axis.title = element_text(size = 10)) +
  scale_fill_manual(values = c("seagreen3", "steelblue3","grey11")) +
  geom_point(aes(color = type), size = 1, show.legend = FALSE) +
  scale_color_manual(values = c("seagreen3", "steelblue3", "grey11"))

morphspace
```

## Comparisons with LDA permutations

### FLF-MLF LDA

Get the observed prediction rate
```{r}
mlfflf <- caps %>% filter(type %in% c("mlf","flf")) %>% filter(corrugations <= 40)

prop.table(table(mlfflf$type))
mlfflf.lda <- lda(type ~ mass + nare + culmen + tarsus + wingchord + wingdepth + tail, data = mlfflf, prior = c(0.796, 0.204), na.action = "na.omit")

mlfflfp.morph = predict(mlfflf.lda)

mlfflfknifed_bird <- lda(type ~ mass + nare + culmen + tarsus + wingchord + wingdepth + tail, data = mlfflf, prior = c(0.796, 0.204), na.action = "na.omit", CV=TRUE)
summary(mlfflfknifed_bird)
mlfflfknifed_bird$posterior[1:5, ]
confuse <- confusion(mlfflf$type, mlfflfknifed_bird$class, prior = c(0.796, 0.204))

observed <- confuse[1,1]*0.796 + confuse[2,2]*0.204
observed
```

Get p-values with a permutation test
```{r}
#set number of permutations
reps <- 5000

#create initial random data
rand.data <- transform(mlfflf, type = sample(type))

#create empty matrix for randomized effect sizes
coefs.rand <- matrix(NA, nrow = reps, ncol = 1)

#fill matrix with results of tests with randomized data (this runs the DFA model 5000 times)
set.seed(1)
for(i in 1:reps){
  rand.data <- transform(rand.data, type = sample(type)) #randomize the sexes
  rand.fit <- lda(type ~ mass + nare + culmen + tarsus + wingchord + wingdepth + tail, data = rand.data, prior = c(0.796, 0.204), na.action = "na.omit") #run DFA with randomized type data
  
  rand.fit.p <- predict(rand.fit) #generate predictions based on DFA
  rand.confuse <- confusion(rand.data$type, rand.fit.p$class, prior = c(0.796, 0.204)) #generate confusion matrix with prediction rates for each type
  
  #insert prediction rates into the matrix	
  coefs.rand[i] <- rand.confuse[1,1]*0.796 + rand.confuse[2,2]*0.204
}

expected <- data.frame(coefs.rand)

#calculate p-values
p <- sum(expected >= observed) / reps # positive effect
p
```

### MLM-MLF LDA

```{r}

mlfmlm <- caps %>% filter(type %in% c("mlf","mlm")) %>% filter(corrugations <= 40)

prop.table(table(mlfmlm$type))
mlfmlm.lda <- lda(type ~ mass + nare + culmen + tarsus + wingchord + wingdepth + tail, data = mlfmlm, prior = c(0.1, 0.9), na.action = "na.omit")

mlfmlm.morph = predict(mlfmlm.lda)

mlfmlmknifed_bird <- lda(type ~ mass + nare + culmen + tarsus + wingchord + wingdepth + tail, data = mlfmlm, prior = c(0.1, 0.9), na.action = "na.omit", CV = TRUE)
summary(mlfmlmknifed_bird)
mlfmlmknifed_bird$posterior[1:5, ]
confuse <- confusion(mlfmlm$type, mlfmlmknifed_bird$class, prior = c(0.1, 0.9))

observed <- confuse[1]*0.1 + confuse[2,2]*0.9
observed
```

```{r}
#set number of permutations
reps <- 5000

#create initial random data
rand.data <- transform(mlfmlm , type = sample(type))

#create empty matrix for randomized effect sizes
coefs.rand <- matrix(NA, nrow = reps, ncol = 1)

#fill matrix with results of tests with randomized data (this runs the DFA model 5000 times)
set.seed(1)
for(i in 1:reps){
  rand.data <- transform(rand.data, type = sample(type)) #randomize the sexes
  rand.fit <- lda(type ~ mass + nare + culmen + tarsus + wingchord + wingdepth + tail, data = rand.data, prior = c(.1, 0.9), na.action = "na.omit") #run DFA with randomized type data
  
  rand.fit.p <- predict(rand.fit) #generate predictions based on DFA
  rand.confuse <- confusion(rand.data$type, rand.fit.p$class, prior = c(0.1, 0.9)) #generate confusion matrix with prediction rates for each type
  
  #insert prediction rates into the matrix from step 4	
  coefs.rand[i] <- rand.confuse[1,1]*0.1 + rand.confuse[2,2]*0.9
}

expected <- data.frame(coefs.rand)

#calculate p-values
p <- sum(expected >= observed) / reps # positive effect
p
```

### MLM-FLF LDA
```{r}

flfmlm <- caps %>% filter(type %in% c("flf","mlm")) %>% filter(corrugations <= 40)

prop.table(table(flfmlm$type))
flfmlm.lda <- lda(type ~ mass + nare + culmen + tarsus + wingchord + wingdepth + tail, data = flfmlm, prior = c(0.302, 0.698), na.action = "na.omit")

flfmlm.morph = predict(flfmlm.lda)

flfmlmknifed_bird <- lda(type ~ mass + nare + culmen + tarsus + wingchord + wingdepth + tail, data = flfmlm, prior = c(0.302, 0.698), na.action = "na.omit", CV = TRUE)
summary(flfmlmknifed_bird)
flfmlmknifed_bird$posterior[1:5, ]
confuse <- confusion(flfmlm$type, flfmlmknifed_bird$class, prior = c(0.302, 0.698))

observed <- confuse[1,1]*0.302 + confuse[2,2]*0.698
observed
```

```{r}
#set number of permutations
reps <- 5000

#create initial random data
rand.data <- transform(flfmlm , type = sample(type))

#create empty matrix for randomized effect sizes
coefs.rand <- matrix(NA, nrow = reps, ncol = 1)

#fill matrix with results of tests with randomized data (this runs the DFA model 5000 times)
set.seed(1)
for(i in 1:reps){
  rand.data <- transform(rand.data, type = sample(type)) #randomize the sexes
  rand.fit <- lda(type ~ mass + nare + culmen + tarsus + wingchord + wingdepth + tail, data = rand.data, prior = c(0.302, 0.698), na.action = "na.omit") #run DFA with randomized type data
  
  rand.fit.p <- predict(rand.fit) #generate predictions based on DFA
  rand.confuse <- confusion(rand.data$type, rand.fit.p$class, prior = c(0.302, 0.698)) #generate confusion matrix with prediction rates for each type
  
  #insert prediction rates into the matrix from step 4	
  coefs.rand[i] <- rand.confuse[1,1]*0.302 + rand.confuse[2,2]*0.698
}

expected <- data.frame(coefs.rand)

#calculate p-values
p <- sum(expected >= observed) / reps # positive effect
p
```

### Genetic Sex LDA

```{r}

adults <- caps %>% filter(corrugations <= 40)

prop.table(table(adults$gensex))
gensex.lda <- lda(gensex ~ mass + nare + culmen + tarsus + wingchord + wingdepth + tail, data = adults, prior = c(0.353, 0.647), na.action = "na.omit")

sexknifed_bird <- lda(gensex ~ mass + nare + culmen + tarsus + wingchord + wingdepth + tail, data = adults, prior = c(0.353, 0.647), na.action = "na.omit", CV = TRUE)
summary(sexknifed_bird)
sexknifed_bird$posterior[1:5, ]
confuse <- confusion(adults$gensex, sexknifed_bird$class, prior = c(0.353, 0.647))

observed <- confuse[1,1]*0.353 + confuse[2,2]*0.647
observed
```

```{r}
#set number of permutations
reps <- 5000

#create initial random data
rand.data <- transform(adults, gensex = sample(gensex))

#create empty matrix for randomized effect sizes
coefs.rand <- matrix(NA, nrow = reps, ncol = 1)

#fill matrix with results of tests with randomized data (this runs the DFA model 5000 times)
set.seed(1)
for(i in 1:reps){
  rand.data <- transform(rand.data, gensex = sample(gensex)) #randomize the sexes
  rand.fit <- lda(gensex ~ mass + nare + culmen + tarsus + wingchord + wingdepth + tail, data = rand.data, prior = c(0.353, 0.647), na.action = "na.omit") #run DFA with randomized type data
  
  rand.fit.p <- predict(rand.fit) #generate predictions based on DFA
  rand.confuse <- confusion(rand.data$gensex, rand.fit.p$class, prior = c(0.353, 0.647)) #generate confusion matrix with prediction rates for each type
  
  #insert prediction rates into the matrix from step 4	
  coefs.rand[i] <- rand.confuse[1,1]*0.353 + rand.confuse[2,2]*0.647
}

expected <- data.frame(coefs.rand)

#calculate p-values
p <- sum(expected >= observed) / reps # positive effect
p
```