setwd('C:/R') G <- read.csv('Garanardcam.csv') SGcam02 <- read.csv('SGcam02.csv') Balmoral <-read.csv('Balmoralcam.csv') Ballywhite <- read.csv('Ballywhitecam.csv') Ballywalter<- read.csv('Ballywaltercam.csv') MS2 <- read.csv('MS02cam.csv') MS3 <- read.csv('MS03cam.csv') MS4<- read.csv('MS4cam.csv') TM2cam<- read.csv('TM2cam.csv') Danesford<- read.csv('Danesfordcam.csv') Stranmillis <- read.csv('Stranmilliscam.csv') MS4<- read.csv('MS4cam.csv') NW1<- read.csv('Nugents01cam.csv') NW2 <- read.csv('Nugents02cam.csv') squirrels<- read.csv('2daysquirrelmaster.csv') Lacefield <-read.csv('Lacefieldcam.csv') UFTM <- read.csv('UFTMcam.csv') NF <- read.csv('Newforgecam.csv') T <- read.csv('Torroschcam.csv') Greys <- read.csv('GreysNoVisits3days.csv') Reds <- read.csv('RedsNoVisits3day.csv') SG03 <- read.csv('SGcam3.csv') TM01 <- read.csv('TM01cam.csv') TM02<- read.csv('TM02cam.csv') DM <- read.csv('Dunmurrycam.csv') BP <- read.csv('Bristowcam.csv') behaviour<- read.csv('behaviour2dayw0.csv') visits<- read.csv("SquirrelsNoVisits2daysPrePost.csv") #master dataset write.csv(a3, file = "squirrelsmaster.csv") squirrels$TIME.AT.FEEDER<- as.numeric(squirrels$TIME.AT.FEEDER) a1 <- group_by(squirrels, SPECIES, VISIT.NUMBER, TREATMENT, SITE) a2 <- select(a1, SITE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) a3 <- summarise(a2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) VISITSc<-summarySE(visits, measurevar="VISITS", groupvars=c("SPECIES","TREATMENT")) DURc<-summarySE(a3, measurevar="TIME.AT.FEEDER", groupvars=c("SPECIES","TREATMENT")) lattice <- grid.arrange(arrangeGrob(Visbar, durbar, vbar, fbar, ncol = 2, nrow = 2), mylegend, nrow = 2, heights=c(10, 1)) png("C:/R/latticefinal02.png", width = 18, height = 21, units = 'cm', res = 300) ggarrange(p1, p2, p3, p4, ncol=2, nrow=2, common.legend = TRUE, legend="bottom") durbar <-ggplot(DURc, aes(SPECIES, TIME.AT.FEEDER, fill = TREATMENT, width = 0.5)) + xlab("Species") + ylab("Duration of visits to feeder (s)")+ geom_bar(stat="identity", position=position_dodge(width=.8), colour = 'black') + geom_errorbar(aes(ymin=TIME.AT.FEEDER, ymax=TIME.AT.FEEDER+se), width=.2, # Width of the error bars position=position_dodge(.8)) + scale_fill_manual(values=c("grey", "black"))+ theme_bw() + scale_y_continuous(expand=c(0,0), limits=c(0, 80), breaks=seq(0, 80, by=10))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text.x=element_text(colour = "black", face = "bold"))+ guides(fill=FALSE)+labs(title = "B") mylegend<-g_legend(durbar) mylegend,nrow=2,heights=c(10, 1) Visbar<- ggplot(VISITSc, aes(SPECIES, VISITS, fill = TREATMENT, width = 0.5)) + xlab("Species") + ylab("Number of vists / day")+ geom_bar(stat="identity", position=position_dodge(width=.8), colour = 'black') + geom_errorbar(aes(ymin=VISITS, ymax=VISITS+se), width=.2, # Width of the error bars position=position_dodge(.8)) + scale_fill_manual(values=c("grey", "black"))+ theme_bw() + scale_y_continuous(expand=c(0,0), limits=c(0, 110), breaks=seq(0, 110, by=10))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text.x=element_text(colour = "black", face = "bold")) + guides(fill=FALSE)+ labs(title = "A") a4 <- subset(a3, SPECIES=="Sciurus vulgaris", select=TREATMENT:TIME.AT.FEEDER) a5 <- subset(a3, SPECIES == 'Sciurus carolinensis', select=TREATMENT:TIME.AT.FEEDER) a7 <- subset(squirrels, SPECIES == 'Sciurus carolinensis', select=TIME.AT.FEEDER) a3 <- summarise(a, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) a9 <- subset(squirrels, SPECIES == 'Sciurus vulgaris', select=TIME.AT.FEEDER) a4c <-summarySE(a4, measurevar="TIME.AT.FEEDER", groupvars=c("SITE","TREATMENT")) a5c <-summarySE(a5, measurevar="TIME.AT.FEEDER", groupvars=c("SITE","TREATMENT")) BFc <-summarySE(behaviour, measurevar="FEEDING", groupvars=c("SPECIES","TREATMENT")) BVc <-summarySE(behaviour, measurevar="VIGILENCE", groupvars=c("SPECIES","TREATMENT")) BVc$TREATMENT<- factor(BVc$TREATMENT, levels=c("Pre","Post")) BFc$TREATMENT<- factor(BFc$TREATMENT, levels=c("Pre","Post")) VISITSc$TREATMENT<- factor(VISITSc$TREATMENT, levels=c("Pre","Post")) DURc$TREATMENT<- factor(DURc$TREATMENT, levels=c("Pre","Post")) sedatcSE$Season <- factor(sedatcSE$Season, levels=c("Winter", "Spring", "Summer", "Autumn")) fbar <-ggplot(BFc, aes(SPECIES, FEEDING, fill = TREATMENT, width = 0.5)) + xlab("Species") + ylab("Percentage feeding (%)")+ geom_bar(stat="identity", position=position_dodge(width=.8), colour = 'black') + geom_errorbar(aes(ymin=FEEDING, ymax=FEEDING+se), width=.2, # Width of the error bars position=position_dodge(.8)) + scale_fill_manual(values=c("grey", "black"))+ theme_bw() + scale_y_continuous(expand=c(0,0), limits=c(0, 100), breaks=seq(0, 100, by=10))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text.x=element_text(colour = "black", face = "bold")) + guides(fill=FALSE)+ labs(title = "D") vbar <-ggplot(BVc, aes(SPECIES, VIGILENCE, fill = TREATMENT, width = 0.5)) + xlab("Species") + ylab("Percentage time vigilant (%)")+ geom_bar(stat="identity", position=position_dodge(width=.8), colour = 'black') + geom_errorbar(aes(ymin=VIGILENCE, ymax=VIGILENCE+se), width=.2, # Width of the error bars position=position_dodge(.8)) + scale_fill_manual(values=c("grey", "black"))+ theme_bw() + scale_y_continuous(expand=c(0,0), limits=c(0, 25), breaks=seq(0, 25, by=5))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text.x=element_text(colour = "black", face = "bold")) + guides(fill=FALSE)+ labs(title = "C") behaviour$VIGILENCE<- as.numeric(behaviour$VIGILENCE) mod1<- lmer(VIGILENCE~TREATMENT+SPECIES+TREATMENT:SPECIES + (1|SITE), data = behaviour) summary(mod1) anova(mod1) hist(behaviour$VIGILENCE) vrs <- subset(behaviour, SPECIES == 'Red Squirrel (Sciurus vulgaris)', select=c("VIGILENCE", "TREATMENT", "SITE")) grs <- subset(behaviour, SPECIES == 'Grey Squirrel (Sciurus carolinensis)', select=c("VIGILENCE", "TREATMENT", "SITE")) frs <- subset(behaviour, SPECIES == 'Red Squirrel (Sciurus vulgaris)', select=c("FEEDING", "TREATMENT", "SITE")) fgs <- subset(behaviour, SPECIES == 'Grey Squirrel (Sciurus carolinensis)', select=c("FEEDING", "TREATMENT", "SITE")) t1<- table(vrs$VIGILENCE, vrs$TREATMENT) chisq.test(t1) #non.signif t.test(VIGILENCE~TREATMENT, data = vrs) library(nlme) m1 <- lme(VIGILENCE~TREATMENT,random=~1|SITE, method = 'REML', data=vrs) anova(m1) m1 <- lme(VIGILENCE~TREATMENT,random=~1|SITE, method = 'REML', data=grs) anova(m1) m1 <- lme(FEEDING~TREATMENT,random=~1|SITE, method = 'REML', data=frs) anova(m1) m1 <- lme(FEEDING~TREATMENT,random=~1|SITE, method = 'REML', data=fgs) anova(m1) m1 <- lme(FEEDING~TREATMENT,random=~1|SITE, method = 'REML', data=fgs) anova(m1) a3 <- summarise(a, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) a9 <- subset(squirrels, SPECIES == 'Sciurus vulgaris', select=TIME.AT.FEEDER) mod1<- glm(FEEDING~TREATMENT+SPECIES+TREATMENT:SPECIES, data = behaviour) summary(mod1) redsbar <-ggplot(a4c, aes(SITE, TIME.AT.FEEDER, fill = TREATMENT)) + geom_bar(stat="identity", position = "dodge") + geom_errorbar(aes(ymin=TIME.AT.FEEDER-a4c$ci, ymax=TIME.AT.FEEDER+a4c$ci), width=.2, # Width of the error bars position=position_dodge(.9)) + scale_fill_brewer(palette = "Set1") + scale_y_continuous(limits = c(0, 200)) + ggtitle('Red Squirrels (Sciurus vulgaris)')+ xlab("Site") + ylab("Duration of visits (Secs)") squirrelbar <- ggplot(a3, aes(SPECIES, TIME.AT.FEEDER, fill = TREATMENT)) + geom_bar(stat="l", position = "dodge") mod1<-lm(TIME.AT.FEEDER~SPECIES+TREATMENT, data = a3) summary(mod1) geom_errorbar(aes(ymin=TIME.AT.FEEDER-a4c$ci, ymax=TIME.AT.FEEDER+a4c$ci), width=.2, # Width of the error bars position=position_dodge(.9)) a3$TIME.AT.FEEDER <- log(a3$TIME.AT.FEEDER) a3$TIME.AT.FEEDER<- as.numeric(a3$TIME.AT.FEEDER) squirrelbox <- ggplot(a3, aes(y = log(TIME.AT.FEEDER), x = SPECIES, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() greysbar <-ggplot(a5c, aes(SITE, TIME.AT.FEEDER, fill = TREATMENT)) + geom_bar(stat="identity", position = "dodge") + geom_errorbar(aes(ymin=TIME.AT.FEEDER-a5c$ci, ymax=TIME.AT.FEEDER+a5c$ci), width=.2, # Width of the error bars position=position_dodge(.9)) + scale_fill_brewer(palette = "Set1") + scale_y_continuous(limits = c(0, 200)) + ggtitle('Grey Squirrels (Sciurus carolinensis)')+ xlab("Site") + ylab("Duration of visits (Secs)") ggplot(a4c, aes(x=dose, y=len, fill=supp)) + geom_bar(position=position_dodge(), stat="identity") + geom_errorbar(aes(ymin=len-ci, ymax=len+ci), width=.2, # Width of the error bars position=position_dodge(.9)) a4$SITE <- as.factor(a4$SITE) a9 <- subset(squirrels, SPECIES == 'Sciurus vulgaris', select=TIME.AT.FEEDER) RS.time.boxplot<- ggplot(a4, aes(y = TIME.AT.FEEDER, x = SITE, fill = TREATMENT, order = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Red Squirrel (Sciurus vulgaris) response to pine marten scent") + xlab("SITE") + ylab("DURATION OF VISITS (s)")+ scale_y_continuous(limits = c(0, 1000)) GS.time.boxplot<- ggplot(a5, aes(y = TIME.AT.FEEDER, x = SITE, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Grey Squirrel (Sciurus carolinensis) response to pine marten scent")+ xlab("SITE") + ylab("DURATION OF VISITS (s)") + scale_y_continuous(limits = c(0, 1000)) #SGcam02 SG02fig<- ggplot(SG3, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") SGcam02$TIME.AT.FEEDER<- as.numeric(SGcam02$TIME.AT.FEEDER) t.test(TIME.AT.FEEDER~TREATMENT, data= a3) SG1 <- group_by(SGcam02, SPECIES, VISIT.NUMBER, TREATMENT) SG2 <- select(SG1, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) SG3 <- summarise(SG2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) #Number of visits a1 <- group_by(SGcam02, DATE, TREATMENT) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #SGcam 03 a1 <- group_by(SG03, DATE, TREATMENT) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) SG03fig<- ggplot(SG033, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") SG03$TIME.AT.FEEDER<- as.numeric(SG03$TIME.AT.FEEDER) t.test(TIME.AT.FEEDER~TREATMENT, data= SG033) SG031 <- group_by(SG03, SPECIES, VISIT.NUMBER, TREATMENT) SG032 <- select(SG031, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) SG033 <- summarise(SG032, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) SG033 <- SG033[-c(1:13),] #TM01cam TM3 <- TM3[-c(1:6),] TM01$TIME.AT.FEEDER<- as.numeric(TM01$TIME.AT.FEEDER) t.test(TIME.AT.FEEDER~TREATMENT, data= a3) TM1 <- group_by(TM01, SPECIES, VISIT.NUMBER, TREATMENT) TM2 <- select(TM1, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) TM3 <- summarise(TM2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) Tollymore01 <-ggplot(TM3, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Tollymore 02") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") a1 <- group_by(TM01, DATE, TREATMENT) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #TM02cam a1 <- group_by(TM02, DATE, TREATMENT) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #MS2cam a1 <- group_by(MS2, DATE, TREATMENT) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #MS3cam a1 <- group_by(MS3, DATE, TREATMENT, SPECIES) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) TM3 <- TM3[-c(1:6),] MS3$TIME.AT.FEEDER<- as.numeric(MS3$TIME.AT.FEEDER) t.test(TIME.AT.FEEDER~TREATMENT, data= MS33) MS31 <- group_by(MS3, SPECIES, VISIT.NUMBER, TREATMENT) MS32 <- select(MS31, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) MS33 <- summarise(MS32, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) MountStewart3 <-ggplot(MS33, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Mount Stewart 03") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") #Nugents01 a1 <- group_by(NW1, DATE, TREATMENT, SPECIES) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) NW1$TIME.AT.FEEDER<- as.numeric(NW1$TIME.AT.FEEDER) a3 <- a3[-c(1:3),] t.test(TIME.AT.FEEDER~TREATMENT, data= NW13) NW11 <- group_by(NW2, SPECIES, VISIT.NUMBER, TREATMENT) NW12 <- select(NW11, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) NW13 <- summarise(NW12, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) Nugents01 <- ggplot(NW13, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Nugent's Woods 02") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") #Nugents02 a1 <- group_by(NW2, DATE, TREATMENT, SPECIES) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) a3 <- a3[-c(16),] #Ballywhite a1 <- group_by(Ballywhite, DATE, TREATMENT, SPECIES) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #MS4 a1 <- group_by(MS4, DATE, TREATMENT, SPECIES) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) a3 <- a3[-c(7, 11, 13),] #Ballywalter Ballywalter$TIME.AT.FEEDER <- as.numeric(Ballywalter$TIME.AT.FEEDER) #turn time into numeric BA1 <- group_by(Ballywalter, SPECIES, VISIT.NUMBER, TREATMENT) #group by visit number BA2 <- select(BA1, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) #select other variables BA3 <- summarise(BA2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) #sum duration time of each visit t.test(TIME.AT.FEEDER~TREATMENT, data= a3) Ballywalterfig <-ggplot(BA3, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle('Ballywalter') + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") #No. visits a1 <- group_by(Ballywalter, DATE, TREATMENT, SPECIES) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) a3 <- a3[-c(7, 11, 13),] #Dunmurry No. Vists a1 <- group_by(DM, DATE, TREATMENT, SPECIES) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #Bristow no. visits a1 <- group_by(BP, DATE, TREATMENT, SPECIES) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #Garranard a1 <- group_by(G, DATE, TREATMENT, SPECIES) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #Torrosch a1 <- group_by(T, DATE, TREATMENT, SPECIES) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #Stranmillis a1 <- group_by(Stranmillis, DATE, TREATMENT, SPECIES) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #Balmoral cam Balmoral$TIME.AT.FEEDER <- as.numeric(Balmoral$TIME.AT.FEEDER) #turn time into numeric BM1 <- group_by(Balmoral, SPECIES, VISIT.NUMBER, TREATMENT) #group by visit number BM2 <- select(BM1, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) #select other variables BM3 <- summarise(BM2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) #sum duration time of each visit BM3 <- BM3[-c(1:3),] #remove badgers from database Balmoralfig <- ggplot(BM3, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Balmoral") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") t.test(TIME.AT.FEEDER~TREATMENT, data= BM3) #Balmoral visits a1 <- group_by(Balmoral, DATE, TREATMENT) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #Ballywhite cam Ballywhite$TIME.AT.FEEDER <- as.numeric(Ballywhite$TIME.AT.FEEDER) #turn time into numeric BW1 <- group_by(Ballywhite, SPECIES, VISIT.NUMBER, TREATMENT) #group by visit number BW2 <- select(BW1, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) #select other variables BW3 <- summarise(BW2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) #sum duration time of each visit t.test(TIME.AT.FEEDER~TREATMENT, data= a3) Ballywhitefig <-ggplot(BW3, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Ballywhite")+ xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") #MountStewart2 MS2$TIME.AT.FEEDER <- as.numeric(MS2$TIME.AT.FEEDER) #turn time into numeric MS21 <- group_by(MS2, SPECIES, VISIT.NUMBER, TREATMENT) #group by visit number MS22 <- select(MS21, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) #select other variables MS23 <- summarise(MS22, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) #sum duration time of each visit MS21 <- group_by(MS2, SPECIES, TREATMENT) #group by visit number MS22 <- select(MS21, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) #select other variables MS23 <- summarise(MS22, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) #sum duration time of each visit t.test(TIME.AT.FEEDER~TREATMENT, data= a3) MountStewart2 <- ggplot(MS23, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Mount Stewart 2") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") #Tollymore 02 cam a3 <- a3[-c(1:36),] TM2cam$TIME.AT.FEEDER<- as.numeric(TM2cam$TIME.AT.FEEDER) t.test(TIME.AT.FEEDER~TREATMENT, data= a3) TM021 <- group_by(TM2cam, SPECIES, VISIT.NUMBER, TREATMENT) TM022 <- select(TM021, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) TM023 <- summarise(TM022, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) Tollymore02 <-ggplot(TM023, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Tollymore 02") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") #Danesford cam Danesford$TIME.AT.FEEDER<- as.numeric(Danesford$TIME.AT.FEEDER) t.test(TIME.AT.FEEDER~TREATMENT, data= D3) D1 <- group_by(Danesford, SPECIES, VISIT.NUMBER, TREATMENT) D2 <- select(D1, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) D3 <- summarise(D2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) Danesfordfig <- ggplot(D3, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle('Danesford') + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") #Danesford no of visits a1 <- group_by(Danesford, DATE, TREATMENT) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #Stranmillis Stranmillis$TIME.AT.FEEDER<- as.numeric(Stranmillis$TIME.AT.FEEDER) t.test(TIME.AT.FEEDER~TREATMENT, data= S3) S1 <- group_by(Stranmillis, SPECIES, VISIT.NUMBER, TREATMENT) S2 <- select(S1, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) S3 <- summarise(S2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) Stranmillisfig <- ggplot(S3, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Stranmillis") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") a1 <- group_by(Stranmillis, DATE, TREATMENT) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #Garranard G$TIME.AT.FEEDER<- as.numeric(G$TIME.AT.FEEDER) t.test(TIME.AT.FEEDER~TREATMENT, data= G3) G1 <- group_by(G, SPECIES, VISIT.NUMBER, TREATMENT) G2 <- select(G1, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) G3 <- summarise(G2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) G3 <- G3[-c(1),] Gfig <- ggplot(G3, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Garranard") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") #MS4 MS4$TIME.AT.FEEDER<- as.numeric(MS4$TIME.AT.FEEDER) a3 <- a3[-c(1:5),] t.test(TIME.AT.FEEDER~TREATMENT, data= a3) MS41 <- group_by(MS4, SPECIES, VISIT.NUMBER, TREATMENT) MS42 <- select(MS41, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) MS43 <- summarise(MS42, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) MS4fig <-ggplot(MS43, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Mount Stewart 04") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") #Nugents02 NW2$TIME.AT.FEEDER<- as.numeric(NW2$TIME.AT.FEEDER) a3 <- a3[-c(1:3),] t.test(TIME.AT.FEEDER~TREATMENT, data= NW2) NW21 <- group_by(NW2, SPECIES, VISIT.NUMBER, TREATMENT) NW22 <- select(NW21, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) NW23 <- summarise(NW22, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) Nugents02 <- ggplot(NW23, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Nugent's Woods 02") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") #Lacefield Lacefield$TIME.AT.FEEDER<- as.numeric(Lacefield$TIME.AT.FEEDER) L3 <- L3[-c(1:2),] t.test(TIME.AT.FEEDER~TREATMENT, data= L3) L1 <- group_by(Lacefield, SPECIES, VISIT.NUMBER, TREATMENT) L2 <- select(L1, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) L3 <- summarise(L2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) Lacefieldfig <- ggplot(L3, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Lacefield") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") #Dunmurry DM$TIME.AT.FEEDER<- as.numeric(DM$TIME.AT.FEEDER) t.test(TIME.AT.FEEDER~TREATMENT, data= DM3) DM1 <- group_by(DM, SPECIES, VISIT.NUMBER, TREATMENT) DM2 <- select(DM1, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) DM3 <- summarise(DM2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) Dunmurry<- ggplot(DM3, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Dunmurry") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") #Nugents02 N02$TIME.AT.FEEDER<- as.numeric(N02$TIME.AT.FEEDER) a3 <- a3[-c(1:3),] t.test(TIME.AT.FEEDER~TREATMENT, data= a3) a1 <- group_by(Lacefield, DATE) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) ggplot(a3, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") #New Forge NF$TIME.AT.FEEDER<- as.numeric(NF$TIME.AT.FEEDER) a3 <- a3[-c(1:3),] t.test(TIME.AT.FEEDER~TREATMENT, data= NF3) NF1 <- group_by(NF, SPECIES, VISIT.NUMBER, TREATMENT) NF2 <- select(NF1, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) NF3 <- summarise(NF2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) NFfig <- ggplot(NF3, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("New Forge") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") a1 <- group_by(NF, DATE, TREATMENT) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #UFTM VISITS UFTM$TIME.AT.FEEDER<- as.numeric(UFTM$TIME.AT.FEEDER) a3 <- a3[-c(1:3),] t.test(TIME.AT.FEEDER~TREATMENT, data= UFTM3) UFTM1 <- group_by(UFTM, SPECIES, VISIT.NUMBER, TREATMENT) UFTM2 <- select(UFTM1, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) UFTM3 <- summarise(UFTM2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) UFTMfig <- ggplot(UFTM3, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Ulster Folk-Transport Museum") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") a1 <- group_by(UFTM, DATE, TREATMENT) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #Torrosch visits T$TIME.AT.FEEDER<- as.numeric(T$TIME.AT.FEEDER) a3 <- a3[-c(1:3),] t.test(TIME.AT.FEEDER~TREATMENT, data= T3) T1 <- group_by(T, SPECIES, VISIT.NUMBER, TREATMENT) T2 <- select(T1, DATE, SPECIES, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) T3 <- summarise(T2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) Tfig <- ggplot(T3, aes(y = TIME.AT.FEEDER, x = TREATMENT, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Torrosch") + xlab("TREATMENT") + ylab("DURATION OF VISITS (s)") a1 <- group_by(Torrosch, DATE, TREATMENT) a2 <- select(a1, DATE, VISITS, TREATMENT) a3 <- summarise(a2, VISITS = sum(VISITS, na.rm = TRUE)) #Boxplot of all durations squirrels$TIME.AT.FEEDER <- as.numeric(squirrels$TIME.AT.FEEDER) a1 <- group_by(squirrels, SPECIES, SITE, VISIT.NUMBER, TREATMENT) a2 <- select(a1, DATE, SPECIES, SITE, TREATMENT, VISIT.NUMBER, TIME.AT.FEEDER) a3 <- summarise(a2, TIME.AT.FEEDER = sum(TIME.AT.FEEDER, na.rm = TRUE)) a3 <- a3[-c(1:80),] ggplot(a3, aes(y = log(TIME.AT.FEEDER, x = SPECIES, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Red Squirrel Response to Pine Marten Scent") + xlab("TREATMENT") + ylab("NUMBER OF VISITS (per day)") ggplot(a3, aes(x = Species, y = TIME.AT.FEEDER, fill = Species)) + geom_boxplot() + facet_wrap(~variable, scales = "free") + theme_classic() mod3 <- lm(TIME.AT.FEEDER~TREATMENT+SPECIES, data = a3) summary(mod3) #Reds No. Visits #boxplots of all visit numbers Red.boxplot <- ggplot(Reds, aes(y = VISITS, x = SITE, fill = TREATMENT, order = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Red Squirrel Response to Pine Marten Scent") + xlab("SITE") + ylab("NUMBER OF VISITS (per day)")+ scale_y_continuous(limits = c(0, 400)) #scatterplots Red.scatter <-ggplot(Reds, aes(x=Day, y=VISITS, color=SITE)) + geom_point() + geom_smooth(method='loess', se=FALSE, fullrange=TRUE) + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Red Squirrel Response to Pine Marten Scent") + xlab("Day Number") + ylab("Number of Visits (per day)") ggplot(Reds, aes(x=Day, y=VISITS, color=SITE, shape=TREATMENT)) + geom_point() + geom_smooth(method=lm, se = FALSE) #bargraphs of number of visits Redsc <-summarySE(Reds, measurevar="VISITS", groupvars=c("Day")) Greysc <-summarySE(Greys, measurevar="VISITS", groupvars=c("Day")) RSbar <-ggplot(Redsc, aes(Day, VISITS, fill = "Red.Squirrels")) + geom_bar(stat="identity", position = "dodge") + geom_errorbar(aes(ymin=VISITS, ymax=VISITS+se), width=.2, # Width of the error bars position=position_dodge(.9)) + scale_fill_manual(values=c("grey", "black"))+ theme_bw() + scale_y_continuous(expand=c(0,0), limits=c(0, 100), breaks=seq(0, 100, by=10))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text.x=element_text(colour = "black", face = "bold")) + guides(fill=FALSE)+ labs(title = "A")+ xlab("Day") + ylab("Number of visits (per day)") RSbar + scale_fill_manual(values=c(Red.Squirrels="red")) GSbar <-ggplot(Greysc, aes(Day, VISITS)) + geom_bar(stat="identity", position = "dodge") + geom_errorbar(aes(ymin=VISITS-se, ymax=VISITS+se), width=.2, # Width of the error bars position=position_dodge(.9)) + xlab("Day") + ylab("Number of visits (per day)") + scale_y_continuous(limits = c(0, 150)) vistsfigs <-grid.arrange(RSbar, GSbar, nrow = 2, ncol = 1) #Greys No. Visits #boxplots of all visit numbers Grey.boxplot <- ggplot(Greys, aes(y = VISITS, x = SITE, fill = TREATMENT)) + #remove fill = site for black and white geom_boxplot() + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Grey Squirrel Response to Pine Marten Scent") + xlab("SITE") + ylab("NUMBER OF VISITS (per day)") + scale_y_continuous(limits = c(0, 400)) #scatterplot of all visit numbers ggplot(Greys, aes(x=Day, y=VISITS, color=SITE, shape=SITE)) + geom_point() + geom_smooth(method=lm) # Remove confidence intervals # Extend the regression lines Grey.scatter <-ggplot(Greys, aes(x=Day, y=VISITS, color=SITE)) + geom_point() + geom_smooth(method='loess', se=FALSE, fullrange=TRUE) + geom_jitter(position = position_jitter(width = 0.05, height = NULL)) + #remove this line if oyu don't want the points of each plotted ggtitle("Grey Squirrel Response to Pine Marten Scent") + xlab("Day Number") + ylab("Number of Visits (per day)") grid.arrange(Red.boxplot, Grey.boxplot, ncol = 1) grid.arrange(Grey.scatter, Red.scatter, ncol = 1) grid.arrange(redsbar, greysbar, ncol = 1) grid.arrange(SG02fig, SG03fig, Tolly more01, Tollymore02, MountStewart2, MountStewart3, Nugents01, Nugents02, Ballywhitefig, MS4fig, Ballywalterfig, ncol = 2) grid.arrange(RS.time.boxplot, GS.time.boxplot, ncol = 1) grid.arrange(RSbar, GSbar, ncol = 1) grid.arrange #cheeky linear model mod1 <- lm(VISITS~TREATMENT, data = Greys) summary(mod1) mod2<- lm(VISITS~TREATMENT, data= Reds) summary(mod2) mod3<- lmer(VISITS~TREATMENT + (1|SITE) + (1|Day), data = Greys) summary(mod3) mod4<- lmer(VISITS~0 + (1|SITE) + (1|Day), data = Greys) anova(mod3, mod4) mod5<- lmer(VISITS~TREATMENT + (1|SITE) + (1|Day), data = Reds) summary(mod5) mod6<- lmer(VISITS~0 + (1|SITE) + (1|Day), data = Reds) anova(mod5) lrtest(mod5,mod6) mod5<- lmer(TIME.AT.FEEDER~TREATMENT + (1|SITE), data = a4) mod1 <- lmer(VISITS~SPECIES*TREATMENT*Day + (1|SITE), data = visits) summary (mod1) mod2<- aov(mod1) #lme library(nlme) m1 <- lme(VISITS~TREATMENT,random=~1|SITE, method = 'REML', data=Greys) anova(m1) m2 <- lme(VISITS~TREATMENT,random=~1|SITE, method = 'REML', data=Reds) anova(m2) m3<- lme(TIME.AT.FEEDER~TREATMENT,random=~1|SITE, data=a4) m5<- lm(TIME.AT.FEEDER~TREATMENT, data=a4) summary(m3) anova(m3,m5) m3<- lmer(TIME.AT.FEEDER~TREATMENT+SPECIES+SPECIES:TREATMENT + (1|SITE), data=a3) a3$TIME.AT.FEEDER <- log(a3$TIME.AT.FEEDER) m4<- glmer(TIME.AT.FEEDER~TREATMENT+SPECIES+SPECIES:TREATMENT + (1|SITE), family = poisson, data=a3) summary(m4) anova(m3, m4) #model 4 has lower AIC value a4 <- a3 a4$TIME.AT.FEEDER <- log(a4$TIME.AT.FEEDER) m4<- glmer(TIME.AT.FEEDER~TREATMENT+SPECIES+SPECIES:TREATMENT + (1|SITE), family = Gamma, link = 'log', data=a3) m5<- glmer(TIME.AT.FEEDER~TREATMENT+SPECIES+SPECIES:TREATMENT + (1|SITE), family = Gamma, data=a3) anova(m4, m5) library('asreml') reml1 <- asreml(fixed = VISITS~TREATMENT,random=~SITE, family.asreml = poisson(identity) , data=Greys) anova(reml1) reml2 <- asreml(fixed = VISITS~TREATMENT,random=~SITE, family.asreml = poisson(identity) , data=Reds) anova(reml2) reml3<- asreml(fixed = TIME.AT.FEEDER~TREATMENT,random=~SITE, family.asreml = inverse.gaussian(identify), data=a4) anova(reml3) reml4<- asreml(fixed = TIME.AT.FEEDER~TREATMENT,random=~SITE, family.asreml = inverse.gaussian(identify), data=a5) anova(reml4) #Exploring the data hist(Greys$VISITS) hist(a4$TIME.AT.FEEDER) hist(a5$TIME.AT.FEEDER) hist(Reds$VISITS) setwd('C:/R') catdata <- read.csv('distance moved.csv') hist(log(NF3$TIME.AT.FEEDER)) #Checking distribution qualtatively #Histogram first hist(a4$TIME.AT.FEEDER) #you are looking for a bell-shaped curve #Then qqplots qqnorm(a4$TIME.AT.FEEDER) #looking for a roughly diagonal line #Checking distribution quantatiatively #shapiro test shapiro.test(a4$TIME.AT.FEEDER) #significant #Transforming non normally distributed data a4$TIME.AT.FEEDER <- log(a4$TIME.AT.FEEDER) qqnorm(a4$TIME.AT.FEEDER) a4$TIME.AT.FEEDER<- log10(a4$TIME.AT.FEEDER) a4$TIME.AT.FEEDER<-exp(a4$TIME.AT.FEEDER) a4$TIME.AT.FEEDER<- log(a4$TIME.AT.FEEDER)/(1-a4$TIME.AT.FEEDER) #tests involving continuous variables, normal distribution = t.test, non normal = wilcox.test #Modelling mod1 <- lm(TIME.AT.FEEDER~TREATMENT + SPECIES, data = a4) t.test(TIME.AT.FEEDER~TREATMENT, data = a4) mod2<- lmer(TIME.AT.FEEDER~TREATMENT * SPECIES + (1|SITE), data= a4) summary(mod1) ggplot(a4, aes(y = TIME.AT.FEEDER, x = SPECIES)) + #remove fill = site for black and white geom_jitter(aes(colour = TREATMENT)) #Model selection catmod2<- update(catmod1,~.-Community:Sex:Multicat) summary(catmod2) catmod3<-update(catmod2,~.-Sex:Multicat) summary(catmod3) #define summarySE function summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE, conf.interval=.95, .drop=TRUE) { library(plyr) # New version of length which can handle NA's: if na.rm==T, don't count them length2 <- function (x, na.rm=FALSE) { if (na.rm) sum(!is.na(x)) else length(x) } # This does the summary. For each group's data frame, return a vector with # N, mean, and sd datac <- ddply(data, groupvars, .drop=.drop, .fun = function(xx, col) { c(N = length2(xx[[col]], na.rm=na.rm), mean = mean (xx[[col]], na.rm=na.rm), sd = sd (xx[[col]], na.rm=na.rm) ) }, measurevar ) # Rename the "mean" column datac <- rename(datac, c("mean" = measurevar)) datac$se <- datac$sd / sqrt(datac$N) # Calculate standard error of the mean # Confidence interval multiplier for standard error # Calculate t-statistic for confidence interval: # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1 ciMult <- qt(conf.interval/2 + .5, datac$N-1) datac$ci <- datac$se * ciMult return(datac) } g_legend<-function(a.gplot){ tmp <- ggplot_gtable(ggplot_build(a.gplot)) leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box") legend <- tmp$grobs[[leg]] return(legend)} #Squirrel stats redone rs <- subset(behaviour, SPECIES == "Red Squirrel") gs <- subset(behaviour, SPECIES == "Grey Squirrel") vgs <-subset(visits, SPECIES == "Grey Squirrel") vrs <-subset(visits, SPECIES == "Red Squirrel") mod1<- lme(VISITS~SPECIES*TREATMENT, random=~1|SITE, data = visits) summary(mod1) anova(mod1) mod1<- lme(VIGILENCE~SPECIES*TREATMENT, random=~1|SITE, data = behaviour) anova(mod1) mod2<- lme(FEEDING~SPECIES*TREATMENT, random=~1|SITE, data = behaviour) anova(mod2) mod3 <- lme(TIME.AT.FEEDER~SPECIES*TREATMENT, random=~1|SITE, data = a3) anova(mod3) mod1<- lme(VIGILENCE~TREATMENT, random=~1|SITE, data = gs) anova(mod1) mod1<- lme(FEEDING~TREATMENT, random=~1|SITE, data = rs) anova(mod1) mod1<- lme(VIGILENCE~TREATMENT, random=~1|SITE, data = rs) anova(mod1) mod1<- lme(MOVEMENT~TREATMENT, random=~1|SITE, data = rs) anova(mod1) mod1<- lme(AGRESSION~TREATMENT, random=~1|SITE, data = rs) anova(mod1) mod1<- lme(SOCIAL~TREATMENT, random=~1|SITE, data = rs) anova(mod1) mod1<- lme(INVESTIGATORY~TREATMENT, random=~1|SITE, data = rs) anova(mod1) mod1<- lme(FEEDING~TREATMENT, random=~1|SITE, data = gs) anova(mod1) mod1<- lme(VIGILENCE~TREATMENT, random=~1|SITE, data = gs) anova(mod1) mod1<- lme(MOVEMENT~TREATMENT, random=~1|SITE, data = gs) anova(mod1) mod1<- lme(AGRESSION~TREATMENT, random=~1|SITE, data = gs) anova(mod1) mod1<- lme(SOCIAL~TREATMENT, random=~1|SITE, data = gs) anova(mod1) mod1<- lme(INVESTIGATORY~TREATMENT, random=~1|SITE, data = gs) anova(mod1) mod1<- lme(VISITS~TREATMENT, random=~1|SITE, data = vrs) anova(mod1) visits<- read.csv("SquirrelsNoVisits3daysPrePost.csv") VISITSc<-summarySE(visits, measurevar="VISITS", groupvars=c("SPECIES", "DAY")) Visbar<- ggplot(VISITSc, aes(DAY, VISITS, fill = SPECIES, width = .5)) + xlab("Species") + ylab("Number of vists per day")+ geom_bar(stat="identity", position=position_dodge(width=.5), colour = 'black') + geom_errorbar(aes(ymin=VISITS, ymax=VISITS+se), width=.2, # Width of the error bars position=position_dodge(.5)) + scale_fill_manual(values=c("grey", "black"))+ geom_vline(xintercept = 3.5, size = 3)+ theme_bw() + scale_y_continuous(expand=c(0,0), limits=c(0, 120), breaks=seq(0, 120, by=20))+ scale_x_discrete(name = "Day", limits = c("1", "2", "3", "4", "5", "6")) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text.x=element_text(colour = "black", face = "bold"))+guides(fill=FALSE)