###The present script contains the R code used for the analyses of my thesis. 
###The name of the csv files may change according to how they are obtained from the Supplementary material (excel files).
###The supplementary material excel files include the main datasets for the 3 objectives for the 4 traits.
###For the analyses of small and large species, the corresponding datasets can be obtained by diving the main datasets based on the column "size".

getwd()
setwd("F:/lincoln uni/tesis/data analysis")

#########################################################################################
#########################################################################################
OBJ. 1 IDENTIFICATION OF INSULAR MORPHOLOGICAL PATTERNS
#########################################################################################
###########
# A) WINGS
###########

wing.island <- read.csv("wings.islands.csv",header=T)
wing.island
str(wing.island)

#########
## Comparing mainland and insular wing lengths at spp level 
#########

t.test(log(wing.island$mainland.mean),log(wing.island$insular.mean),paired=T)
#calculating Si, Sm and SR
wing.island$Si.sp <- wing.island$insular.mean
wing.island$Si.sp 
wing.island$Sm.sp <- wing.island$mainland.mean
wing.island$Sm.sp
wing.island$SR.sp <- wing.island$Si.sp/wing.island$Sm.sp
wing.island$SR.sp 
length(wing.island$mainland.sp[wing.island$SR.sp>1])
length(wing.island$mainland.sp[wing.island$SR.sp<1])

#regressing Si against Sm at spp level
summary(m1 <- lm(log(wing.island$Si.sp)~(log(wing.island$Sm.sp))+offset(log(wing.island$Sm.sp))))
plot(log(wing.island$Si.sp)~log(wing.island$Sm.sp),main="Regression of Si against Sm for wing length at species level",xlab="Mainland wing length(Log Sm)",ylab="Insular wing length(Log Si)")
abline(m1)
plot(((wing.island$SR.sp)*100)~log(wing.island$Sm.sp),main="Insular wing length as a function of mainland wing length at species level",xlab="Mainland wing length (Log Sm)",ylab="Insular wing length (as % Sm)")
abline(100,0,lty=2)

#########
## Comparing mainland and insular wing lengths at family level
#########

result_wing.island <- as.data.frame(matrix(nrow=length(levels(wing.island$family)),ncol=16))
names(result_wing.island) <- c("family","N","t.value","p.value","mn.of.diffs","df","lower.ci","upper.ci","Si","Sm","Sr","percentSR","intercept","slope","Pvalue.against.1","colour")

x11(9,7)
plot((wing.island$SR.sp*100)~log(wing.island$Sm.sp),type="n",main="Insular wing length as function of mainland wing length per Family",xlab="Mainland wing length (Log Sm)",ylab="Insular wing length (as % Sm)",cex.axis=1.3,cex.lab=1.3)
  abline(100,0,lty=2)

colours()
ccc <- c("green3","yellowgreen","darkolivegreen","purple","brown","brown","palevioletred","blue","green","orange","gold","khaki","darkolivegreen","darkgoldenrod","beige")

col.count <- 1

for (i in 1:length(levels(wing.island$family))) {
  dat <- wing.island[wing.island$family==levels(wing.island$family)[i],]
  N <- dim(dat)[1]
  family <- as.character(dat$family[1])
  result_wing.island$family[i] <- family
  result_wing.island$N[i] <- N
  
  result_wing.island$Si[i] <- (mean(dat$insular.mean))
  result_wing.island$Sm[i] <- (mean(dat$mainland.mean))
  result_wing.island$Sr[i] <- round((result_wing.island$Si[i]/result_wing.island$Sm[i]),3)
  result_wing.island$percentSR[i] <- round((result_wing.island$Si[i]/result_wing.island$Sm[i])*100,1)
  
  if(N<=10) {i = i+1}
  
  if(N>10) {
    # t-test
    tt <- t.test(log(dat$mainland.mean),log(dat$insular.mean),paired=T)
    result_wing.island$t.value[i] <- round(tt$statistic,2) # t-value
    result_wing.island$p.value[i] <- round(tt$p.value,3) # t-test P-value
    result_wing.island$mn.of.diffs[i] <- round(tt$estimate,2) # mean of the differences
    result_wing.island$df[i] <- tt$parameter # degrees of freedom
    result_wing.island$lower.ci[i] <- round(tt$conf.int[1],2) # lower bound of c.i.
    result_wing.island$upper.ci[i] <- round(tt$conf.int[2],2) # upper bound of c.i.
    
    # regression
    rr <- lm(log(dat$Si.sp)~log(dat$Sm.sp))
    or <- lm(log(dat$Si.sp)~log(dat$Sm.sp)+offset(log(dat$Sm.sp)))
    result_wing.island$intercept[i] <- round(rr$coef[1],2)
    result_wing.island$slope[i] <- round(rr$coef[2],2)
    result_wing.island$Pvalue.against.1[i] <- round(summary(or)$coefficients[2,4],3) # P-value is from offset regression testing if slope differs from one (not zero)
    
    # graph  
    yy <- (dat$SR.sp)*100
    xx <- log(dat$Sm.sp)
    points(yy~xx,pch=21,bg=ccc[col.count],cex=1.2)
    result_wing.island$colour[i] <- ccc[col.count]
    col.count <- col.count+1
  }
}

head(result_wing.island)
result_wing.island
write.table(result_wing.island,"wing.island_result.csv",sep=",",row.names=F)
leg <- result_wing.island[is.na(result_wing.island$t.value)==FALSE,c("family","colour")]
legend(locator(1),leg$family,fill=leg$colour,bty="n")


#########
## Comparing mainland and insular wing lengths at order level
#########

result_wing.island2 <- as.data.frame(matrix(nrow=length(levels(wing.island$order)),ncol=16))
names(result_wing.island2) <- c("order","N","t.value","p.value","mn.of.diffs","df","lower.ci","upper.ci","Si","Sm","Sr","percentSR","intercept","slope","Pvalue.against.1","colour")

x11(9,7)
plot((wing.island$SR.sp*100)~log(wing.island$Sm.sp),type="n",main="Insular wing length as function of mainland wing length per Order",xlab="Mainland wing length (Log Sm)",ylab="Insular wing length (as % Sm)",cex.axis=1.3,cex.lab=1.3)
  abline(100,0,lty=2)

colours()
ccc <- c("blue","green","purple","orange","brown","orange","purple","brown","grey","wheat","violet","khaki","darkolivegreen","darkgoldenrod","beige")

col.count <- 1

for (i in 1:length(levels(wing.island$order))) {
  dat <- wing.island[wing.island$order==levels(wing.island$order)[i],]
  N <- dim(dat)[1]
  order <- as.character(dat$order[1])
  result_wing.island2$order[i] <- order
  result_wing.island2$N[i] <- N
  
  result_wing.island2$Si[i] <- (mean(dat$insular.mean))
  result_wing.island2$Sm[i] <- (mean(dat$mainland.mean))
  result_wing.island2$Sr[i] <- round((result_wing.island2$Si[i]/result_wing.island2$Sm[i]),3)
  result_wing.island2$percentSR[i] <- round((result_wing.island2$Si[i]/result_wing.island2$Sm[i])*100,1)
  
  if(N<=10) {i = i+1}
  
  if(N>10) {
    # t-test
    tt <- t.test(log(dat$mainland.mean),log(dat$insular.mean),paired=T)
    result_wing.island2$t.value[i] <- round(tt$statistic,2) # t-value
    result_wing.island2$p.value[i] <- round(tt$p.value,3) # t-test P-value
    result_wing.island2$mn.of.diffs[i] <- round(tt$estimate,2) # mean of the differences
    result_wing.island2$df[i] <- tt$parameter # degrees of freedom
    result_wing.island2$lower.ci[i] <- round(tt$conf.int[1],2) # lower bound of c.i.
    result_wing.island2$upper.ci[i] <- round(tt$conf.int[2],2) # upper bound of c.i.
    
    # regression
    rr <- lm(log(dat$Si.sp)~log(dat$Sm.sp))
    or <- lm(log(dat$Si.sp)~log(dat$Sm.sp)+offset(log(dat$Sm.sp)))
    result_wing.island2$intercept[i] <- round(rr$coef[1],2)
    result_wing.island2$slope[i] <- round(rr$coef[2],2)
    result_wing.island2$Pvalue.against.1[i] <- round(summary(or)$coefficients[2,4],3)# P-value is from offset regression testing if slope differs from one (not zero)
    
    # graph  
    yy <- (dat$SR.sp)*100
    xx <- log(dat$Sm.sp)
    points(yy~xx,pch=21,bg=ccc[col.count],cex=1.2)
    result_wing.island2$colour[i] <- ccc[col.count]
    col.count <- col.count+1
  }
}

