#running R version 3.5.1 #Code written by Dani Crain, ddiancrain@gmail.com library(tidyverse) library(nlme) ICprogpreg=read.csv("IceSealClaws_prog_updatepregstats_W0only_noNAs.csv", check.names=F) model_abs = aov(FinalProg ~ PregnancyStat, data=ICprogpreg) anova(model_abs) #Analysis of Variance Table #Response: FinalProg #Df Sum Sq Mean Sq F value Pr(>F) #PregnancyStat 3 18054 6018.1 3.809 0.02041 * # Residuals 29 45819 1580.0 ggplot(ICprogpreg, aes(x=PregnancyStat, y=FinalProg))+ geom_boxplot() model_Z = aov(ProgZscores ~ PregnancyStat, data=ICprogpreg) anova(model_Z) #Analysis of Variance Table #Response: ProgZscores #Df Sum Sq Mean Sq F value Pr(>F) #PregnancyStat 3 52.70 17.5665 6.7163 0.001404 ** # Residuals 29 75.85 2.6155 TukeyHSD(model_abs) #pp-dia_ov p = 0.060, preg vs pp p = 0.051 TukeyHSD(model_Z) #nonpreg vs dia_ov significanly dif, pp vs nonpreg significantly different plot(model_Z, 1) #install.packages("car") library(car) leveneTest(ProgZscores ~ PregnancyStat, data = ICprogpreg) #Levene's Test for Homogeneity of Variance (center = median) # Df F value Pr(>F) #group 3 0.8722 0.4668 # 29 png("ReproductiveStatus_iceseals.png", width = 6, height = 4, units = 'in', res = 300) ggplot(ICprogpreg, aes(x=PregnancyStat, y=ProgZscores))+ geom_boxplot()+ geom_jitter(shape=21, position=position_jitter(0.2),size=2, aes(fill = Species))+ scale_x_discrete(limits=c("nonpreg", "dia_ov", "preg","pp"), labels=c("Non-pregnant","Unimplanted","Implanted","Post-partum"))+ theme_classic()+ scale_fill_manual(values=c("gray","black"), labels=c("Bearded seal","Ringed seal"))+ xlab("Reproductive Status")+ ylab("Progesterone Z-scores") dev.off() ICprogseasons=read.csv("IceSealClaws_prog_updatepregstats_W0only_Season_noNAs.csv", check.names=F) model_season = aov(ProgZscores ~ Season, data=ICprogseasons) anova(model_season) #Analysis of Variance Table #Response: ProgZscores #Df Sum Sq Mean Sq F value Pr(>F) #Season 3 35.215 11.7382 3.8932 0.01768 * # Residuals 32 96.482 3.0151 TukeyHSD(model_season) #spring vs fall p = 0.045, summer vs fall p = 0.024 png("ReproductiveStatus_Seasons_iceseals.png", width = 6, height = 4, units = 'in', res = 300) ggplot(ICprogseasons, aes(x=Season, y=ProgZscores))+ geom_boxplot()+ geom_jitter(shape=21, position=position_jitter(0.2),size=2, aes(fill = Species))+ scale_x_discrete(limits=c("Fall", "Winter", "Spring", "Summer"))+ scale_fill_manual(values=c("gray","black"), labels=c("Bearded seal","Ringed seal"))+ theme_classic()+ xlab("Season")+ ylab("Progesterone Z-scores") dev.off() #Stable Isotopes as compared to pregnancy status EB=filter(ICprogpreg, Species=="bearded seal") PH=filter(ICprogpreg, Species=="ringed seal") png("ReproductiveStatus_iceseals_carbonSI_EB.png", width = 6, height = 4, units = 'in', res = 300) ggplot(ICprogpreg, aes(x=PregnancyStat, y=d13c.Suesscor))+ geom_boxplot()+ geom_jitter(shape=21, position=position_jitter(0.2),size=2, aes(fill = Species))+ scale_x_discrete(limits=c("nonpreg", "dia_ov", "preg","pp"), labels=c("Non-pregnant","Unimplanted","Implanted","Post-partum"))+ theme_classic()+ scale_fill_manual(values=c("gray","black"))+ labs(x="Reproductive Status", y=expression(paste(delta^{13}, "C"))) dev.off() model_Z_SId13c = aov(d13c.Suesscor ~ PregnancyStat, data=ICprogpreg) anova(model_Z_SId13c) #Analysis of Variance Table #Response: d13c.Suesscor #Df Sum Sq Mean Sq F value Pr(>F) #PregnancyStat 3 32.133 10.7111 20.136 2.988e-07 *** # Residuals 29 15.426 0.5319 TukeyHSD(model_Z_SId13c) #nonpreg to dia_ov, pp to nonpreg, preg to nonpreg, preg to pp #Tukey multiple comparisons of means #95% family-wise confidence level #Fit: aov(formula = d13c.Suesscor ~ PregnancyStat, data = ICprogpreg) #$`PregnancyStat` #diff lwr upr p adj #nonpreg-dia_ov -2.099833 -3.3030844 -0.8965823 0.0002793 #pp-dia_ov 0.572750 -0.4849667 1.6304667 0.4648095 #preg-dia_ov -0.923000 -2.0113815 0.1653815 0.1189464 #pp-nonpreg 2.672583 1.6790315 3.6661352 0.0000003 #preg-nonpreg 1.176833 0.1506974 2.2029693 0.0198335 #preg-pp -1.495750 -2.3465770 -0.6449230 0.0002537 png("ReproductiveStatus_iceseals_nitrogenSI_EB.png", width = 6, height = 4, units = 'in', res = 300) ggplot(ICprogpreg, aes(x=PregnancyStat, y=d15N))+ geom_boxplot()+ geom_jitter(shape=21, position=position_jitter(0.2),size=2, aes(fill = Species))+ scale_x_discrete(limits=c("nonpreg", "dia_ov", "preg","pp"), labels=c("Non-pregnant","Unimplanted","Implanted","Post-partum"))+ theme_classic()+ scale_fill_manual(values=c("gray","black"))+ labs(x="Reproductive Status", y=expression(paste(delta^{15}, "N"))) dev.off() EB=filter(ICprogpreg,Species=="bearded seal") PH=filter(ICprogpreg,Species=="ringed seal") EB_model_abs = aov(FinalProg ~ PregnancyStat, data=EB) anova(model) summary(model) ggplot(EB, aes(x=PregnancyStat, y=FinalProg))+ geom_boxplot() EB_model_Z = aov(ProgZscores ~ PregnancyStat, data=EB) anova(model) summary(model) plot(EB_model_Z, 1) #install.packages("car") #library(car) leveneTest(ProgZscores ~ PregnancyStat, data = EB) ggplot(EB, aes(x=PregnancyStat, y=ProgZscores))+ geom_boxplot() TukeyHSD(EB_model_abs) #pp vs dia_ov sig diff, preg vs dia_ov sig diff TukeyHSD(EB_model_Z) #no sig diff PH_model_abs = aov(FinalProg ~ PregnancyStat, data=PH) anova(model) summary(model) ggplot(PH, aes(x=PregnancyStat, y=FinalProg))+ geom_boxplot() PH_model_Z = aov(ProgZscores ~ PregnancyStat, data=PH) anova(model) summary(model) plot(PH_model_Z, 1) #install.packages("car") #library(car) leveneTest(ProgZscores ~ PregnancyStat, data = PH) ggplot(PH, aes(x=PregnancyStat, y=ProgZscores))+ geom_boxplot() TukeyHSD(PH_model_abs) #preg vs non preg sif diff TukeyHSD(PH_model_Z) #no sig diff