#############################################
### Feather analaysis preliminary stats
### Cheek R and Alza L
### Created 25-Jan-2017. Last update 21-Oct-2017
###Analysis of feather data
#############################################

                
## load libraries
library(lme4)
library(lmerTest)
library(pwr)
library(simpleboot)
library(car)
library(MASS)


#############################################
##Power analysis

#power of t-test down length
#?pwr.t.test
#d=Effect Size Calculator for T-Test. For the independent samples T-test, Cohen's d is determined by calculating the mean difference between your two groups, and then dividing the result by the pooled standard deviation.
pwr.t.test (n=6, sig.level=0.05, d=1.605, type= c("two.sample"))
#power=0.71

#power of ANOVA down length
#?pwr.anova.test
#f=Effect Size needs the values from the ANOVA test
pwr.anova.test(k= 2, n=6, f= 1.211986, sig.level= .05)
#power = 0.97

#power of t-test down barbs
pwr.t.test(n=6, sig.level=0.05, d=1.289, type=c("two.sample"))
#power = 0.52

#power of ANOVA down barbs. calculate f using population
pwr.anova.test(k= 2, n=6, f= 3.4536, sig.level= .05)
#power = 1
#Avoid t-test and will use a muli variable model instead because its easier to build


#############################################
## set working directory
setwd("C:/Users/lmcheek/Desktop/M.armata data")


#############################################
#Down Feathers Analysis
## load data
down<-read.csv("M. armata_Down_LM1.csv")
names(down)

##Two-way ANOVA with repeated measures down length 
down.mean<- aggregate(down$Length, by=list(down$CatalogueNumber, down$Location), FUN='mean')
colnames(down.mean) <- c("CatalogueNumber", "Location", "Length")
down.mean <- down.mean[order(down.mean$Location),]
head(down.mean)
downlength.aov <- with(down.mean, aov(Length~Location))
summary(downlength.aov)
?aov
                                        
#normality test Barbs
Barbs <- down$Barbs
shapiro.test(Barbs)

#normality test Length
Length <- down$Length
shapiro.test(Length)
#Normality assumption met. Move onto tests

#######
##mixed model to control for repeated measures
## Down Length model
?lmer
class(down$CatalogueNumber)
downlength <- lmer(Length~Location1+(1|CatalogueNumber),na.action=na.omit, data=down)
plot(downlength)
summary(downlength)
shapiro.test(resid(downlength))
difflsmeans(downlength, test.effs = "Location1")#for linear mixed-effects models
boxplot(Length ~ Location1, data=down, main= "Down Feather Length",ylab="Length (mm)", ylim= c(13,21), names= c("Highland", "Lowland"), col= (c("lightgrey", "darkgrey")), boxwex= .5)
anova(downlength)
?confint
confint(downlength, level=0.95, oldNames = FALSE)

#Down Barb number model Poisson distribution because count data
?glmer
downbarb <- glmer(Barbs~Location1+(1|CatalogueNumber), data=down, family= "poisson")
plot(downbarb)
summary(downbarb)
shapiro.test(resid(downbarb))
boxplot(Barbs~Location1, data=down, main= "Down Barb Number",ylab="Number of Barbs",ylim=c(26, 47), names= c("Highland", "Lowland"), col= (c("light grey", "dark grey")), boxwex= 0.5)
anova(downbarb)
confint(downbarb, level=0.95, oldNames = FALSE)#for count data


#############################################
#Contour Feather Analysis 
## load data
contour<-read.csv("M. armata_Contour_LM1.csv")
names(contour)

##Two-way ANOVA with repeated measures down length 

contourlength.aov <- aov(Length~Location1+CatalogueNumber, data=contour)
summary(downlength.aov)

#Prop plum model
propplummod<-lmer(nplumulaceous~Location1+(1|CatalogueNumber), na.action=na.omit, data=contour)
plot(propplummod)
shapiro.test(resid(propplummod))
summary(propplummod)
difflsmeans(propplummod, test.effs="Location1")
boxplot(nplumulaceous~Location, data=contour, main= "Proportion Plumulaceous",ylab="Proportion Plumulaceous", names= c("Highland", "Lowland"), col= (c("light grey", "dark grey")), boxwex= 0.5)
anova(propplummod)
confint(propplummod, oldNames= FALSE)

#Contour length model 
Contour.length<- lmer(Length~Location1+(1|CatalogueNumber), na.action=na.omit, data=contour)
plot(Contour.length)
shapiro.test(resid(Contour.length))
summary(Contour.length)
difflsmeans(Contour.length, test.effs="Location1")#for linear models
boxplot(Length~Location, data=contour, main= "Contour Feather Length", ylab="Length (mm)", names= c("Highland", "Lowland"), col= (c("light grey", "dark grey")), boxwex= 0.5)
anova(Contour.length)
confint(Contour.length, oldNames=FALSE)

#Barb Plum model 
Barbs.Plum<-glmer(BarbsPlum~Location1+(1|CatalogueNumber), na.action = na.omit, data=contour, family="poisson")
plot(Barbs.Plum)
shapiro.test(resid(Barbs.Plum))
summary(Barbs.Plum)
boxplot(BarbsPlum~Location, data=contour, main= "Barb Density Plumulaceous", ylab="Barb Number Plumulaceous", names= c("Highland", "Lowland"), col= (c("light grey", "dark grey")), boxwex= 0.5)
anova(Barbs.Plum)
confint(Barbs.Plum, oldNames= FALSE)

#Barb Penn model
Barbs.Penn <- glmer(BarbsPenn~Location1+(1|CatalogueNumber), na.action=na.omit, data=contour, family="poisson")
plot(Barbs.Penn)
shapiro.test(resid(Barbs.Penn))
hist(resid(Barbs.Penn))
boxplot(resid(Barbs.Penn))

bar<-resid(Barbs.Penn)
qqnorm(bar)
qqline(bar)

summary(Barbs.Penn)
boxplot(BarbsPenn~Location, data=contour, main= "Barb Density Pennaceous", ylab="Barb Number Pennaceous", names= c("Highland", "Lowland"), col= (c("light grey", "dark grey")), boxwex= 0.5)
anova(Barbs.Penn)
confint(Barbs.Penn, oldNames=FALSE)


#############################################
#Barb Penn model No Outlier
outlier.contour <- read.csv("M. armata_Contour_LM_No Outlier.csv")
Barb.Penn.out <- glmer(BarbsPennaceous~Location1+(1|CatalogueNumber), na.action=na.omit, data= outlier.contour, family="poisson")
shapiro.test(Barb.Penn.out)
plot(Barb.Penn.out)

qqnorm(resid(Barb.Penn.out))
qqline(resid(Barb.Penn.out))

shapiro.test(resid(Barb.Penn.out))
hist(resid(Barb.Penn.out))
difflsmeans(Barb.Penn.out, test.effs="Location")
boxplot(resid(Barb.Penn.out))
anova(Barb.Penn.out)
confint(Barb.Penn.out)