head(result_wing.island2)
result_wing.island2
write.table(result_wing.island2,"wing.island2_result.csv",sep=",",row.names=F)
leg <- result_wing.island2[is.na(result_wing.island2$t.value)==FALSE,c("order","colour")]
legend(locator(1),leg$order,fill=leg$colour,bty="n")

###########
# B) BILL
###########

bill.island <- read.csv("bills.islands.csv",header=T)
bill.island
str(bill.island)

#########
## Comparing mainland and insular bill lengths at spp level
#########

t.test(log(bill.island$mainland.mean),log(bill.island$insular.mean),paired=T)
#calculating Si, Sm and SR
bill.island$Si.sp <- (bill.island$insular.mean)
bill.island$Si.sp 
bill.island$Sm.sp <- (bill.island$mainland.mean)
bill.island$Sm.sp
bill.island$SR.sp <- bill.island$Si.sp/bill.island$Sm.sp
bill.island$SR.sp 
length(bill.island$mainland.sp[bill.island$SR.sp>1])
length(bill.island$mainland.sp[bill.island$SR.sp<1])

#regresing Si against Sm per sp
summary(m2 <- lm(log(bill.island$Si.sp)~(log(bill.island$Sm.sp))+offset(log(bill.island$Sm.sp))))
plot(log(bill.island$Si.sp)~log(bill.island$Sm.sp),main="Regression of Si against Sm for bill length at species level",xlab="Mainland bill length(Log Sm)",ylab="Insular bill length(Log Si)")
abline(m2)
plot(((bill.island$SR.sp)*100)~log(bill.island$Sm.sp),main="Insular bill length as a function of mainland bill length at species level",xlab="Mainland bill length (Log Sm)",ylab="Insular bill length (as % Sm)")
abline(100,0,lty=2)

#########
## Comparing mainland and insular bill lengths at family level
#########

result_bill.island <- as.data.frame(matrix(nrow=length(levels(bill.island$family)),ncol=16))
names(result_bill.island) <- c("family","N","t.value","p.value","mn.of.diffs","df","lower.ci","upper.ci","Si","Sm","Sr","percentSR","intercept","slope","Pvalue.against.1","colour")

x11(9,7)
plot((bill.island$SR.sp*100)~log(bill.island$Sm.sp),type="n",main="Insular bill length as function of mainland bill length per Family",xlab="Mainland bill length (Log Sm)",ylab="Insular bill length (as % Sm)",cex.axis=1.3,cex.lab=1.3)
  abline(100,0,lty=2)

colours()
ccc <- c("yellowgreen","purple","yellow1","brown","blue","darkolivegreen","brown","yellow","green","darkgoldenrod","beige","pink","palevioletred")

col.count <- 1

for (i in 1:length(levels(bill.island$family))) {
  dat <- bill.island[bill.island$family==levels(bill.island$family)[i],]
  N <- dim(dat)[1]
  family <- as.character(dat$family[1])
  result_bill.island$family[i] <- family
  result_bill.island$N[i] <- N
  
  result_bill.island$Si[i] <- (mean(dat$insular.mean))
  result_bill.island$Sm[i] <- (mean(dat$mainland.mean))
  result_bill.island$Sr[i] <- round((result_bill.island$Si[i]/result_bill.island$Sm[i]),3)
  result_bill.island$percentSR[i] <- round((result_bill.island$Si[i]/result_bill.island$Sm[i])*100,1)
  
  if(N<=10) {i = i+1}
  
  if(N>10) {
    # t-test
    tt <- t.test(log(dat$mainland.mean),log(dat$insular.mean),paired=T)
    result_bill.island$t.value[i] <- round(tt$statistic,2) # t-value
    result_bill.island$p.value[i] <- round(tt$p.value,3) # t-test P-value
    result_bill.island$mn.of.diffs[i] <- round(tt$estimate,2) # mean of the differences
    result_bill.island$df[i] <- tt$parameter # degrees of freedom
    result_bill.island$lower.ci[i] <- round(tt$conf.int[1],2) # lower bound of c.i.
    result_bill.island$upper.ci[i] <- round(tt$conf.int[2],2) # upper bound of c.i.
    
    # regression
    rr <- lm(log(dat$Si.sp)~log(dat$Sm.sp))
    or <- lm(log(dat$Si.sp)~log(dat$Sm.sp)+offset(log(dat$Sm.sp)))
    result_bill.island$intercept[i] <- round(rr$coef[1],2)
    result_bill.island$slope[i] <- round(rr$coef[2],2)
    result_bill.island$Pvalue.against.1[i] <- round(summary(or)$coefficients[2,4],3)# P-value is from offset regression testing if slope differs from one (not zero)
    
    # graph  
    yy <- (dat$SR.sp)*100
    xx <- log(dat$Sm.sp)
    points(yy~xx,pch=21,bg=ccc[col.count],cex=1.2)
    result_bill.island$colour[i] <- ccc[col.count]
    col.count <- col.count+1
  }
}

head(result_bill.island)
result_bill.island
write.table(result_bill.island,"bill.island_result.csv",sep=",",row.names=F)
leg <- result_bill.island[is.na(result_bill.island$t.value)==FALSE,c("family","colour")]
legend(locator(1),leg$family,fill=leg$colour,bty="n")

#########
## Comparing mainland and insular bill lengths at order level
#########

result_bill.island2 <- as.data.frame(matrix(nrow=length(levels(bill.island$order)),ncol=16))
names(result_bill.island2) <- c("order","N","t.value","p.value","mn.of.diffs","df","lower.ci","upper.ci","Si","Sm","Sr","percentSR","intercept","slope","Pvalue.against.1","colour")

x11(9,7)
plot((bill.island$SR.sp*100)~log(bill.island$Sm.sp),type="n",main="Insular bill length as function of mainland bill length per Order",xlab="Mainland bill length (Log Sm)",ylab="Insular bill length (as % Sm)",cex.axis=1.3,cex.lab=1.3)
  abline(100,0,lty=2)

colours()
ccc <- c("blue","yellow1","green","purple","orange","brown","purple","brown","grey","wheat","violet","khaki","darkolivegreen","darkgoldenrod","beige")

col.count <- 1

