## ### the following r script adapted from SAS code originally developed by Laura McFarlane Tranquilla, published as Appendix 1 to Marine Ornithology 48: 263–272 (2020) ######################################################################## library(lubridate) # to combine date/times in POSIX format library(summarytools) library(stringr) # to right subset variables library(fitdistrplus) library(psych) # general stats such as means and SEs library(doBy) # table summary stats by year - also produce means for GGplot at end library(dplyr) # create lagged variables for calculating foraging trip length library(lme4) library(car) # for Anova function (capital A) library(jtools) # parameter estimation such as summ() library(lsmeans) # for post-hoc comparisons library(lmerTest) # to estimate DF for lmer but nothing works for glmer apparently library(lmtest) # for lrtest and inclusion of individual as a random variable library(ggplot2) library(ggstance) #horizontal boxplot extension to GGPLOT2 library(gtable) # for ggplot grid.newpage and grid.draw functions to combine plots library(grid) #same as above # 15-minute categories for Triangle DCC2000arr_dep files # data input setwd("C:/your_pathway") obs<-read.csv("CaauData.csv") ### ARRIVAL TIMES processing ####################################### # round arrival minutes up to 15 minutes increments obs$arrMn<-ifelse(obs$min_arr>14, ifelse(obs$min_arr<30,.25, ifelse(obs$min_arr<45,.5,.75) ) , .0) dfSummary(obs$arrMn);dfSummary(obs$min_arr) obs$hr_arr<-as.numeric(obs$hr_arr) # View(obs) obs$arrTm<-obs$hr_arr+obs$arrMn dfSummary(obs$arrTm) # # arrival time re-centre at midnight and +/-12 #first isolate outliers and take a look at whats going on then redo arrHR above # there is a problem with the extraction of hrs and minutes from arrival timee ### problem solved and verified - hrs/mins check out new var - arrHr,arrMn, and arrTm created ## arrival time manipulation to linear scale - obs$zArrTm<-ifelse(obs$arrTm > 12, obs$arrTm-12, obs$arrTm+12) dfSummary(obs$arrTm); dfSummary(obs$zArrTm) # check for consistency # View(obs) #normal<-fitdist(obs$zArrTm, "norm") #gamma <- fitdist(obs$zArrTm, "gamma",method="mge", gof="ADL") #lnorm<-fitdist(obs$zArrTm, "lnorm") #weibull<-fitdist(obs$zArrTm, "weibull") #plot(normal); title(main="Normal", cex.main=2) #plot(gamma); title(main="Gamma", cex.main=2) #plot(lnorm); title(main="LogNormal", cex.main=2) #plot(weibull); title(main="weibull", cex.main=2)## ## gamma fits best with data zeroed at noon but ## with a really sharp peak and long tail that's truncated on both ## ends by 12 ### DEPARTURE TIMES processing ####################################### # round arrival minutes up to 15 minutes increments obs$depMn<-ifelse(obs$min_dep>14, ifelse(obs$min_dep<30,.25, ifelse(obs$min_dep<45,.5,.75) ) , .0) dfSummary(obs$depMn);dfSummary(obs$min_dep) obs$hr_dep<-as.numeric(obs$hr_dep) obs$depTm<-obs$hr_dep+obs$depMn dfSummary(obs$depTm) # probably similar problems with distribution fits here, but scaling similarly to arrival times for ease of comparison/graphs ## departure time manipulation to linear scale - obs$zdepTm<-ifelse(obs$depTm > 12, obs$depTm-12, obs$depTm+12) dfSummary(obs$depTm); dfSummary(obs$zdepTm) # check for consistency #View(obs) ########################################## SUMMARY STATS FOR PUB ######################################################################### a<-as.character(obs$arrTm); b<-as.character(obs$depTm) table(a); table(b) describe(obs$zArrTm); describe(obs$zdepTm) # this provides general stats overall use zatt and zdep because of zero hours describeBy(obs$zArrTm, obs$year); describeBy(obs$zdepTm, obs$year) # stats by year # summaryBy(arrTm + depTm ~ year, data=obs, FUN=c(m="mean",s="sd") ) # this is from the doBy package nrow(obs[c(obs$zArrTm>11 & obs$zArrTm <12.25),])/1625 # obs between 2315 and 0015 - check with summary(arrTm) nrow(obs[c(obs$zdepTm>15 & obs$zdepTm <16.25),])/1625 # obs between 315 and 415 - check with summary(depTm) ########################################################### ################################################### ######## trip/attendance duration data processing ################################################### # # create date time composite variable for attendance/foraging trip length obs$trip_dur<-as.numeric(obs$trip_dur) obs$ArrDt<-ifelse(obs$year=="1999",as.Date(obs$jdat_arr, origin="1999-01-01"), ifelse(obs$year=="2000",as.Date(obs$jdat_arr, origin="2000-01-01"), as.Date(obs$jdat_arr, origin="2001-01-01") ) ) obs$ArrDt<-as.Date(obs$ArrDt, origin=as.Date("1970-01-01"), format='%Y-%m-%d') obs$DepDt<-ifelse(obs$year=="1999",as.Date(obs$jdat_dep, origin="1999-01-01"), ifelse(obs$year=="2000",as.Date(obs$jdat_dep, origin="2000-01-01"), as.Date(obs$jdat_dep, origin="2001-01-01") ) ) obs$DepDt<-as.Date(obs$DepDt, origin=as.Date("1970-01-01"), format='%Y-%m-%d') # create datetime variables ########################## obs$ArrDtime <- with(obs, as.POSIXct(paste(obs$ArrDt, paste(obs$hr_arr,":",obs$min_arr,":00",sep="")), format="%Y-%m-%d %H:%M:%S")) obs$DepDtime <- with(obs, as.POSIXct(paste(obs$DepDt, paste(obs$hr_dep,":",obs$min_dep,":00",sep="")), format="%Y-%m-%d %H:%M:%S")) obsDT<-obs[,c("birdID","sex", "year","ArrDtime","DepDtime","shift_dur","trip_dur")] sort<-obsDT[order(obsDT$birdID,obsDT$ArrDtime),] sort$l.birdID<-lead(sort$birdID) sort$l.arrDtime<-as.POSIXct(lead(sort$ArrDtime)) obsDT<-sort obsDT$forg<-ifelse(obsDT$birdID==obsDT$l.birdID, difftime(obsDT$l.arrDtime, obsDT$DepDtime, units="hours"), 99999) #obsDT$forg<-replace(obsDT$forg,obsDT$forg==99999, NA) obsDT$forg<-as.numeric(obsDT$forg) obsDT$att<-as.numeric(difftime(obsDT$DepDtime, obsDT$ArrDtime, units="hours")) check<-obsDT[abs(obsDT$trip_dur-obsDT$forg)>0.1,] #outputs lines with all NAs when trip_dur =NA -looks like all these lines are due to NA in trip_dur check2<-obsDT[abs(obsDT$shift_dur-obsDT$att)>0.1,] # looks like everything checks out - #rem NAs for the following dist because fitdistr can't handle them there is a way around - make sure last line eliminated ... obsDT<-obsDT[obsDT$forg<9000,] # fitdist cant handle NA as above so do this way then reinsert NA obsDT<-obsDT[!is.na(obsDT$l.arrDtime),] # this eliminates last line - artefact of lag ## check for distributions again - #normal<-fitdist(obsDT$forg, "norm") #gamma <- fitdist(obsDT$forg, "gamma",method="mge", gof="ADL") #lnorm<-fitdist(obsDT$forg, "lnorm") #weibull<-fitdist(obsDT$forg, "weibull") #plot(normal); title(main="Normal", cex.main=2) #plot(gamma); title(main="Gamma", cex.main=2) #plot(lnorm); title(main="LogNormal", cex.main=2) #plot(weibull); title(main="weibull", cex.main=2)## #normal<-fitdist(obsDT$att, "norm") #gamma <- fitdist(obsDT$att, "gamma",method="mge", gof="ADL") #lnorm<-fitdist(obsDT$att, "lnorm") #weibull<-fitdist(obsDT$att, "weibull") #plot(normal); title(main="Normal", cex.main=2) #plot(gamma); title(main="Gamma", cex.main=2) #plot(lnorm); title(main="LogNormal", cex.main=2) #plot(weibull); title(main="weibull", cex.main=2) table(obsDT$att) summary(obsDT$att) # there are 63 zeros? ## checked with original data and this is problem at original dataset level - removed zeros<-obsDT[obsDT$att==0,] ## zeros for attendance analyses obsDT<-obsDT[obsDT$att>0,] gt24<-obsDT[obsDT$forg>25,] table(gt24$year) ### remove wild outliers ... obsDT<-obsDT[obsDT$att>0 & obsDT$forg<25,] # # redo fit distribution above ... # remove NAS for fitdistrplus above #normal<-fitdist(obsDT$forg, "norm") #gamma <- fitdist(obsDT$forg, "gamma",method="mge", gof="ADL") #lnorm<-fitdist(obsDT$forg, "lnorm") #weibull<-fitdist(obsDT$forg, "weibull") #plot(normal); title(main="Normal", cex.main=2) #plot(gamma); title(main="Gamma", cex.main=2) #plot(lnorm); title(main="LogNormal", cex.main=2) #plot(weibull); title(main="weibull", cex.main=2)## gamma best again #normal<-fitdist(obsDT$att, "norm") #gamma <- fitdist(obsDT$att, "gamma",method="mge", gof="ADL") #lnorm<-fitdist(obsDT$att, "lnorm") #weibull<-fitdist(obsDT$att, "weibull") #plot(normal); title(main="Normal", cex.main=2) #plot(gamma); title(main="Gamma", cex.main=2) #plot(weibull); title(main="weibull", cex.main=2)## gamma best again # ## redo obsDT file from sort step above - remove zeros and outliers as above, then do analyses below with NAs included # obsDT<-sort obsDT$year<-as.character(obsDT$year) obsDT$forg<-ifelse(obsDT$birdID==obsDT$l.birdID, difftime(obsDT$l.arrDtime, obsDT$DepDtime, units="hours"), 99999) obsDT$forg<-replace(obsDT$forg,obsDT$forg==99999, NA) obsDT$forg<-as.numeric(obsDT$forg) obsDT$att<-as.numeric(difftime(obsDT$DepDtime, obsDT$ArrDtime, units="hours")) check<-obsDT[abs(obsDT$trip_dur-obsDT$forg)>0.1,] #outputs lines with all NAs when trip_dur =NA -looks like all these lines are due to NA in trip_dur check2<-obsDT[abs(obsDT$shift_dur-obsDT$att)>0.1,] # looks like everything checks out - obsDT<-obsDT[obsDT$att>0, ] # there should not be any attedance times of zero - there should be some max time # & obsDT$forg<25,] # ie. likelihood that there some missed detections while arriving? obsDT2<-obsDT[obsDT$forg<55, ] # drop extremes summary(obsDT2$forg) ########################################## SUMMARY STATS FOR PUB ################## summary(obsDT2$forg);summary(obsDT$att) # with reduced data removing zeros (maybe consider removing some of the other outliers) # summary(obs$trip_dur);summary(obs$shift_dur) describe(obsDT2$forg);describe(obsDT$att) # general stats describeBy(obsDT$att,obsDT$year); describeBy(obsDT2$forg,obsDT2$year) # stats year nrow(obsDT[obsDT$att>3 & obsDT$att<5.25,])/1625; nrow(obsDT2[obsDT2$forg>19 & obsDT2$forg<21.25,])/1625 ################################################################################### ####################################################################### ### analyses ###### ######################################################################## ## all models are nearly fully saturated and may not represent the optimal parsimonious model (as published), which was achieved through # backwards step-wise reduction # Arrival time ######################## obs$year<-as.character(obs$year) obs$day<-as.character(obs$jdat_arr) # character variable to look at within year variation glmer.Arr<-glmer(zArrTm ~ year + jdat_arr + sex + jdat_arr:sex + jdat_arr:year + (1|year/birdID), data=obs, family=Gamma, na.action = na.exclude) # + sex:year #summ(glmer.Arr, scale=TRUE) #fixef(glmer.Arr,type="response") #summary(glmer.Arr) Anova(glmer.Arr,type='III') dfSummary(obs$jdat_arr)# find best jdate to be "intercept" to estimate emm means below at that date emm.ArrY<-emmeans(glmer.Arr, ~ jdat_arr:year, at=list(jdat_arr=145), transform="response") %>% summary(infer=TRUE); emT.ArrY<-emtrends(glmer.Arr, "year", var="jdat_arr", transform="response") %>% summary(infer=TRUE); emm.ArrS<-emmeans(glmer.Arr, ~ jdat_arr:sex, at=list(jdat_arr=145), transform="response") %>% summary(infer=TRUE); emT.ArrS<-emtrends(glmer.Arr, "sex", var="jdat_arr", transform="response") %>% summary(infer=TRUE); emm.ArrY; emT.ArrY emm.ArrS; emT.ArrS # pairwise difference at beginning of study emm.ArrY<-emmeans(glmer.Arr, ~ jdat_arr:year, at=list(jdat_arr=145), transform="response") # pairs doesn't work with pipe emT and emm files - redo these before running pairs emT.ArrY<-emtrends(glmer.Arr, "year", var="jdat_arr", transform="response") pairs(emm.ArrY); pairs(emT.ArrY) # pairwise difference at end of study emm.ArrY<-emmeans(glmer.Arr, ~ jdat_arr:year, at=list(jdat_arr=187), transform="response") # pairs doesn't work with pipe emT and emm files - redo these before running pairs pairs(emm.ArrY); emm.ArrS<-emmeans(glmer.Arr, ~ jdat_arr:sex, at=list(jdat_arr=187), transform="response") pairs(emm.ArrS) # visualize slopes emmeans(glmer.Arr, ~year, transform="response") emtrends(glmer.Arr, ~ sex, var="jdat_arr", transform="response") # test for effect of individual as random factor glm.Arr<-glmer.Arr<-glm(zArrTm ~ year + jdat_arr + sex + jdat_arr:sex + jdat_arr:year , data=obs, family=Gamma, na.action = na.exclude) glm.Arr$aic lrtest(glmer.Arr,glm.Arr) # Departure time ######################## obs$year<-as.character(obs$year) glmer.Dep<-glmer(zdepTm ~ year + jdat_dep + jdat_dep:year + (1|year/birdID), data=obs, family=Gamma(link="log"), na.action=na.exclude) # + sex:year + jdat_dep:sex + sex #summary(glmer.Dep) Anova(glmer.Dep,type='III') emm.DepY<-emmeans(glmer.Dep, ~ jdat_dep:year, at=list(jdat_dep=145), transform="response") %>% summary(infer=TRUE); emT.DepY<-emtrends(glmer.Dep, "year", var="jdat_dep", transform="response") %>% summary(infer=TRUE); emm.DepY; emT.DepY # pairwise comparisons at beginning of study emm.DepY<-emmeans(glmer.Dep, ~ jdat_dep:year, at=list(jdat_dep=145), transform="response") # pairs doesn't work with pipe emT and emm files - redo these before running pairs emT.DepY<-emtrends(glmer.Dep, "year", var="jdat_dep", transform="response") pairs(emm.DepY); pairs(emT.DepY) # pairwise comparisons at end of study emm.DepY<-emmeans(glmer.Dep, ~ jdat_dep:year, at=list(jdat_dep=187), transform="response") # pairs doesn't work with pipe emT and emm files - redo these before running pairs pairs(emm.DepY) # Nest Attendance ######################## obsDT$jdatA<-yday(obsDT$ArrDtime) obsDT$jdatD<-yday(obsDT$DepDtime) # create a new jdate for obsDT glmer.att<-glmer(att ~ year + sex + jdatA +jdatA:year + jdatA:sex + (1|year/birdID), data=obsDT, family=Gamma(link="inverse")) # + sex:year summary(glmer.att) Anova(glmer.att,type='III') emm.attY<-emmeans(glmer.att, ~ jdatA:year, at=list(jdatA=145), transform="response") %>% summary(infer=TRUE); emT.attY<-emtrends(glmer.att, "year", var="jdatA", transform="response") %>% summary(infer=TRUE); emm.attS<-emmeans(glmer.att, ~ jdatA:sex, at=list(jdatA=145), transform="response") %>% summary(infer=TRUE); emT.attS<-emtrends(glmer.att, "sex", var="jdatA", transform="response") %>% summary(infer=TRUE); emm.attY; emT.attY emm.attS; emT.attS # pairwise comparison at beginning of year emm.attY<-emmeans(glmer.att, ~ jdatA:year, at=list(jdatA=145), transform="response") # pairs doesn't work with pipe emT and emm files - redo these before running pairs emT.attY<-emtrends(glmer.att, "year", var="jdatA", transform="response") pairs(emm.attY); pairs(emT.attY) # pairwise comparison at end of year emm.attY<-emmeans(glmer.att, ~ jdatA:year, at=list(jdatA=187), transform="response") emm.attS<-emmeans(glmer.att, ~ jdatA:sex, at=list(jdatA=187), transform="response") pairs(emm.attY); pairs(emm.attS) # foraging trip duration ################### glmer.forg<-glmer(forg ~ sex + jdatD + year + jdatD:sex + jdatD:year + (1|year/birdID), data=obsDT, family=Gamma(link="inverse")) # + sex:year summary(glmer.forg) Anova(glmer.forg,type='III') glm.forg<-glm(forg ~ year , data=obsDT, family=Gamma(link="inverse")) # + sex:year glm.forg$aic lrtest(glmer.forg,glm.forg) emm.forgY<-emmeans(glmer.forg, ~ jdatD:year, at=list(jdatD=145), transform="response") %>% summary(infer=TRUE); emT.forgY<-emtrends(glmer.forg, "year", var="jdatD", transform="response") %>% summary(infer=TRUE); emm.forgS<-emmeans(glmer.forg, ~ jdatD:sex, at=list(jdatD=145), transform="response") %>% summary(infer=TRUE); emT.forgS<-emtrends(glmer.forg, "sex", var="jdatD", transform="response") %>% summary(infer=TRUE); emm.forgY; emT.forgY emm.forgS; emT.forgS # pairwise comparisons at beginning of study period emm.forgY<-emmeans(glmer.forg, ~ jdatD:year, at=list(jdatD=145), transform="response") # pairs doesn't work with pipe emT and emm files - redo these before running pairs emT.forgY<-emtrends(glmer.forg, "year", var="jdatD", transform="response") pairs(emm.forgY); pairs(emT.forgY) # pairwise comparisons at end of study period emm.forgY<-emmeans(glmer.forg, ~ jdatD:year, at=list(jdatD=187), transform="response") # pairs doesn't work with pipe emT and emm files - redo these before running pairs emm.forgS<-emmeans(glmer.forg, ~ jdatD:sex, at=list(jdatD=187), transform="response") pairs(emm.forgY); pairs(emm.forgS) ## reporting N size for analyses tables in MS describeBy(obs$arrTm,obs$year); describeBy(obs$depTm,obs$year) # stats year describeBy(obs$arrTm,obs$sex); describeBy(obs$depTm,obs$sex) # stats sex describeBy(obsDT$att,obsDT$year); describeBy(obsDT$forg,obsDT$year) # stats year describeBy(obsDT$att,obsDT$sex); describeBy(obsDT$forg,obsDT$sex) # stats sex # figure ################# ########### # combine arrival/departure times for boxplots below at<-as.data.frame(obs$zArrTm) at$type<-"Arrival" at$year<-obs$year names(at)<-c("Time","Type","Year") dt<-as.data.frame(obs$zdepTm) dt$type<-"Departure" dt$year<-obs$year names(dt)<-c("Time","Type","Year") times<-rbind(at,dt) p<-ggplot(data = times, aes(x = Time, y = Year, fill=Type)) + geom_boxploth() + scale_y_discrete(name="Year") + scale_x_continuous(name="Time of Day (24Hr)", breaks = c(5,10,15,20), labels = c("15","20","5","10")) + scale_fill_grey(drop=FALSE, start = 0.3) + theme_bw() p ggsave("ArrivalDeparture_Times.jpg", plot = last_plot(), device = NULL, dpi = 300) # combine foraging trip and attendance times for boxplots below at<-as.data.frame(obsDT2$forg) at$type<-"Trip Duration" at$year<-obsDT2$year names(at)<-c("Time","Type","Year") dt<-as.data.frame(obsDT$att) dt$type<-"Nest Attendance" dt$year<-obsDT$year names(dt)<-c("Time","Type","Year") times<-rbind(at,dt) p<-ggplot(data = times, aes(x = Year, y = Time, fill=Type)) + geom_boxplot() + scale_y_continuous(name="Time (Hr)", breaks = c(5,10,15,20,25,30,35,40,45,50)) + scale_x_discrete(name="Year") + scale_fill_grey(drop=FALSE, start = 0.3) + theme_bw() p ggsave("Foraging_NestAttendance_Times.jpg", plot = last_plot(), device = NULL, dpi = 300) plot(p) # Figures II ###################### # combine arrival/departure times for boxplots below at<-as.data.frame(obs$zArrTm) at$type<-"Arrival" at$year<-obs$year at$sex<-obs$sex at$day<-obs$jdat_arr names(at)<-c("Time","Type","Year","Sex","Year.Day") dt<-as.data.frame(obs$zdepTm) dt$type<-"Departure" dt$year<-obs$year dt$sex<-obs$sex dt$day<-obs$jdat_arr names(dt)<-c("Time","Type","Year","Sex","Year.Day") times<-rbind(at,dt) yr<-summaryBy(Time ~ Year.Day + Year + Type, data=times, FUN=c(length,mean,sd)) yr$se<-yr$Time.sd/sqrt(yr$Time.length); yr<-yr[yr$Time.mean>6,] sx<-summaryBy(Time ~ Year.Day + Sex + Type, data=times, FUN=c(length,mean,sd)) sx$se<-sx$Time.sd/sqrt(sx$Time.length); sx<-sx[sx$Time.mean>6,] yrA<-yr[yr$Type=="Arrival",] yrD<-yr[yr$Type=="Departure",] sxA<-sx[sx$Type=="Arrival" & !is.na(sx$Sex),] sxD<-sx[sx$Type=="Departure" & !is.na(sx$Sex),] p<-ggplot(data = yrA, aes(x = Year.Day, y = Time.mean), colour=Year, fill=Year) + scale_y_continuous(name="Time of Day (24Hr)", limits=c(9, 17), breaks = c(10,12,14,16), labels = c("22","00","02","04")) + scale_x_continuous(limits=c(140,190)) + scale_colour_grey(start=.1, end=.8) + scale_fill_grey(start=0.99, end=0.99) + stat_smooth(data=yrA, aes(group=Year, linetype=Year), method=lm,formula =y~x, colour="black", size=.75) + stat_smooth(data=yrD, aes(group=Year, linetype=Year), method=lm,formula =y~x, colour="black", size=.75) + geom_errorbar(data=yrA, aes(ymin=Time.mean-se, ymax=Time.mean+se), width=0.01) + geom_point(data=yrA, aes(shape=Year), size=3) + geom_errorbar(data=yrD, aes(ymin=Time.mean-se, ymax=Time.mean+se), width=0.01) + geom_point(data=yrD, aes(shape=Year), size=3) + theme_bw() + theme(legend.key.size = unit(1, "cm"), text=element_text(size=16), axis.title.x=element_blank(), axis.text.x=element_blank(), panel.border = element_blank(), axis.line = element_line(size=1)) p # ggsave("ArrivalDeparture_TimesYrs.jpg", plot = last_plot(), device = NULL, dpi = 600, height=17.5, width=30, units="cm") q<-ggplot(data = sxA, aes(x = Year.Day, y = Time.mean), colour=Sex, fill=Sex) + scale_y_continuous(name="Time of Day (24Hr)", limits=c(9, 17), breaks = c(10,12,14,16), labels = c("22","00","02","04")) + scale_x_continuous(name="Year-Day (from 1 Jan.)", limits=c(140,190)) + scale_colour_grey(start=.1, end=.8) + scale_fill_grey(start=0.99, end=0.99) + stat_smooth(data=sxA, aes(group=Sex, linetype=Sex), method=lm,formula =y~x, colour="black", size=.75) + stat_smooth(data=sxD, aes(group=Sex, linetype=Sex), method=lm,formula =y~x, colour="black", size=.75) + geom_errorbar(data=sxA, aes(ymin=Time.mean-se, ymax=Time.mean+se), width=0.01) + geom_point(data=sxA, aes(shape=Sex), size=3) + geom_errorbar(data=sxD, aes(ymin=Time.mean-se, ymax=Time.mean+se), width=0.01) + geom_point(data=sxD, aes(shape=Sex), size=3) + theme_bw() + theme(legend.key.size = unit(1, "cm"), text=element_text(size=16), panel.border = element_blank(), axis.line = element_line(size=1)) q grid.newpage() # this step is for combining plots and aligning x-axes grid.draw(rbind(ggplotGrob(p), ggplotGrob(q), size="last")) ggsave("ArrivalDeparture_Times.jpg", grid.draw(rbind(ggplotGrob(p), ggplotGrob(q), size="last")), device = NULL, dpi = 600, height=20, width=25, units="cm")