library(readxl)
hh1 <- read_excel("Skylark flight data (HH site 2)(CSV).xlsx", 
    sheet = "Sheet2")
hh2 <- read_excel("SKylark flight data (HH site 3).xlsx", 
    sheet = "Sheet2")
md <- read_excel("Skylark flight data (MD site 1)(CSV).xlsx", 
    sheet = "Sheet2")
hh1$site<-"hh1"
hh2$site<-"hh2"
md$site<-"md"
names(md)<-names(hh1)
d<-rbind(hh1,md)
d<-rbind(hh2,d)
d %>% ggplot(aes(x=site,y=Song_Length)) -> g0
g0+geom_boxplot()

aqm::ci(g0)

mod<-lm(data=d, Song_Length~site)
anova(mod)
## Analysis of Variance Table
## 
## Response: Song_Length
##           Df  Sum Sq Mean Sq F value Pr(>F)
## site       2   64503   32252  1.5387 0.2205
## Residuals 86 1802526   20960
summary(mod)
## 
## Call:
## lm(formula = Song_Length ~ site, data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -254.36  -72.26  -13.26   68.64  746.64 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   158.43      38.69   4.095 9.52e-05 ***
## sitehh2        45.83      42.90   1.068   0.2884    
## sitemd         95.93      54.72   1.753   0.0831 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 144.8 on 86 degrees of freedom
## Multiple R-squared:  0.03455,    Adjusted R-squared:  0.0121 
## F-statistic: 1.539 on 2 and 86 DF,  p-value: 0.2205
performance::check_model(mod)