for (i in 1:length(levels(bill.island$order))) {
  dat <- bill.island[bill.island$order==levels(bill.island$order)[i],]
  N <- dim(dat)[1]
  order <- as.character(dat$order[1])
  result_bill.island2$order[i] <- order
  result_bill.island2$N[i] <- N
  
  result_bill.island2$Si[i] <- (mean(dat$insular.mean))
  result_bill.island2$Sm[i] <- (mean(dat$mainland.mean))
  result_bill.island2$Sr[i] <- round((result_bill.island2$Si[i]/result_bill.island2$Sm[i]),3)
  result_bill.island2$percentSR[i] <- round((result_bill.island2$Si[i]/result_bill.island2$Sm[i])*100,1)
  
  if(N<=10) {i = i+1}
  
  if(N>10) {
    # t-test
    tt <- t.test(log(dat$mainland.mean),log(dat$insular.mean),paired=T)
    result_bill.island2$t.value[i] <- round(tt$statistic,2) # t-value
    result_bill.island2$p.value[i] <- round(tt$p.value,3) # t-test P-value
    result_bill.island2$mn.of.diffs[i] <- round(tt$estimate,2) # mean of the differences
    result_bill.island2$df[i] <- tt$parameter # degrees of freedom
    result_bill.island2$lower.ci[i] <- round(tt$conf.int[1],2) # lower bound of c.i.
    result_bill.island2$upper.ci[i] <- round(tt$conf.int[2],2) # upper bound of c.i.
    
    # regression
    rr <- lm(log(dat$Si.sp)~log(dat$Sm.sp))
    or <- lm(log(dat$Si.sp)~log(dat$Sm.sp)+offset(log(dat$Sm.sp)))
    result_bill.island2$intercept[i] <- round(rr$coef[1],2)
    result_bill.island2$slope[i] <- round(rr$coef[2],2)
    result_bill.island2$Pvalue.against.1[i] <- round(summary(or)$coefficients[2,4],3)# P-value is from offset regression testing if slope differs from one (not zero)
    
    # graph  
    yy <- (dat$SR.sp)*100
    xx <- log(dat$Sm.sp)
    points(yy~xx,pch=21,bg=ccc[col.count],cex=1.2)
    result_bill.island2$colour[i] <- ccc[col.count]
    col.count <- col.count+1
  }
}

head(result_bill.island2)
result_bill.island2
write.table(result_bill.island2,"bill.island2_result.csv",sep=",",row.names=F)
leg <- result_bill.island2[is.na(result_bill.island2$t.value)==FALSE,c("order","colour")]
legend(locator(1),leg$order,fill=leg$colour,bty="n")

###########
# C) TARSUS
###########

tarsus.island <- read.csv("tarsuss.islands.csv",header=T)
tarsus.island
str(tarsus.island)

#########
## Comparing mainland and insular tarsus lengths at spp level
#########

t.test(log(tarsus.island$mainland.mean),log(tarsus.island$insular.mean),paired=T)
#calculating Si, Sm and SR
tarsus.island$Si.sp <- (tarsus.island$insular.mean)
tarsus.island$Si.sp 
tarsus.island$Sm.sp <- (tarsus.island$mainland.mean)
tarsus.island$Sm.sp
tarsus.island$SR.sp <- tarsus.island$Si.sp/tarsus.island$Sm.sp
tarsus.island$SR.sp 
length(tarsus.island$mainland.sp[tarsus.island$SR.sp>1])
length(tarsus.island$mainland.sp[tarsus.island$SR.sp<1])

#regresing Si against Sm per sp
summary(m3 <- lm(log(tarsus.island$Si.sp)~(log(tarsus.island$Sm.sp))+offset(log(tarsus.island$Sm.sp))))
plot(log(tarsus.island$Si.sp)~log(tarsus.island$Sm.sp),main="Regression of Si against Sm for tarsus length at species level",xlab="Mainland tarsus length(Log Sm)",ylab="Insular tarsus length(Log Si)")
abline(m3)
plot(((tarsus.island$SR.sp)*100)~log(tarsus.island$Sm.sp),main="Insular tarsus length as a function of mainland tarsus length at species level",xlab="Mainland tarsus length (Log Sm)",ylab="Insular tarsus length (as % Sm)")
abline(100,0,lty=2)

#########
## Comparing mainland and insular tarsus lengths at family level
#########

result_tarsus.island <- as.data.frame(matrix(nrow=length(levels(tarsus.island$family)),ncol=16))
names(result_tarsus.island) <- c("family","N","t.value","p.value","mn.of.diffs","df","lower.ci","upper.ci","Si","Sm","Sr","percentSR","intercept","slope","Pvalue.against.1","colour")

x11(9,7)
plot((tarsus.island$SR.sp*100)~log(tarsus.island$Sm.sp),type="n",main="Insular tarsus length as function of mainland tarsus length per Family",xlab="Mainland tarsus length (Log Sm)",ylab="Insular tarsus length (as % Sm)",cex.axis=1.3,cex.lab=1.3)
  abline(100,0,lty=2)

colours()
ccc <- c("green3","yellowgreen","darkolivegreen","lawngreen","purple","yellow1","purple","brown","grey","wheat","violet","khaki","darkolivegreen","darkgoldenrod","beige")

col.count <- 1

for (i in 1:length(levels(tarsus.island$family))) {
  dat <- tarsus.island[tarsus.island$family==levels(tarsus.island$family)[i],]
  N <- dim(dat)[1]
  family <- as.character(dat$family[1])
  result_tarsus.island$family[i] <- family
  result_tarsus.island$N[i] <- N
  
  result_tarsus.island$Si[i] <- (mean(dat$insular.mean))
  result_tarsus.island$Sm[i] <- (mean(dat$mainland.mean))
  result_tarsus.island$Sr[i] <- round((result_tarsus.island$Si[i]/result_tarsus.island$Sm[i]),3)
  result_tarsus.island$percentSR[i] <- round((result_tarsus.island$Si[i]/result_tarsus.island$Sm[i])*100,1)
  
  if(N<=10) {i = i+1}
  
  if(N>10) {
    # t-test
    tt <- t.test(log(dat$mainland.mean),log(dat$insular.mean),paired=T)
    result_tarsus.island$t.value[i] <- round(tt$statistic,2) # t-value
    result_tarsus.island$p.value[i] <- round(tt$p.value,3) # t-test P-value
    result_tarsus.island$mn.of.diffs[i] <- round(tt$estimate,2) # mean of the differences
    result_tarsus.island$df[i] <- tt$parameter # degrees of freedom
    result_tarsus.island$lower.ci[i] <- round(tt$conf.int[1],2) # lower bound of c.i.
    result_tarsus.island$upper.ci[i] <- round(tt$conf.int[2],2) # upper bound of c.i.
    
    # regression
    rr <- lm(log(dat$Si.sp)~log(dat$Sm.sp))
    or <- lm(log(dat$Si.sp)~log(dat$Sm.sp)+offset(log(dat$Sm.sp)))
    result_tarsus.island$intercept[i] <- round(rr$coef[1],2)
    result_tarsus.island$slope[i] <- round(rr$coef[2],2)
    result_tarsus.island$Pvalue.against.1[i] <- round(summary(or)$coefficients[2,4],3)# P-value is from offset regression testing if slope differs from one (not zero)
    
    # graph  
    yy <- (dat$SR.sp)*100
    xx <- log(dat$Sm.sp)
    points(yy~xx,pch=21,bg=ccc[col.count],cex=1.2)
    result_tarsus.island$colour[i] <- ccc[col.count]
    col.count <- col.count+1
  }
}

head(result_tarsus.island)
result_tarsus.island
write.table(result_tarsus.island,"tarsus.island_result.csv",sep=",",row.names=F)
leg <- result_tarsus.island[is.na(result_tarsus.island$t.value)==FALSE,c("family","colour")]
legend(locator(1),leg$family,fill=leg$colour,bty="n")

#########
## Comparing mainland and insular tarsus lengths at order level
#########

result_tarsus.island2 <- as.data.frame(matrix(nrow=length(levels(tarsus.island$order)),ncol=16))
names(result_tarsus.island2) <- c("order","N","t.value","p.value","mn.of.diffs","df","lower.ci","upper.ci","Si","Sm","Sr","percentSR","intercept","slope","Pvalue.against.1","colour")

x11(9,7)
plot((tarsus.island$SR.sp*100)~log(tarsus.island$Sm.sp),type="n",main="Insular tarsus length as function of mainland tarsus length per Order",xlab="Mainland tarsus length (Log Sm)",ylab="Insular tarsus length (as % Sm)",cex.axis=1.3,cex.lab=1.3)
  abline(100,0,lty=2)

colours()
ccc <- c("yellow1","green","purple","pink","green","orange","purple","brown","grey","wheat","violet","khaki","darkolivegreen","darkgoldenrod","beige")

col.count <- 1

for (i in 1:length(levels(tarsus.island$order))) {
  dat <- tarsus.island[tarsus.island$order==levels(tarsus.island$order)[i],]
  N <- dim(dat)[1]
  order <- as.character(dat$order[1])
  result_tarsus.island2$order[i] <- order
  result_tarsus.island2$N[i] <- N
  
  result_tarsus.island2$Si[i] <- (mean(dat$insular.mean))
  result_tarsus.island2$Sm[i] <- (mean(dat$mainland.mean))
  result_tarsus.island2$Sr[i] <- round((result_tarsus.island2$Si[i]/result_tarsus.island2$Sm[i]),3)
  result_tarsus.island2$percentSR[i] <- round((result_tarsus.island2$Si[i]/result_tarsus.island2$Sm[i])*100,1)
  
  if(N<=10) {i = i+1}
  
  if(N>10) {
    # t-test
    tt <- t.test(log(dat$mainland.mean),log(dat$insular.mean),paired=T)
    result_tarsus.island2$t.value[i] <- round(tt$statistic,2) # t-value
    result_tarsus.island2$p.value[i] <- round(tt$p.value,3) # t-test P-value
    result_tarsus.island2$mn.of.diffs[i] <- round(tt$estimate,2) # mean of the differences
    result_tarsus.island2$df[i] <- tt$parameter # degrees of freedom
    result_tarsus.island2$lower.ci[i] <- round(tt$conf.int[1],2) # lower bound of c.i.
    result_tarsus.island2$upper.ci[i] <- round(tt$conf.int[2],2) # upper bound of c.i.
    
    # regression
    rr <- lm(log(dat$Si.sp)~log(dat$Sm.sp))
    or <- lm(log(dat$Si.sp)~log(dat$Sm.sp)+offset(log(dat$Sm.sp)))
    result_tarsus.island2$intercept[i] <- round(rr$coef[1],2)
    result_tarsus.island2$slope[i] <- round(rr$coef[2],2)
    result_tarsus.island2$Pvalue.against.1[i] <- round(summary(or)$coefficients[2,4],3)# P-value is from offset regression testing if slope differs from one (not zero)
    
    # graph  
    yy <- (dat$SR.sp)*100
    xx <- log(dat$Sm.sp)
    points(yy~xx,pch=21,bg=ccc[col.count],cex=1.2)
    result_tarsus.island2$colour[i] <- ccc[col.count]
    col.count <- col.count+1
  }
}

head(result_tarsus.island2)
result_tarsus.island2
write.table(result_tarsus.island2,"tarsus.island2_result.csv",sep=",",row.names=F)
leg <- result_tarsus.island2[is.na(result_tarsus.island2$t.value)==FALSE,c("order","colour")]
legend(locator(1),leg$order,fill=leg$colour,bty="n")

###########
# D) WEIGHT
###########

weight.island <- read.csv("weights.islands.csv",header=T)
weight.island
str(weight.island)

#########
## Comparing mainland and insular weight measurements at spp level
#########

t.test(log(weight.island$mainland.mean),log(weight.island$insular.mean),paired=T)
#calculating Si, Sm  and SR
weight.island$Si.sp <- weight.island$insular.mean
weight.island$Si.sp 
weight.island$Sm.sp <- weight.island$mainland.mean
weight.island$Sm.sp
weight.island$SR.sp <- weight.island$Si.sp/weight.island$Sm.sp
weight.island$SR.sp 
length(weight.island$mainland.sp[weight.island$SR.sp>1])
length(weight.island$mainland.sp[weight.island$SR.sp<1])

#regresing Si against Sm per sp
summary(m4 <- lm(log(weight.island$Si.sp)~(log(weight.island$Sm.sp))+offset(log(weight.island$Sm.sp))))
plot(log(weight.island$Si.sp)~log(weight.island$Sm.sp),main="Regression of Si against Sm for weight at species level",xlab="Mainland weight(Log Sm)",ylab="Insular weight(Log Si)")
abline(m4)
plot(((weight.island$SR.sp)*100)~log(weight.island$Sm.sp),main="Insular weight as function of mainland weight at species level",xlab="Mainland weight (Log Sm)",ylab="Insular weight (as % Sm)")
abline(100,0,lty=2)

#########
## Comparing mainland and insular weight measurements at family level
#########

result_weight.island <- as.data.frame(matrix(nrow=length(levels(weight.island$family)),ncol=16))
names(result_weight.island) <- c("family","N","t.value","p.value","mn.of.diffs","df","lower.ci","upper.ci","Si","Sm","Sr","percentSR","intercept","slope","Pvalue.against.1","colour")

x11(9,7)
plot((weight.island$SR.sp*100)~log(weight.island$Sm.sp),type="n",main="Insular weight as function of mainland weight per Family",xlab="Mainland weight (Log Sm)",ylab="Insular weight (as % Sm)",cex.axis=1.3,cex.lab=1.3)
  abline(100,0,lty=2)

colours()
ccc <- c("yellowgreen","darkolivegreen","purple","brown","green","orange","purple","brown","grey","wheat","violet","khaki","darkolivegreen","darkgoldenrod","beige")

col.count <- 1

for (i in 1:length(levels(weight.island$family))) {
  dat <- weight.island[weight.island$family==levels(weight.island$family)[i],]
  N <- dim(dat)[1]
  family <- as.character(dat$family[1])
  result_weight.island$family[i] <- family
  result_weight.island$N[i] <- N
  
  result_weight.island$Si[i] <- (mean(dat$insular.mean))
  result_weight.island$Sm[i] <- (mean(dat$mainland.mean))
  result_weight.island$Sr[i] <- round((result_weight.island$Si[i]/result_weight.island$Sm[i]),3)
  result_weight.island$percentSR[i] <- round((result_weight.island$Si[i]/result_weight.island$Sm[i])*100,1)
  
  if(N<=10) {i = i+1}
  
  if(N>10) {
    # t-test
    tt <- t.test(log(dat$mainland.mean),log(dat$insular.mean),paired=T)
    result_weight.island$t.value[i] <- round(tt$statistic,2) # t-value
    result_weight.island$p.value[i] <- round(tt$p.value,3) # t-test P-value
    result_weight.island$mn.of.diffs[i] <- round(tt$estimate,2) # mean of the differences
    result_weight.island$df[i] <- tt$parameter # degrees of freedom
    result_weight.island$lower.ci[i] <- round(tt$conf.int[1],2) # lower bound of c.i.
    result_weight.island$upper.ci[i] <- round(tt$conf.int[2],2) # upper bound of c.i.
    
    # regression
    rr <- lm(log(dat$Si.sp)~log(dat$Sm.sp))
    or <- lm(log(dat$Si.sp)~log(dat$Sm.sp)+offset(log(dat$Sm.sp)))
    result_weight.island$intercept[i] <- round(rr$coef[1],2)
    result_weight.island$slope[i] <- round(rr$coef[2],2)
    result_weight.island$Pvalue.against.1[i] <- round(summary(or)$coefficients[2,4],3)# P-value is from offset regression testing if slope differs from one (not zero)
    
    # graph  
    yy <- (dat$SR.sp)*100
    xx <- log(dat$Sm.sp)
    points(yy~xx,pch=21,bg=ccc[col.count],cex=1.2)
    result_weight.island$colour[i] <- ccc[col.count]
    col.count <- col.count+1
  }
}

head(result_weight.island)
result_weight.island
write.table(result_weight.island,"weight.island_result.csv",sep=",",row.names=F)
leg <- result_weight.island[is.na(result_weight.island$t.value)==FALSE,c("family","colour")]
legend(locator(1),leg$family,fill=leg$colour,bty="n")


#########
## Comparing mainland and insular weight measurements at order level
#########

result_weight.island2 <- as.data.frame(matrix(nrow=length(levels(weight.island$order)),ncol=16))
names(result_weight.island2) <- c("order","N","t.value","p.value","mn.of.diffs","df","lower.ci","upper.ci","Si","Sm","Sr","percentSR","intercept","slope","Pvalue.against.1","colour")

x11(9,7)
plot((weight.island$SR.sp*100)~log(weight.island$Sm.sp),type="n",main="Insular weight as function of mainland weight per Order",xlab="Mainland weight (Log Sm)",ylab="Insular weight (as % Sm)",cex.axis=1.3,cex.lab=1.3)
  abline(100,0,lty=2)

colours()
ccc <- c("green","purple","orange","brown","green","orange","purple","brown","grey","wheat","violet","khaki","darkolivegreen","darkgoldenrod","beige")

col.count <- 1

for (i in 1:length(levels(weight.island$order))) {
  dat <- weight.island[weight.island$order==levels(weight.island$order)[i],]
  N <- dim(dat)[1]
  order <- as.character(dat$order[1])
  result_weight.island2$order[i] <- order
  result_weight.island2$N[i] <- N
  
  result_weight.island2$Si[i] <- (mean(dat$insular.mean))
  result_weight.island2$Sm[i] <- (mean(dat$mainland.mean))
  result_weight.island2$Sr[i] <- round((result_weight.island2$Si[i]/result_weight.island2$Sm[i]),3)
  result_weight.island2$percentSR[i] <- round((result_weight.island2$Si[i]/result_weight.island2$Sm[i])*100,1)
  
  if(N<=10) {i = i+1}
  
  if(N>10) {
    # t-test
    tt <- t.test(log(dat$mainland.mean),log(dat$insular.mean),paired=T)
    result_weight.island2$t.value[i] <- round(tt$statistic,2) # t-value
    result_weight.island2$p.value[i] <- round(tt$p.value,3) # t-test P-value
    result_weight.island2$mn.of.diffs[i] <- round(tt$estimate,2) # mean of the differences
    result_weight.island2$df[i] <- tt$parameter # degrees of freedom
    result_weight.island2$lower.ci[i] <- round(tt$conf.int[1],2) # lower bound of c.i.
    result_weight.island2$upper.ci[i] <- round(tt$conf.int[2],2) # upper bound of c.i.
    
    # regression
    rr <- lm(log(dat$Si.sp)~log(dat$Sm.sp))
    or <- lm(log(dat$Si.sp)~log(dat$Sm.sp)+offset(log(dat$Sm.sp)))
    result_weight.island2$intercept[i] <- round(rr$coef[1],2)
    result_weight.island2$slope[i] <- round(rr$coef[2],2)
    result_weight.island2$Pvalue.against.1[i] <- round(summary(or)$coefficients[2,4],3)# P-value is from offset regression testing if slope differs from one (not zero)
    
    # graph  
    yy <- (dat$SR.sp)*100
    xx <- log(dat$Sm.sp)
    points(yy~xx,pch=21,bg=ccc[col.count],cex=1.2)
    result_weight.island2$colour[i] <- ccc[col.count]
    col.count <- col.count+1
  }
}

head(result_weight.island2)
result_weight.island2
write.table(result_weight.island2,"weight.island2_result.csv",sep=",",row.names=F)
leg <- result_weight.island2[is.na(result_weight.island2$t.value)==FALSE,c("order","colour")]
legend(locator(1),leg$order,fill=leg$colour,bty="n")

############################################################################################################
# OBJ. 2 ASSESSING DIFFERENCES BETWEEN MAINLAND AND INSULAR SSD
############################################################################################################
##########
# A) WINGS
##########

wings.is.sd <- read.csv("wings.is.sd.csv",header=T)
wings.is.sd
str(wings.is.sd)

#Identifying M/F spp pairs
wings.is.sd$pair <- rep(c(1:(dim(wings.is.sd)[1]/2)),each=2)

#########
## Comparing mainland and insular SSD for wing lengths at spp level
#########

res_wings.is.sd <- as.data.frame(matrix(nrow=(dim(wings.is.sd)[1]/2),ncol=4))
names(res_wings.is.sd) <- c("mainland.sp","mainland.SSD","insular.sp","insular.SSD")

for(i in 1:max(wings.is.sd$pair)) {
  a1 <- wings.is.sd[wings.is.sd$pair==i,]
  res_wings.is.sd$mainland.sp[i] <- as.character(a1$mainland.sp[1])
  res_wings.is.sd$mainland.SSD[i] <- a1$mainland.mean[a1$sex=="M"]/a1$mainland.mean[a1$sex=="F"]
  res_wings.is.sd$insular.sp[i] <- as.character(a1$insular.sp[1])
  res_wings.is.sd$insular.SSD[i] <- a1$insular.mean[a1$sex=="M"]/a1$insular.mean[a1$sex=="F"] 
}

res_wings.is.sd
write.table(res_wings.is.sd,"wings.is.sd_res.csv",sep=",",row.names=F)

t.test(log(res_wings.is.sd$mainland.SSD),log(res_wings.is.sd$insular.SSD),alternative="less",paired=T)

#########
## Comparing mainland and insular SSD for wing lengths at family and order levels
#########

df <- data.frame(order="order",family="fam",N=0,pair=0,mssd=0,issd=0)
wings.is.sd <- wings.is.sd[order(wings.is.sd$family),] # reorder data frame by family

for (i in 1:length(levels(wings.is.sd$family))) {
  dat <- wings.is.sd[wings.is.sd$family==levels(wings.is.sd$family)[i],]
  df1 <- data.frame(order=unique(dat$order),family=unique(dat$family),N=1:length(unique(dat$pair)),pair=0,mssd=0,issd=0)
  for(j in 1:length(unique(dat$pair))) {
    dat$pair <- factor(dat$pair)
    pp <- dat[dat$pair==levels(dat$pair)[j],]
    family <- as.character(unique(dat$family))
    order <- as.character(unique(dat$order))
    N <- length(unique(dat$pair))
    pair <- as.numeric(as.character(unique(pp$pair)))
    mssd <- pp$mainland.mean[pp$sex=="M"]/pp$mainland.mean[pp$sex=="F"]
    issd <- pp$insular.mean[pp$sex=="M"]/pp$insular.mean[pp$sex=="F"]
    df1[j,3:6] <- c(N,pair,mssd,issd)
  }
  df <- rbind(df,df1)
}
  res.df <- df[-1,]
  row.names(res.df) <- 1:length(res.df$N)
  res.df$family <- factor(res.df$family) # remove unecessary level
  res.df$order <- factor(res.df$order) # remove unecessary level
  res.df
 
  res.df <- res.df[order(res.df$family),] # reorder dataframe by family before loop
  res.fam <- data.frame(family=unique(res.df$family),N=NA,t.value=NA,p.value=NA)
for(i in 1:length(res.fam$N))  {
  dat <- res.df[res.df$family==levels(res.df$family)[i],] 
  res.fam$N[i] <- dat$N[1] # record the number of pairs in the family
  if(res.fam$N[i]>5) {
    # t-test
    tt <- t.test(log(dat$mssd),log(dat$issd),alternative="less",paired=T)
    res.fam$t.value[i] <- round(tt$statistic,2) # t-value
    res.fam$p.value[i] <- round(tt$p.value,3) # t-test P-value
  }  
}
  
res.df <- res.df[order(res.df$order),] # reorder dataframe by order before loop
res.ord <- data.frame(order=unique(res.df$order),N=NA,t.value=NA,p.value=NA,no.pairs=NA)
for(i in 1:length(res.ord$N))  {
  dat <- res.df[res.df$order==levels(res.df$order)[i],]
  res.ord$N[i] <- length(unique(dat$family)) # count the number of families in the order
  res.ord$no.pairs[i] <- length(dat$N) # number of pairs in that order
  if(length(dat$N)>5){
    # t-test
    tt <- t.test(log(dat$mssd),log(dat$issd),alternative="less",paired=T)
    res.ord$t.value[i] <- round(tt$statistic,2) # t-value
    res.ord$p.value[i] <- round(tt$p.value,3) # t-test P-value
  }  
} 

res.fam
res.ord


##########
# B) BILL
##########

bills.is.sd <- read.csv("bills.is.sd.csv",header=T)
bills.is.sd
str(bills.is.sd)

#Identifying M/F spp pairs
bills.is.sd$pair <- rep(c(1:(dim(bills.is.sd)[1]/2)),each=2)

#########
## Comparing mainland and insular SSD for bill lengths at spp level
#########

res_bills.is.sd <- as.data.frame(matrix(nrow=(dim(bills.is.sd)[1]/2),ncol=4))
names(res_bills.is.sd) <- c("mainland.sp","mainland.SSD","insular.sp","insular.SSD")

for(i in 1:max(bills.is.sd$pair)) {
  a1 <- bills.is.sd[bills.is.sd$pair==i,]
  res_bills.is.sd$mainland.sp[i] <- as.character(a1$mainland.sp[1])
  res_bills.is.sd$mainland.SSD[i] <- a1$mainland.mean[a1$sex=="M"]/a1$mainland.mean[a1$sex=="F"]
  res_bills.is.sd$insular.sp[i] <- as.character(a1$insular.sp[1])
  res_bills.is.sd$insular.SSD[i] <- a1$insular.mean[a1$sex=="M"]/a1$insular.mean[a1$sex=="F"] 
}

res_bills.is.sd
write.table(res_bills.is.sd,"bills.is.sd_res.csv",sep=",",row.names=F)

t.test(log(res_bills.is.sd$mainland.SSD),log(res_bills.is.sd$insular.SSD),alternative="less",paired=T)

#########
## Comparing mainland and insular SSD for bill lengths at family and order level
#########

df <- data.frame(order="order",family="fam",N=0,pair=0,mssd=0,issd=0)
bills.is.sd <- bills.is.sd[order(bills.is.sd$family),] # reorder data frame by family

for (i in 1:length(levels(bills.is.sd$family))) {
  dat <- bills.is.sd[bills.is.sd$family==levels(bills.is.sd$family)[i],]
  df1 <- data.frame(order=unique(dat$order),family=unique(dat$family),N=1:length(unique(dat$pair)),pair=0,mssd=0,issd=0)
  for(j in 1:length(unique(dat$pair))) {
    dat$pair <- factor(dat$pair)
    pp <- dat[dat$pair==levels(dat$pair)[j],]
    family <- as.character(unique(dat$family))
    order <- as.character(unique(dat$order))
    N <- length(unique(dat$pair))
    pair <- as.numeric(as.character(unique(pp$pair)))
    mssd <- pp$mainland.mean[pp$sex=="M"]/pp$mainland.mean[pp$sex=="F"]
    issd <- pp$insular.mean[pp$sex=="M"]/pp$insular.mean[pp$sex=="F"]
    df1[j,3:6] <- c(N,pair,mssd,issd)
  }
  df <- rbind(df,df1)
}
  res.df <- df[-1,]
  row.names(res.df) <- 1:length(res.df$N)
  res.df$family <- factor(res.df$family) # remove unecessary level
  res.df$order <- factor(res.df$order) # remove unecessary level
  res.df
 
  res.df <- res.df[order(res.df$family),] # reorder dataframe by family before loop
  res.fam <- data.frame(family=unique(res.df$family),N=NA,t.value=NA,p.value=NA)
for(i in 1:length(res.fam$N))  {
  dat <- res.df[res.df$family==levels(res.df$family)[i],] 
  res.fam$N[i] <- dat$N[1] # record the number of pairs in the family
  if(res.fam$N[i]>5) {
    # t-test
    tt <- t.test(log(dat$mssd),log(dat$issd),alternative="less",paired=T)
    res.fam$t.value[i] <- round(tt$statistic,2) # t-value
    res.fam$p.value[i] <- round(tt$p.value,3) # t-test P-value
  }  
}
  
res.df <- res.df[order(res.df$order),] # reorder dataframe by order before loop
res.ord <- data.frame(order=unique(res.df$order),N=NA,t.value=NA,p.value=NA,no.pairs=NA)
for(i in 1:length(res.ord$N))  {
  dat <- res.df[res.df$order==levels(res.df$order)[i],]
  res.ord$N[i] <- length(unique(dat$family)) # count the number of families in the order
  res.ord$no.pairs[i] <- length(dat$N) # number of pairs in that order
  if(length(dat$N)>5){
    # t-test
    tt <- t.test(log(dat$mssd),log(dat$issd),alternative="less",paired=T)
    res.ord$t.value[i] <- round(tt$statistic,2) # t-value
    res.ord$p.value[i] <- round(tt$p.value,3) # t-test P-value
  }  
} 

res.fam
res.ord


##########
# C) TARSUS
##########

tarsuss.is.sd <- read.csv("tarsuss.is.sd.csv",header=T)
tarsuss.is.sd
str(tarsuss.is.sd)

#Identifying M/F spp pairs
tarsuss.is.sd$pair <- rep(c(1:(dim(tarsuss.is.sd)[1]/2)),each=2)

#########
## Comparing mainland and insular SSD for tarsus lengths at spp level
#########

res_tarsuss.is.sd <- as.data.frame(matrix(nrow=(dim(tarsuss.is.sd)[1]/2),ncol=4))
names(res_tarsuss.is.sd) <- c("mainland.sp","mainland.SSD","insular.sp","insular.SSD")

for(i in 1:max(tarsuss.is.sd$pair)) {
  a1 <- tarsuss.is.sd[tarsuss.is.sd$pair==i,]
  res_tarsuss.is.sd$mainland.sp[i] <- as.character(a1$mainland.sp[1])
  res_tarsuss.is.sd$mainland.SSD[i] <- a1$mainland.mean[a1$sex=="M"]/a1$mainland.mean[a1$sex=="F"]
  res_tarsuss.is.sd$insular.sp[i] <- as.character(a1$insular.sp[1])
  res_tarsuss.is.sd$insular.SSD[i] <- a1$insular.mean[a1$sex=="M"]/a1$insular.mean[a1$sex=="F"] 
}

res_tarsuss.is.sd
write.table(res_tarsuss.is.sd,"tarsuss.is.sd_res.csv",sep=",",row.names=F)

t.test(log(res_tarsuss.is.sd$mainland.SSD),log(res_tarsuss.is.sd$insular.SSD),alternative="less",paired=T)

#########
## Comparing mainland and insular SSD for tarsus lengths at family and order level
#########

df <- data.frame(order="order",family="fam",N=0,pair=0,mssd=0,issd=0)
tarsuss.is.sd <- tarsuss.is.sd[order(tarsuss.is.sd$family),] # reorder data frame by family

for (i in 1:length(levels(tarsuss.is.sd$family))) {
  dat <- tarsuss.is.sd[tarsuss.is.sd$family==levels(tarsuss.is.sd$family)[i],]
  df1 <- data.frame(order=unique(dat$order),family=unique(dat$family),N=1:length(unique(dat$pair)),pair=0,mssd=0,issd=0)
  for(j in 1:length(unique(dat$pair))) {
    dat$pair <- factor(dat$pair)
    pp <- dat[dat$pair==levels(dat$pair)[j],]
    family <- as.character(unique(dat$family))
    order <- as.character(unique(dat$order))
    N <- length(unique(dat$pair))
    pair <- as.numeric(as.character(unique(pp$pair)))
    mssd <- pp$mainland.mean[pp$sex=="M"]/pp$mainland.mean[pp$sex=="F"]
    issd <- pp$insular.mean[pp$sex=="M"]/pp$insular.mean[pp$sex=="F"]
    df1[j,3:6] <- c(N,pair,mssd,issd)
  }
  df <- rbind(df,df1)
}
  res.df <- df[-1,]
  row.names(res.df) <- 1:length(res.df$N)
  res.df$family <- factor(res.df$family) # remove unecessary level
  res.df$order <- factor(res.df$order) # remove unecessary level
  res.df
 
  res.df <- res.df[order(res.df$family),] # reorder dataframe by family before loop
  res.fam <- data.frame(family=unique(res.df$family),N=NA,t.value=NA,p.value=NA)
for(i in 1:length(res.fam$N))  {
  dat <- res.df[res.df$family==levels(res.df$family)[i],] 
  res.fam$N[i] <- dat$N[1] # record the number of pairs in the family
  if(res.fam$N[i]>5) {
    # t-test
    tt <- t.test(log(dat$mssd),log(dat$issd),alternative="less",paired=T)
    res.fam$t.value[i] <- round(tt$statistic,2) # t-value
    res.fam$p.value[i] <- round(tt$p.value,3) # t-test P-value
  }  
}
  
res.df <- res.df[order(res.df$order),] # reorder dataframe by order before loop
res.ord <- data.frame(order=unique(res.df$order),N=NA,t.value=NA,p.value=NA,no.pairs=NA)
for(i in 1:length(res.ord$N))  {
  dat <- res.df[res.df$order==levels(res.df$order)[i],]
  res.ord$N[i] <- length(unique(dat$family)) # count the number of families in the order
  res.ord$no.pairs[i] <- length(dat$N) # number of pairs in that order
  if(length(dat$N)>5){
    # t-test
    tt <- t.test(log(dat$mssd),log(dat$issd),alternative="less",paired=T)
    res.ord$t.value[i] <- round(tt$statistic,2) # t-value
    res.ord$p.value[i] <- round(tt$p.value,3) # t-test P-value
  }  
} 

res.fam
res.ord


##########
# D) WEIGHT
##########

weights.is.sd <- read.csv("weights.is.sd.csv",header=T)
weights.is.sd
str(weights.is.sd)

#Identifying M/F spp pairs
weights.is.sd$pair <- rep(c(1:(dim(weights.is.sd)[1]/2)),each=2)

#########
## Comparing mainland and insular SSD for weight measurements at spp level
#########

res_weights.is.sd <- as.data.frame(matrix(nrow=(dim(weights.is.sd)[1]/2),ncol=4))
names(res_weights.is.sd) <- c("mainland.sp","mainland.SSD","insular.sp","insular.SSD")

for(i in 1:max(weights.is.sd$pair)) {
  a1 <- weights.is.sd[weights.is.sd$pair==i,]
  res_weights.is.sd$mainland.sp[i] <- as.character(a1$mainland.sp[1])
  res_weights.is.sd$mainland.SSD[i] <- a1$mainland.mean[a1$sex=="M"]/a1$mainland.mean[a1$sex=="F"]
  res_weights.is.sd$insular.sp[i] <- as.character(a1$insular.sp[1])
  res_weights.is.sd$insular.SSD[i] <- a1$insular.mean[a1$sex=="M"]/a1$insular.mean[a1$sex=="F"] 
}

res_weights.is.sd
write.table(res_weights.is.sd,"weights.is.sd_res.csv",sep=",",row.names=F)

t.test(log(res_weights.is.sd$mainland.SSD),log(res_weights.is.sd$insular.SSD),alternative="less",paired=T)

#########
## Comparing mainland and insular SSD for weight measurements at family and order level
#########

df <- data.frame(order="order",family="fam",N=0,pair=0,mssd=0,issd=0)
weights.is.sd <- weights.is.sd[order(weights.is.sd$family),] # reorder data frame by family

for (i in 1:length(levels(weights.is.sd$family))) {
  dat <- weights.is.sd[weights.is.sd$family==levels(weights.is.sd$family)[i],]
  df1 <- data.frame(order=unique(dat$order),family=unique(dat$family),N=1:length(unique(dat$pair)),pair=0,mssd=0,issd=0)
  for(j in 1:length(unique(dat$pair))) {
    dat$pair <- factor(dat$pair)
    pp <- dat[dat$pair==levels(dat$pair)[j],]
    family <- as.character(unique(dat$family))
    order <- as.character(unique(dat$order))
    N <- length(unique(dat$pair))
    pair <- as.numeric(as.character(unique(pp$pair)))
    mssd <- pp$mainland.mean[pp$sex=="M"]/pp$mainland.mean[pp$sex=="F"]
    issd <- pp$insular.mean[pp$sex=="M"]/pp$insular.mean[pp$sex=="F"]
    df1[j,3:6] <- c(N,pair,mssd,issd)
  }
  df <- rbind(df,df1)
}
  res.df <- df[-1,]
  row.names(res.df) <- 1:length(res.df$N)
  res.df$family <- factor(res.df$family) # remove unecessary level
  res.df$order <- factor(res.df$order) # remove unecessary level
  res.df
 
  res.df <- res.df[order(res.df$family),] # reorder dataframe by family before loop
  res.fam <- data.frame(family=unique(res.df$family),N=NA,t.value=NA,p.value=NA)
for(i in 1:length(res.fam$N))  {
  dat <- res.df[res.df$family==levels(res.df$family)[i],] 
  res.fam$N[i] <- dat$N[1] # record the number of pairs in the family
  if(res.fam$N[i]>5) {
    # t-test
    tt <- t.test(log(dat$mssd),log(dat$issd),alternative="less",paired=T)
    res.fam$t.value[i] <- round(tt$statistic,2) # t-value
    res.fam$p.value[i] <- round(tt$p.value,3) # t-test P-value
  }  
}
  
res.df <- res.df[order(res.df$order),] # reorder dataframe by order before loop
res.ord <- data.frame(order=unique(res.df$order),N=NA,t.value=NA,p.value=NA,no.pairs=NA)
for(i in 1:length(res.ord$N))  {
  dat <- res.df[res.df$order==levels(res.df$order)[i],]
  res.ord$N[i] <- length(unique(dat$family)) # count the number of families in the order
  res.ord$no.pairs[i] <- length(dat$N) # number of pairs in that order
  if(length(dat$N)>5){
    # t-test
    tt <- t.test(log(dat$mssd),log(dat$issd),alternative="less",paired=T)
    res.ord$t.value[i] <- round(tt$statistic,2) # t-value
    res.ord$p.value[i] <- round(tt$p.value,3) # t-test P-value
  }  
} 

res.fam
res.ord


############################################################################################################
# OBJ. 3 META-ANALYSES (CALCULATING MEAN EFFECT SIZE, TESTING MODERATORS AND FITTING MODELS)
############################################################################################################

library(metafor)

library(vegan)

library(AICcmodavg)
search()

##########
# A) WINGS
##########

wingsmeta.islands <- read.csv("wingsmeta.islands.csv",header=T)
wingsmeta.islands
str(wingsmeta.islands)

#Standardizing data (different units)
wingsmeta.islands.stand <- decostand(wingsmeta.islands[,c("island.area","distance.mainland","sea.T")], method ="standardize")
wingsmeta.islands.stand
colnames(wingsmeta.islands.stand) <- c("island.area.stand","distance.mainland.stand","sea.T.stand")
str(wingsmeta.islands.stand)

#Mixing stand. data frame with original data
wingsmeta.isl <- cbind(wingsmeta.islands,wingsmeta.islands.stand)
str(wingsmeta.isl)

#Calculating outcome measurement
wingsmeta.islands2 <- escalc(n1i = insular.sample, n2i = mainland.sample, m1i = insular.mean, m2i = mainland.mean, sd1i = insular.std.dev, sd2i = mainland.std.dev, data = wingsmeta.isl, measure = "ROM", append = TRUE)
wingsmeta.islands2

#Accounting for lack of independence by including random effects argument
wingsmeta.islands2$species_pair <- paste(wingsmeta.islands2$mainland.sp, wingsmeta.islands2$insular.sp)
str(wingsmeta.islands2)

null.wings.1 <- rma.mv(yi, vi, random = ~ 1|family/species_pair, data = wingsmeta.islands2)
summary(null.wings.1)

area.wings.1 <- rma.mv(yi, vi, mods = cbind(island.area.stand), random = ~ 1|family/species_pair, data = wingsmeta.islands2)
summary(area.wings.1)

distance.wings.1 <- rma.mv(yi, vi, mods = cbind(distance.mainland.stand), random = ~ 1|family/species_pair, data = wingsmeta.islands2)
summary(distance.wings.1)

seaT.wings.1 <- rma.mv(yi, vi, mods = cbind(sea.T.stand), random = ~ 1|family/species_pair, data = wingsmeta.islands2)
summary(seaT.wings.1)

full.wings.1 <- rma.mv(yi, vi, mods = cbind(island.area.stand, distance.mainland.stand, sea.T.stand), random = ~ 1|family/species_pair, data = wingsmeta.islands2)
summary(full.wings.1)

##########
# B) BILLS
##########

billsmeta.islands <- read.csv("billsmeta.islands.csv",header=T)
billsmeta.islands
str(billsmeta.islands)

#Standardizing data (different units)
billsmeta.islands.stand <- decostand(billsmeta.islands[,c("island.area","distance.mainland","sea.T")], method ="standardize")
billsmeta.islands.stand
colnames(billsmeta.islands.stand) <- c("island.area.stand","distance.mainland.stand","sea.T.stand")
str(billsmeta.islands.stand)

#Mixing stand. data frame with original data
billsmeta.isl <- cbind(billsmeta.islands,billsmeta.islands.stand)
str(billsmeta.isl)

#Calculating outcome measurement
billsmeta.islands2 <- escalc(n1i = insular.sample, n2i = mainland.sample, m1i = insular.mean, m2i = mainland.mean, sd1i = insular.std.dev, sd2i = mainland.std.dev, data = billsmeta.isl, measure = "ROM", append = TRUE)
billsmeta.islands2

#Accounting for lack of independence by including random effects argument
billsmeta.islands2$species_pair <- paste(billsmeta.islands2$mainland.sp, billsmeta.islands2$insular.sp)
billsmeta.islands2

null.bills.1 <- rma.mv(yi, vi, random = ~ 1|family/species_pair, data = billsmeta.islands2)
summary(null.bills.1)

area.bills.1 <- rma.mv(yi, vi, mods = cbind(island.area.stand), random = ~ 1|family/species_pair, data = billsmeta.islands2)
summary(area.bills.1)

distance.bills.1 <- rma.mv(yi, vi, mods = cbind(distance.mainland.stand), random = ~ 1|family/species_pair, data = billsmeta.islands2)
summary(distance.bills.1)

seaT.bills.1 <- rma.mv(yi, vi, mods = cbind(sea.T.stand), random = ~ 1|family/species_pair, data = billsmeta.islands2)
summary(seaT.bills.1)

full.bills.1 <- rma.mv(yi, vi, mods = cbind(island.area.stand, distance.mainland.stand, sea.T.stand), random = ~ 1|family/species_pair, data = billsmeta.islands2)
summary(full.bills.1)

##########
# C) TARSUS
##########

tarsusmeta.islands <- read.csv("tarsussmeta.islands.csv",header=T)
tarsusmeta.islands
str(tarsusmeta.islands)

#Standardizing data (different units)
tarsusmeta.islands.stand <- decostand(tarsusmeta.islands[,c("island.area","distance.mainland","sea.T")], method ="standardize")
tarsusmeta.islands.stand
colnames(tarsusmeta.islands.stand) <- c("island.area.stand","distance.mainland.stand","sea.T.stand")
str(tarsusmeta.islands.stand)

#Mixing stand. data frame with original data
tarsusmeta.isl <- cbind(tarsusmeta.islands,tarsusmeta.islands.stand)
str(tarsusmeta.isl)

#Calculating outcome measurement
tarsusmeta.islands2 <- escalc(n1i = insular.sample, n2i = mainland.sample, m1i = insular.mean, m2i = mainland.mean, sd1i = insular.std.dev, sd2i = mainland.std.dev, data = tarsusmeta.isl, measure = "ROM", append = TRUE)
tarsusmeta.islands2

#Accounting for lack of independence by including random effects argument
tarsusmeta.islands2$species_pair <- paste(tarsusmeta.islands2$mainland.sp, tarsusmeta.islands2$insular.sp)
tarsusmeta.islands2

null.tarsus.1 <- rma.mv(yi, vi, random = ~ 1|family/species_pair, data = tarsusmeta.islands2)
summary(null.tarsus.1)

area.tarsus.1 <- rma.mv(yi, vi, mods = cbind(island.area.stand), random = ~ 1|family/species_pair, data = tarsusmeta.islands2)
summary(area.tarsus.1)

distance.tarsus.1 <- rma.mv(yi, vi, mods = cbind(distance.mainland.stand), random = ~ 1|family/species_pair, data = tarsusmeta.islands2)
summary(distance.tarsus.1)

seaT.tarsus.1 <- rma.mv(yi, vi, mods = cbind(sea.T.stand), random = ~ 1|family/species_pair, data = tarsusmeta.islands2)
summary(seaT.tarsus.1)

full.tarsus.1 <- rma.mv(yi, vi, mods = cbind(island.area.stand, distance.mainland.stand, sea.T.stand), random = ~ 1|family/species_pair, data = tarsusmeta.islands2)
summary(full.tarsus.1)

##########
# D) WEIGHT
##########

weightsmeta.islands <- read.csv("weightsmeta.islands.csv",header=T)
weightsmeta.islands
str(weightsmeta.islands)

#Standardizing data (different units)
weightsmeta.islands.stand <- decostand(weightsmeta.islands[,c("island.area","distance.mainland","sea.T")], method ="standardize")
weightsmeta.islands.stand
colnames(weightsmeta.islands.stand) <- c("island.area.stand","distance.mainland.stand","sea.T.stand")
str(weightsmeta.islands.stand)

#Mixing stand. data frame with original data
weightsmeta.isl <- cbind(weightsmeta.islands,weightsmeta.islands.stand)
str(weightsmeta.isl)

#Calculating outcome measurement
weightsmeta.islands2 <- escalc(n1i = insular.sample, n2i = mainland.sample, m1i = insular.mean, m2i = mainland.mean, sd1i = insular.std.dev, sd2i = mainland.std.dev, data = weightsmeta.isl, measure = "ROM", append = TRUE)
weightsmeta.islands2

#Meta-analytic models without random effects argument (for weight data of small species use the 'rma.mv' code like for the other traits)
null.weight.1 <- rma(yi, vi, data = weightsmeta.islands2)
summary(null.weight.1)

area.weight.1 <- rma(yi, vi, mods = cbind(island.area.stand), data = weightsmeta.islands2)
summary(area.weight.1)

distance.weight.1 <- rma(yi, vi, mods = cbind(distance.mainland.stand), data = weightsmeta.islands2)
summary(distance.weight.1)

seaT.weight.1 <- rma(yi, vi, mods = cbind(sea.T.stand), data = weightsmeta.islands2)
summary(seaT.weight.1)

full.weight.1 <- rma(yi, vi, mods = cbind(island.area.stand, distance.mainland.stand, sea.T.stand), data = weightsmeta.islands2)
summary(full.weight.1)