From b46c6d741e6d1906a8307169ce3df083a5e10d6f Mon Sep 17 00:00:00 2001 From: Sam Gurr Date: Wed, 17 Jan 2024 11:46:39 -0500 Subject: [PATCH] resp and survival scripts --- RAnalysis/Output/Survival/F3_survival.pdf | Bin 7049 -> 7031 bytes RAnalysis/Scripts/.Rhistory | 1002 ++++++++--------- .../RespirationRates_analysis_larvae.Rmd | 179 ++- RAnalysis/Scripts/RespirationRates_calc.Rmd | 2 + RAnalysis/Scripts/Survival.Rmd | 114 +- 5 files changed, 761 insertions(+), 536 deletions(-) diff --git a/RAnalysis/Output/Survival/F3_survival.pdf b/RAnalysis/Output/Survival/F3_survival.pdf index cd217adb0d43103a4f7f420d479be0ad775507e0..376ccc0a359f2e4b5079e7fb289da46e9c09b64d 100644 GIT binary patch delta 3005 zcmZWlc{~&DA9qy@NsishLVZg<@UBqEs&#qr_dQ1jE~_vE#kEYWTp`T$ zv31Nn#;ShEWWoL(WoIVZ=Jr<+r<-vIO!3=W8JneDW zyW4DY02RS9(Km$u)#x9w@AhLI5oH2e< zFrmS(EodMEpb_AZg<`HdR>f4uQbkyD-VU+M)vYa zH)`q>FC~qC0vk!T!ukg!ijrKP%8MqQ>t`dBJ~cBqOXRm2>P26m3Kc|WTwPk2Hs~8y zEtz$AF~$Bv4=9BTb9UV*JW;b6kCcHTw+omaSqs1(L}UZ zNfzx99V<oU+wj^la@LF~(!({J`{P#UOSum7$1%_QWPNhQFq#_-mN-x1gHAS|;4ehs z`1%%f9yGfSoea$`^X#vO_*Yy>>eKpGx22OnzLQcjK%WtX>%`?CTu=n0M2-~LiYc19 zXax7qcCv8CXB8>IgsZlKE@0^C`=N$bM~O=0#v6 zxRuGR(GNZkrJTMeitwv)+v4zz^E1sr+ouabswefamO?aOmEzibR;87Ih^PlnI%lPS zRgfp9pY>`;M+l8M(obJi!H}IO7o#h1^Q!q;QH*s8#`9Y+H7g|gx=1PZb2K)N+TB>8 z`8E4Wt!51_m`R>TM*wQoFOX$@c7|0zme7zU{x$1(!&XGOlAT@AeIBk>yy{(aYo%$d zh5fDK3%80v#qvC*FU5mHXsB9vhf1>48C(VjVR<$w(R|DSDZEZpml+L4{`sQ&JzU4b z#%>n(IH^Zb)}r#BbCXS&?rONamY#n@){dyBRrAXmg;~=BOE>2au#_A+ExYPY7;d8D z)eXE%w@|_AP}XlOw!^9N{pweAT1UF&+=(1@+wIjC>qVCkkzL;sNZYvsvOY6d@N$^L zhf`8jnY=u4keRO1cGGEoquVi8;0on2(bpf|PbPx|H;LuzpBI1o(Vq>b2~n^dd4px3 z@)`wy)$zMRDvv$<<&zd;e1;;ITY5xG2bgJjSbtKhRBisR7jDNEm^qU3 zWhm3z24q*-yG}nBliH^0!xzddYwDg+vJ7H)C9PcQnl_k-kaH}(ydS_8m##bu^qHK- zc@0{(7272vr9fU4l~}ZaEfQ4@Ju4zwv(>TDJw;(gJnFbebC|X~-oU4Jor1l~pL_jB zNs$taP=ROd9cSsCB?QY&v5R8%`C-RZ&64yMm~Cj0m|JrvNo>k*5#pY+sE1m4q++qs z#q7Z}5#ly=WTu|z#f5R7f(vD@lo#?VMx=c&hIWmW`w~1t(c5`D4zv2V2Qj`S(wnrh zUj?Pnk@Q-AKtGiid)n|M}(X7VS3a*J9Vy_-ugJ)X{!aoKQY zYmlg-Qt%P6we%z>+)3>GTDLDExh0nyYo|e})}DR8=Qn3eq^Z;G4bj$h zm)7a4zI-o40vQ_3+0olQ?8wQWjRV3oSHKB^3mVPWPaedZ8_;6vED=(&i5@oJRQGBL zhhVc@w|EF(u{gIKw)(JU2AFhsCH(w#xKrDnUo7}4K6!3rb7b1X-aylmx1hGmi)AT5 z!)enzv3`Ss;f*J6$W{Kebfu*XWcs4I@QqxBB&8kuxu1$vGL1v;=BH1&wKqqi70d-} znprVh=iC){kcI+x%Q0WQym zMsRMY^d>yxP;%WPR#)(fu_KY_{E;Tf*esAyp2*F3D}h^za^1SutW5y#V5lOLi6qPh zQ2*Bi!0RE&Xwl)+4IbiGT9qHv^NefiYE^*tojQ`z6G!CC2^YTL@4QWC0e@N*0?x+9bVaf;VEuXR63R9XP9|6y&Ry?ESrv!b|4~=gU3OP}N z-rVjJhSQ^J**!f?;?`p@@jfx+<`+e9P5P9Ar-n;8NZv7Dv5d0(gYI9zX9gnmj|4`X z=@(Ta2j5HfiC4=Y!tH*otGp=WNYrlm%fiJ$ht@A}YM#`BKqNr_EQxy(fPW7lWw4UU z5l&ede3XMh(w&6W(NMMjmtbJ1@)1r63OdTE9i3H%s2=ets~w31gFs;a`lSAo9R^W@ Z9vOy0An6(+`T(e!GC)S=+6@E1{{UAhvWWly delta 3023 zcmZXQc{J4f8^^=7jeXzCR(8!WW@L?&WZwoE3}cruG4|zai^x!vAvAWSO+$8y5g}Y= zJT+G2n4JShN>%oE}92<8oA@V zB;}2?)WB;0@g6vTrN5+^qfI9Zn1+Zq@U(Bd&YL44&(h(;F%SB}ZU`)1Ma*Arv8lV- zYxkhk5c54aF1DGD-Bo*1f}QY&3%QKGc=mT_+3=aiF%jCEvvZpp1JiB`J}4RVEIIzz z7-Dj)gXB6E*9&yjJ(?>RBWGoaWxn&gLZkhdlZ^|}lhFfqgnp@vkuMei{Ts}%)``to zm7|zRto~p_K~?Hk_Scj@1jJqnCrjXk;|m^qG0=nuzs07fjt8Y2I5EP4rLcN8s-=Dr z?ausiLnLS@`octE*VnJfq_*Jg%fNiCc&i!qjB6k17Ql|8IW=8hdSgFLvWP>#c8MGA zoS28Zu~xb5(;mNY@tZNg&PqXuUQFIykFwnEopTv)&4l!{oL6knHGefIg!+L8UEbUz z@P?FhVslhE5Ko&6adab$Z_kV*GB^#Gv?GPm(Dw`pOi9+m7Pr?`x)qe{PpYJ(lyTMb zh^drf5*)AVs@Y@RI=m(~UjIEmsg%CK zc%Z?gB)=^LOv3Y4hx4)zMRSPH6{Sao7ljjVZ%30;`!=r59h7cI+gU%FggZ_M>zBbe z1bGiBj0XZbUB6fFU49p<^}Y3{D45n5k#&?C=6}osP6>P8#xIxtHS9M0);moxRXt35 z^JeThOn*~FH;TRFYxmNTg_MIowY~v@%1p2pmI%ZFQ!+`sP_^c~JdwO5oTYL09o#q% zzUGO=?!X>St1Z`pu8ti+@$)uI!1?dvu#VT|6z35X)&Qh7rT3>uek&i^_n=O_XzjhN zs(Hyob5buZoz-|`CCMew5G&L~d_;Sj^;t29j-FhZEZ!?E2x2X6G~R#Gf%0Nm=>$uM zaRDKcgBZ9Cc0AwtlePJ~DifU2;C>BBYUP~i8*S=V!Hp7&w1mzUrBOey(Ci=gsm~x@ zqZYXKKj6v=cT4KzxXY}*dvYC4kXilo5$~H%fW1q5!*J`y%RWYq2C?+2&D6LUY z*vrSBj5p1CT2zH7>9%jnmIJK*$KI-Os3T zv0KYrDlolBWo0WgQMP8B{=BPi0xkBf568!G`q9D%w})ROfGD}{eZux=4`8zTB!Iv? z+!MXu+C#4sqwiGQv9x3L(>AGXXw{oz>5gDT+=_~CF9^sfZ|38dqRDrTvR6|l99QpG zIck-J!E#a80>P5K)}d9`vIjwhBr%LdQM6FZO8%hGj$g~pVpmm}N`MOSazCI;RIiP5 z-FErubKqV%lXIYb&{YMLMt|6gNyotRemsme4eo2Npgi<0l0MY?YqXaX^kLUG@+*j+ z$r8DzkayM0K{X}u`A8TiC6r(ZJmiCtbs4wg{y|VURtCzy>I&{umW}buq@u+9!&)xc&$oAs=-T(JCT!q zIa(SEoyoVfpV^`Xf=07^rcDJ#bapirs3`1_zns#a9#&GsB7T1`L3OhILilqCr;k$! zvBLSO;qQna*LSIT(rpKy{a`AUH4gcy@P&Imc+an*KJ|()g4eW!AQsPv^rDk+)yifR zW#7`Vbu8{eb^{x7Y*gV?$}3WQwQ_uto#4|qK#0F`u_&9PvvYA*SEIfXzo7hdaK z{pB>Sn=`a84bmj>%=8wYaLyNX<=wa*_d_w;Tbu`%fcnEID@sf4X1K#|JO>)5LvFkB z>6_ube@z8mKHPDO$9>nPMU*^FP$ctb3R#O{ z6xmou?q)xfRS)DP%$)M@o zv~WJjIOe_yn|vn}6;;c*&FX#W)%69EPN#Ae2wXJJ@sJvasI*U|iW7MF8 zlo2WA`4evG5~=5M$=8>E^)FR(OXFqXhvtUY_XL9*AUMfn0=v`}!tbUbpBY1?^E3G@ zaQ`xsUTX0_NCeU`OAd1GvHZ`IE3qN=Pfo|xfIZ>L6$8k`Z#)+wf!6tgZ+tPs;>l^& zk3*%55-JV|_EJU)lIsq4xV>To4yw6JAtE#WE-b^6E4QN!@8;H@ZIny)WjKlow`UoC zRTALDLenqwt)93CD=iEr{^8b3Uj`a;lU*^PyON(#tUtjrV6X>co(-ZHhK%YqORgH+ z!ZaHbkk&<|@VJsgU^$$M8UD^8fINTHXPv^r4y9?p|xGPvj4XTey*FpZwk2Buy|aO_Gd<1JB`Mk1tssV zg#Gx1JtmKGY1U4?t_{)RHWEY_u(?5GHBvy6(zgjVjm|&xZP-GEt02_x=KL(?UL;

-yMi?$QuaylpLFY1fs3Zo^z9NUJC7|XPj%elm6j{4M6XWh zrZ7C8^8w5kZK1Cs2JVf@Kj~ih;+_I?9o{#3bYL)^9K35+;LCTxjzzUB@F)R)H1la| zzOIQe_N+=Bz5be4oxxmS?dUgVznoo%Kpn~j$QUZPsbQe4%tWU!q#>_)N&NWj6g651a$uVRm*Fj{{hH) BgEasE diff --git a/RAnalysis/Scripts/.Rhistory b/RAnalysis/Scripts/.Rhistory index 9d78f25..90ce9c0 100644 --- a/RAnalysis/Scripts/.Rhistory +++ b/RAnalysis/Scripts/.Rhistory @@ -1,512 +1,512 @@ -binary_df_F3_MEANS -# create surv object -# without tank included - means wihtin treatment SEPCIFICALLY FOR PLOTTING second for loop above -surv_F3_MEANS <- Surv(time = as.numeric(binary_df_F3_MEANS$Age), -event = as.numeric(binary_df_F3_MEANS$Count_dead), type = "right") -ggsurvplot(survfit(surv_F3_MEANS ~ Larvae_pCO2, binary_df_F3_MEANS), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple")), -# c(0.7,0.5, 0.3, 0.1)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS ~ Larvae_pCO2, binary_df_F3_MEANS), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple")), -c(0.7,0.5, 0.3), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS ~ Larvae_pCO2, binary_df_F3_MEANS), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple")), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS ~ Larvae_pCO2, binary_df_F3_MEANS), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple")), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -# call datafarme binary for plotting -binary_df_F3_MEANS <- data.frame() # start dataframe -binary_df_F3_MEANS_high <- data.frame() # start dataframe -binary_df_F3_MEANS_mod <- data.frame() # start dataframe -binary_df_F3_MEANS_low <- data.frame() # start dataframe -# run it -for (i in 1:nrow(df_F3_MEANS)) { -count_alive <- round(df_F3_MEANS[i,]$mean_Count_alive/1000) # divide all counts by 1000, round to nearest -count_dead <- round(df_F3_MEANS[i,]$mean_Count_dead/1000) # divide all counts by 1000, round to nearest -loopDF <- data.frame(matrix(0,ncol = 4, nrow = (count_alive + count_dead))) -colnames(loopDF) <- (c('Parent_pCO2','Larvae_pCO2','Age','Count_dead')) -loopDF <- loopDF %>% -dplyr::mutate( -Parent_pCO2 = df_F3_MEANS[i,]$Parent_pCO2, -Larvae_pCO2 = df_F3_MEANS[i,]$Larvae_pCO2, -Age = df_F3_MEANS[i,]$Age, -) -loopDF[c(1:count_dead),4] = 1 # 1 for dead -loopDF[c((count_dead + 1):nrow(loopDF)),4] = 0 # 0 for alive -loopDF <- data.frame(loopDF) # name dataframe for this single row -binary_df_F3_MEANS <- rbind(binary_df_F3_MEANS,loopDF) #bind to a cumulative list dataframe -if (loopDF$Parent_pCO2[1] == 'High pCO2') { -binary_df_F3_MEANS_high <- rbind(binary_df_F3_MEANS_high,loopDF) #bind to a cumulative list dataframe -} else if (loopDF$Parent_pCO2[1] == 'Moderate pCO2') { -binary_df_F3_MEANS_mod <- rbind(binary_df_F3_MEANS_mod,loopDF) #bind to a cumulative list dataframe -} else { -binary_df_F3_MEANS_low <- rbind(binary_df_F3_MEANS_low,loopDF) #bind to a cumulative list dataframe -} -} -binary_df_F3_MEANS_low -unique(binary_df_F3_MEANS_low$Parent_pCO2) -is.na(binary_df_F3_MEANS_low$Parent_pCO2) <- binary_df_F3_MEANS_low$Parent_pCO2 == "None" -binary_df_F3_MEANS_low$Parent_pCO2 <- factor(binary_df_F3_MEANS_low$Parent_pCO2) -str(binary_df_F3_MEANS_low) # youll see now its fixed -unique(binary_df_F3_MEANS_low$Parent_pCO2) -surv_F3_MEANS_low <- Surv(time = as.numeric(binary_df_F3_MEANS_low$Age), -event = as.numeric(binary_df_F3_MEANS_low$Count_dead), type = "right") -ggsurvplot(survfit(surv_F3_MEANS_low ~ Larvae_pCO2, binary_df_F3_MEANS_low), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple")), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -View(binary_df_F3_MEANS_low) -is.na(binary_df_F3_MEANS_mod$Parent_pCO2) <- binary_df_F3_MEANS_mod$Parent_pCO2 == "None" -binary_df_F3_MEANS_mod$Parent_pCO2 <- factor(binary_df_F3_MEANS_mod$Parent_pCO2) -str(binary_df_F3_MEANS_mod) # youll see now its fixed -surv_F3_MEANS_mod <- Surv(time = as.numeric(binary_df_F3_MEANS_mod$Age), -event = as.numeric(binary_df_F3_MEANS_mod$Count_dead), type = "right") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple")), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_low ~ Larvae_pCO2, binary_df_F3_MEANS_low), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple")), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -is.na(binary_df_F3_MEANS_high$Parent_pCO2) <- binary_df_F3_MEANS_high$Parent_pCO2 == "None" -binary_df_F3_MEANS_high$Parent_pCO2 <- factor(binary_df_F3_MEANS_high$Parent_pCO2) -str(binary_df_F3_MEANS_high) # youll see now its fixed -surv_F3_MEANS_high <- Surv(time = as.numeric(binary_df_F3_MEANS_high$Age), -event = as.numeric(binary_df_F3_MEANS_high$Count_dead), type = "right") -ggsurvplot(survfit(surv_F3_MEANS_high ~ Larvae_pCO2, binary_df_F3_MEANS_high), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple")), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_low ~ Larvae_pCO2, binary_df_F3_low), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple")), -# c(0.7,0.5, 0.3, 0.1)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_high ~ Larvae_pCO2, binary_df_F3_high), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple")), -# c(0.7,0.5, 0.3, 0.1)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_high ~ Larvae_pCO2, binary_df_F3_MEANS_high), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple")), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_low ~ Larvae_pCO2, binary_df_F3_MEANS_low), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple")), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple")), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(0.8,0.6,0.5)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(0.6,0.6,0.6)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(0.1,0.6,0.6)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.6,0.6)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.6,0.1)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.1)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.3)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_high ~ Larvae_pCO2, binary_df_F3_MEANS_high), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.3)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_low ~ Larvae_pCO2, binary_df_F3_MEANS_low), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.3)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -pdf("C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Survival/F3_low_history.pdf") -ggsurvplot(survfit(surv_F3_MEANS_low ~ Larvae_pCO2, binary_df_F3_MEANS_low), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.3)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") +width = 0.6, +size=0.4, +linetype = "dashed", +position = position_dodge(preserve = "single")) + +facet_wrap(~pCO2_parents) +F2_Dinge_boxpots +F2_Dinge_boxpots <- RR_master_F2_Dhinge %>% # omits a single dramatic outlier +dplyr::filter(!resp_umol_hr_indiv > 0.0004) %>% +ggplot(aes(x = factor(pCO2_offspring), +y = resp_umol_hr_indiv, +fill = pCO2_offspring)) + +geom_boxplot(alpha = 0.5, # color hue +width=0.6, # boxplot width +outlier.size=0, # make outliers small +position = position_dodge(preserve = "single")) + +geom_point(pch = 19, +position = position_jitterdodge(0.5), +size=1) + +scale_fill_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), +values=c("forestgreen","orange", "purple")) + +#scale_colour_manual(values=c("forestgreen","orange")) + +theme_classic() + +ggtitle("F2 Dhinge 20230407 (parentpCO2)") + +theme(axis.title.y=element_text(size=7), +axis.title.x=element_text(size=7), +axis.text.x=element_text(size=7), +legend.position = "none") + +#ylim(0, 0.2) + +stat_summary(fun.y=mean, +geom = "errorbar", +aes(ymax = ..y.., ymin = ..y..), +width = 0.6, +size=0.4, +linetype = "dashed", +position = position_dodge(preserve = "single")) + +facet_wrap(~pCO2_parents) +F2_Dinge_boxpots +F2_Dinge_boxpots <- RR_master_F2_Dhinge %>% # omits a single dramatic outlier +dplyr::filter(!resp_umol_hr_indiv > 0.002) %>% +ggplot(aes(x = factor(pCO2_offspring), +y = resp_umol_hr_indiv, +fill = pCO2_offspring)) + +geom_boxplot(alpha = 0.5, # color hue +width=0.6, # boxplot width +outlier.size=0, # make outliers small +position = position_dodge(preserve = "single")) + +geom_point(pch = 19, +position = position_jitterdodge(0.5), +size=1) + +scale_fill_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), +values=c("forestgreen","orange", "purple")) + +#scale_colour_manual(values=c("forestgreen","orange")) + +theme_classic() + +ggtitle("F2 Dhinge 20230407 (parentpCO2)") + +theme(axis.title.y=element_text(size=7), +axis.title.x=element_text(size=7), +axis.text.x=element_text(size=7), +legend.position = "none") + +#ylim(0, 0.2) + +stat_summary(fun.y=mean, +geom = "errorbar", +aes(ymax = ..y.., ymin = ..y..), +width = 0.6, +size=0.4, +linetype = "dashed", +position = position_dodge(preserve = "single")) + +facet_wrap(~pCO2_parents) +F2_Dinge_boxpots +F2_Dinge_boxpots <- RR_master_F2_Dhinge %>% # omits a single dramatic outlier +dplyr::filter(!resp_umol_hr_indiv > 0.0002) %>% +ggplot(aes(x = factor(pCO2_offspring), +y = resp_umol_hr_indiv, +fill = pCO2_offspring)) + +geom_boxplot(alpha = 0.5, # color hue +width=0.6, # boxplot width +outlier.size=0, # make outliers small +position = position_dodge(preserve = "single")) + +geom_point(pch = 19, +position = position_jitterdodge(0.5), +size=1) + +scale_fill_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), +values=c("forestgreen","orange", "purple")) + +#scale_colour_manual(values=c("forestgreen","orange")) + +theme_classic() + +ggtitle("F2 Dhinge 20230407 (parentpCO2)") + +theme(axis.title.y=element_text(size=7), +axis.title.x=element_text(size=7), +axis.text.x=element_text(size=7), +legend.position = "none") + +#ylim(0, 0.2) + +stat_summary(fun.y=mean, +geom = "errorbar", +aes(ymax = ..y.., ymin = ..y..), +width = 0.6, +size=0.4, +linetype = "dashed", +position = position_dodge(preserve = "single")) + +facet_wrap(~pCO2_parents) +F2_Dinge_boxpots +# Summarise Percent Deformities for plotting +F2_DF_mean_tank <- RR_master_F2_Dhinge %>% +dplyr::filter(!resp_umol_hr_indiv > 0.0002) %>% +dplyr::select(Date, pCO2_parents, pCO2_offspring,Chamber_tank, resp_umol_hr_indiv) %>% +dplyr::group_by(Date, pCO2_parents,pCO2_offspring, Chamber_tank) %>% +dplyr::summarise(mean_RR = mean(resp_umol_hr_indiv), +n = n(), +sd_RR = sd(mean_RR), +se_RR = sd_RR/(sqrt(n))) +F2_DF_mean_overall <- F2_DF_mean_tank %>% +dplyr::select(Date, pCO2_parents, pCO2_offspring, mean_RR) %>% +dplyr::group_by(Date, pCO2_parents,pCO2_offspring) %>% +dplyr::summarise(mean_RR_overall = mean(mean_RR), +n = n(), +sd_RR = sd(mean_RR), +se_RR = sd_RR/(sqrt(n))) +F2_Dinge_barplot <- ggplot(F2_DF_mean_tank) + +geom_errorbar(aes(x=pCO2_offspring, +ymin=mean_RR-se_RR, +ymax=mean_RR+se_RR), +width=0, # removes the horizontal line +colour="black", +size=1) + +geom_bar(aes(x=pCO2_offspring, y=mean_RR,fill=factor(pCO2_offspring)), +stat="identity", +width = 0.75, +alpha = 0.5) + +labs(title="F2 Dhinge RR - Mean +- SE", +x ="pCO2 Offspring Exposure", +y = "RR") + +scale_fill_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), +values=c("forestgreen","orange", "purple")) + +theme_classic() + +theme(panel.grid.major = element_blank(), +panel.grid.minor = element_blank(), +axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), +axis.text=element_text(size=12), +legend.position="none") + +facet_wrap(~pCO2_parents) # scales = "free_y") +F2_Dinge_MeanSE <- F2_DF_mean_tank %>% +ggplot(aes(x=pCO2_offspring, +y=mean_RR, +colour=pCO2_offspring)) + +# scale_linetype(c("dotted","solid")) + +scale_colour_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), +values=c("forestgreen","orange", "purple")) + +geom_point(aes(colour = pCO2_offspring), +position = position_dodge2(width = 0.75)) + +stat_summary(fun.y="mean", size = 0.8, +position = position_dodge2(width = 1)) + +stat_summary(fun.min = function(x) mean(x) - sd(x)/sqrt(length(x)), +fun.max = function(x) mean(x) + sd(x)/sqrt(length(x)), +geom = 'errorbar', +position = position_dodge2(width = 1)) + +labs(title="F2 Dhinge RR - Mean +- SE", x ="pCO2 offspring", y = "RR") + +theme_classic() + +theme(panel.grid.major = element_blank(), +panel.grid.minor = element_blank(), +axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), +axis.text=element_text(size=12), +legend.position="none") + +facet_wrap(~pCO2_parents) +F2_Dinge_MeanSE +F2_Dinge_MeanSE <- F2_DF_mean_tank %>% +ggplot(aes(x=pCO2_offspring, +y=mean_RR, +colour=pCO2_offspring)) + +# scale_linetype(c("dotted","solid")) + +scale_colour_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), +values=c("forestgreen","orange", "purple")) + +geom_point(aes(colour = pCO2_offspring), +position = position_dodge2(width = 0.75)) + +stat_summary(fun.y="mean", size = 0.8, +position = position_dodge2(width = 1)) + +stat_summary(fun.min = function(x) mean(x) - sd(x)/sqrt(length(x)), +fun.max = function(x) mean(x) + sd(x)/sqrt(length(x)), +geom = 'errorbar', +position = position_dodge2(width = 1)) + +labs(title="F2 Dhinge RR - Mean +- SE", x ="pCO2 offspring", y = "RR") + +theme_classic() + +theme(panel.grid.major = element_blank(), +panel.grid.minor = element_blank(), +axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), +axis.text=element_text(size=12), +legend.position="none") + +facet_wrap(~pCO2_parents) +F2_Dinge_boxpots <- RR_master_F2_Dhinge %>% # omits a single dramatic outlier +dplyr::filter(!resp_umol_hr_indiv > 0.0002) %>% +ggplot(aes(x = factor(pCO2_offspring), +y = resp_umol_hr_indiv, +fill = pCO2_offspring)) + +geom_boxplot(alpha = 0.5, # color hue +width=0.6, # boxplot width +outlier.size=0, # make outliers small +position = position_dodge(preserve = "single")) + +geom_point(pch = 19, +position = position_jitterdodge(0.5), +size=1) + +scale_fill_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), +values=c("forestgreen","orange", "purple")) + +#scale_colour_manual(values=c("forestgreen","orange")) + +theme_classic() + +ggtitle("F2 Dhinge 20230407 (parentpCO2)") + +theme(axis.title.y=element_text(size=7), +axis.title.x=element_text(size=7), +axis.text.x=element_text(size=7), +legend.position = "none") + +#ylim(0, 0.2) + +stat_summary(fun.y=mean, +geom = "errorbar", +aes(ymax = ..y.., ymin = ..y..), +width = 0.6, +size=0.4, +linetype = "dashed", +position = position_dodge(preserve = "single")) + +facet_wrap(~pCO2_parents) +F2_Dinge_boxpots +F2_Dinge_MeanSE +pdf(paste0(filename = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Respiration/plots/F2_RR_larvae_boxplots.pdf"), width = 8, height = 8) +print(F2_Dinge_boxpots) dev.off() -pdf("C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Survival/F3_moderate_history.pdf") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.3)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") +# OUTPUT THE PLOT +pdf(paste0(filename = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Respiration/plots/F2_RR_larvae.pdf"), width = 5, height = 5) +print(F2_Dinge_MeanSE) dev.off() -pdf("C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Survival/F3_high_history.pdf") -ggsurvplot(survfit(surv_F3_MEANS_high ~ Larvae_pCO2, binary_df_F3_MEANS_high), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.3)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") +pdf(paste0(filename = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Respiration/plots/F2_RR_larvae_boxplots.pdf"), width = 5, height = 5) +print(F2_Dinge_boxpots) dev.off() -survFit_high <- coxph(surv_F3_high ~ Larvae_pCO2 + (1|Tank), binary_df_F3_high) -survFitME_all_tank <- coxme(surv_F3_all ~ Larvae_pCO2*Parent_pCO2 + (1|Tank), binary_df_F3_all) -# remove "none" and the level - for some reason this lingered, dont know why -str(binary_df_F3_all) # before running, youll see PArent_pCO2 has 4 levels - we only want 3! -is.na(binary_df_F3_all$Parent_pCO2) <- binary_df_F3_all$Parent_pCO2 == "None" -binary_df_F3_all$Parent_pCO2 <- factor(binary_df_F3_all$Parent_pCO2) -str(binary_df_F3_all) # youll see now its fixed -survFitME_all_tank <- coxme(surv_F3_all ~ Larvae_pCO2*Parent_pCO2 + (1|Tank), binary_df_F3_all) -binary_df_F3_all -# for loop for stats! -# difference here is that tank is included - use this as a random effect in the model! -# Call the cumulative dataframe that we will write to in the for loop below -binary_df_F3_all <- data.frame() # start dataframe -binary_df_F3_high <- data.frame() # start dataframe -binary_df_F3_mod <- data.frame() # start dataframe -binary_df_F3_low <- data.frame() # start dataframe -# run it -for (i in 1:nrow(df_F3)) { -count_alive <- round(df_F3[i,]$Count_alive/1000) # divide all counts by 1000, round to nearest -count_dead <- round(df_F3[i,]$Count_dead/1000) # divide all counts by 1000, round to nearest -loopDF <- data.frame(matrix(0,ncol = 6, nrow = (count_alive + count_dead))) -colnames(loopDF) <- (c('Parent_pCO2','Larvae_pCO2','Generation','Tank','Age','Count_dead')) -loopDF <- loopDF %>% -dplyr::mutate( -Parent_pCO2 = df_F3[i,]$Parent_pCO2, -Larvae_pCO2 = df_F3[i,]$Larvae_pCO2, -Generation = df_F3[i,]$Generation, -Tank = df_F3[i,]$Tank, -Age = df_F3[i,]$Age, -) -loopDF[c(1:count_dead),6] = 1 # 1 for dead -loopDF[c((count_dead + 1):nrow(loopDF)),6] = 0 # 0 for alive -loopDF <- data.frame(loopDF) # name dataframe for this single row -binary_df_F3_all <- rbind(binary_df_F3_all,loopDF) #bind to a cumulative list dataframe -# print(binary_df_F3) # print to monitor progress -if (loopDF$Parent_pCO2[1] == 'High pCO2') { -binary_df_F3_high <- rbind(binary_df_F3_high,loopDF) #bind to a cumulative list dataframe -} else if (loopDF$Parent_pCO2[1] == 'Moderate pCO2') { -binary_df_F3_mod <- rbind(binary_df_F3_mod,loopDF) #bind to a cumulative list dataframe -} else { -binary_df_F3_low <- rbind(binary_df_F3_low,loopDF) #bind to a cumulative list dataframe -} -} -# create surv object -# with tank included in the data frame (first for loop above) -surv_F3_all <- Surv(time = as.numeric(binary_df_F3_all$Age), -event = as.numeric(binary_df_F3_all$Count_dead), type = "right") -surv_F3_high <- Surv(time = as.numeric(binary_df_F3_high$Age), -event = as.numeric(binary_df_F3_high$Count_dead), type = "right") -surv_F3_mod <- Surv(time = as.numeric(binary_df_F3_mod$Age), -event = as.numeric(binary_df_F3_mod$Count_dead), type = "right") -surv_F3_low <- Surv(time = as.numeric(binary_df_F3_low$Age), -event = as.numeric(binary_df_F3_low$Count_dead), type = "right") -# remove "none" and the level - for some reason this lingered, dont know why -str(binary_df_F3_all) # before running, youll see PArent_pCO2 has 4 levels - we only want 3! -is.na(binary_df_F3_all$Parent_pCO2) <- binary_df_F3_all$Parent_pCO2 == "None" -binary_df_F3_all$Parent_pCO2 <- factor(binary_df_F3_all$Parent_pCO2) -str(binary_df_F3_all) # youll see now its fixed -survFitME_all_tank <- coxme(surv_F3_all ~ Larvae_pCO2*Parent_pCO2 + (1|Tank), binary_df_F3_all) -summary(survFitME_all_tank) -glht(survFitME_all_tank, linfct = mcp(Parent_pCO2 = "Dunnett"),alternative = "greater") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.3)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_low ~ Larvae_pCO2, binary_df_F3_MEANS_low), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.3)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.3)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.5)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.2)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") -pdf("C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Survival/F3_low_history.pdf") -ggsurvplot(survfit(surv_F3_MEANS_low ~ Larvae_pCO2, binary_df_F3_MEANS_low), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.3)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") +F2_Dinge_boxpots <- RR_master_F2_Dhinge %>% # omits a single dramatic outlier +dplyr::filter(!resp_umol_hr_indiv > 0.0002) %>% +ggplot(aes(x = factor(pCO2_offspring), +y = resp_umol_hr_indiv, +fill = pCO2_offspring)) + +geom_boxplot(alpha = 0.5, # color hue +width=0.6, # boxplot width +outlier.size=0, # make outliers small +position = position_dodge(preserve = "single")) + +geom_point(pch = 19, +position = position_jitterdodge(0.5), +size=1) + +scale_fill_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), +values=c("forestgreen","orange", "purple")) + +#scale_colour_manual(values=c("forestgreen","orange")) + +theme_classic() + +ggtitle("F2 Dhinge") + +theme(axis.title.y=element_text(size=7), +axis.title.x=element_text(size=7), +axis.text.x=element_text(size=7), +legend.position = "none") + +#ylim(0, 0.2) + +stat_summary(fun.y=mean, +geom = "errorbar", +aes(ymax = ..y.., ymin = ..y..), +width = 0.6, +size=0.4, +linetype = "dashed", +position = position_dodge(preserve = "single")) + +facet_wrap(~pCO2_parents) +pdf(paste0(filename = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Respiration/plots/F2_RR_larvae_boxplots.pdf"), width = 5, height = 5) +print(F2_Dinge_boxpots) dev.off() -pdf("C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Survival/F3_moderate_history.pdf") -ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.2)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") +F2_Dinge_MeanSE <- F2_DF_mean_tank %>% +ggplot(aes(x=pCO2_offspring, +y=mean_RR, +colour=pCO2_offspring)) + +# scale_linetype(c("dotted","solid")) + +scale_colour_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), +values=c("forestgreen","orange", "purple")) + +geom_point(aes(colour = pCO2_offspring), +position = position_dodge2(width = 0.75)) + +stat_summary(fun.y="mean", size = 0.8, +position = position_dodge2(width = 1)) + +stat_summary(fun.min = function(x) mean(x) - sd(x)/sqrt(length(x)), +fun.max = function(x) mean(x) + sd(x)/sqrt(length(x)), +geom = 'errorbar', +position = position_dodge2(width = 1)) + +labs(title="F2 Dhinge RR - Mean +- SE", x ="pCO2 offspring", y = "RR") + +theme_classic() + +theme(axis.title.y=element_text(size=7), +axis.title.x=element_text(size=7), +axis.text.x=element_text(size=7), +legend.position = "none") + +facet_wrap(~pCO2_parents) +F2_Dinge_MeanSE +pdf(paste0(filename = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Respiration/plots/F2_RR_larvae.pdf"), width = 5, height = 5) +print(F2_Dinge_MeanSE) dev.off() -pdf("C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Survival/F3_high_history.pdf") -ggsurvplot(survfit(surv_F3_MEANS_high ~ Larvae_pCO2, binary_df_F3_MEANS_high), -risk.table = F, pval = F , conf.int = TRUE, -palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.3)), -break.x.by = 2 ,xlab = "Age (days post fertilization)", -legend.title = "treatment") +pdf(paste0(filename = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Respiration/plots/F2_RR_larvae_boxplots.pdf"), width = 5, height = 5) +print(F2_Dinge_boxpots) dev.off() -survFitME_all_tank <- coxme(surv_F3_all ~ Larvae_pCO2*Parent_pCO2 + (1|Tank), binary_df_F3_all) -summary(survFitME_all_tank) -survFitME_all_tank <- coxme(surv_F3_all ~ Larvae_pCO2*Parent_pCO2, binary_df_F3_all) -survFitME_all_tank <- coxme(surv_F3_all ~ Larvae_pCO2*Parent_pCO2 + (1|Tank), binary_df_F3_all) -summary(survFitME_all_tank) -knitr::opts_chunk$set(echo = TRUE, -warning = FALSE, -message = FALSE, -cache = TRUE) -knitr::opts_knit$set(root.dir = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB_EPA_LISS_Disease_Surveillance/Sonde_Data") -# load libraries - notes show the install command needed to install (pre installed) -# Plotting -library(ggplot2) +print(F2_Dinge_boxpots) +# OUTPUT THE PLOT +pdf(paste0(filename = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Respiration/plots/F2_RR_larvae.pdf"), width = 6, height = 4) +print(F2_Dinge_MeanSE) +dev.off() +pdf(paste0(filename = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Respiration/plots/F2_RR_larvae_boxplots.pdf"), width = 6, height = 4) +print(F2_Dinge_boxpots) +dev.off() +RR_master_F2_Dhinge_OM <- RR_master_F2_Dhinge %>% dplyr::filter(!resp_umol_hr_indiv > 0.0002) %>% +``` +RR_master_F2_Dhinge_OM <- RR_master_F2_Dhinge %>% dplyr::filter(!resp_umol_hr_indiv > 0.0002) +RR_master_F2_Dhinge_OM +F2_DF_mean_tank <- RR_master_F2_Dhinge %>% +dplyr::filter(!resp_umol_hr_indiv > 0.0002) %>% +dplyr::select(Date, pCO2_parents, pCO2_offspring,Chamber_tank, resp_umol_hr_indiv) %>% +dplyr::group_by(Date, pCO2_parents,pCO2_offspring, Chamber_tank) %>% +dplyr::summarise(mean_RR = mean(resp_umol_hr_indiv), +n = n(), +sd_RR = sd(mean_RR), +se_RR = sd_RR/(sqrt(n))) +F2_TwoWayANOVA <- lm(mean_RR~pCO2_parents*pCO2_offspring, data=F2_DF_mean_tank) +shapiro.test(resid(F2_TwoWayANOVA)) +leveneTest(F2_TwoWayANOVA) +summary(aov(F2_TwoWayANOVA)) +F3_DF_mean_tank <- F3_DF %>% +dplyr::select(Date, pCO2_parents, pCO2_offspring, Chamber_tank, resp_umol_hr_indiv) %>% +na.omit() %>% +dplyr::group_by(Date, pCO2_parents,pCO2_offspring, Chamber_tank) %>% +dplyr::summarise(mean_RR = mean(resp_umol_hr_indiv), +n = n(), +sd_RR = sd(resp_umol_hr_indiv), +se_RR = sd_RR/(sqrt(n))) +F3_DF_mean_tank +F3_DF_dhinge <- F3_DF_mean_tank %>% dplyr::filter(Date %in%"4/7/2023") +F3_DF_day7 <- F3_DF_mean_tank %>% dplyr::filter(Date %in%"4/12/2023") +F3_DF_day16 <- F3_DF_mean_tank %>% dplyr::filter(Date %in%"4/21/2023") +F3_TwoWayANOVA_dhinge <- lm(mean_RR~pCO2_parents*pCO2_offspring, data=F3_DF_dhinge) +shapiro.test(resid(F3_TwoWayANOVA_dhinge)) # 0.3896 +leveneTest(F3_TwoWayANOVA_dhinge) # 0.8975 +summary(aov(F3_TwoWayANOVA_dhinge)) +F3_TwoWayANOVA_day7 <- lm(mean_RR~pCO2_parents*pCO2_offspring, data=F3_DF_day7) +shapiro.test(resid(F3_TwoWayANOVA_day7)) # 0.1789 +leveneTest(F3_TwoWayANOVA_day7) # 0.7962 +F3_SRH_day7 <- scheirerRayHare(mean_RR~pCO2_parents*pCO2_offspring, data=F3_DF_day7) +summary(F3_SRH_day7) +F3_SRH_day7 +F3_TwoWayANOVA_day16 <- lm(mean_RR~pCO2_parents*pCO2_offspring, data=F3_DF_day16) +shapiro.test(resid(F3_TwoWayANOVA_day16)) # 0.1789 +leveneTest(F3_TwoWayANOVA_day16) # 0.7962 +summary(aov(F3_TwoWayANOVA_day16)) +F3_TwoWayANOVA_dhinge <- lm(mean_RR~pCO2_parents*pCO2_offspring, data=F3_DF_dhinge) +shapiro.test(resid(F3_TwoWayANOVA_dhinge)) # 0.1789 +leveneTest(F3_TwoWayANOVA_dhinge) # 0.7962 +summary(aov(F3_TwoWayANOVA_dhinge)) +F3_Dinge_boxpots <- F3_DF %>% dplyr::filter(Date %in%"4/7/2023") %>% +ggplot(aes(x = factor(pCO2_offspring), +y = resp_umol_hr_indiv, +fill = pCO2_offspring)) + +geom_boxplot(alpha = 0.5, # color hue +width=0.6, # boxplot width +outlier.size=0, # make outliers small +position = position_dodge(preserve = "single")) + +geom_point(pch = 19, +position = position_jitterdodge(0.5), +size=1) + +scale_fill_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), +values=c("forestgreen","orange", "purple")) + +#scale_colour_manual(values=c("forestgreen","orange")) + +theme_classic() + +ggtitle("F2 Dhinge") + +theme(axis.title.y=element_text(size=7), +axis.title.x=element_text(size=7), +axis.text.x=element_text(size=7), +legend.position = "none") + +#ylim(0, 0.2) + +stat_summary(fun.y=mean, +geom = "errorbar", +aes(ymax = ..y.., ymin = ..y..), +width = 0.6, +size=0.4, +linetype = "dashed", +position = position_dodge(preserve = "single")) + +facet_wrap(~pCO2_parents) +F3_Dinge_boxpots +F3_Dinge_boxpots <- F3_DF %>% dplyr::filter(Date %in%"4/7/2023") %>% +ggplot(aes(x = factor(pCO2_offspring), +y = resp_umol_hr_indiv, +fill = pCO2_offspring)) + +geom_boxplot(alpha = 0.5, # color hue +width=0.6, # boxplot width +outlier.size=0, # make outliers small +position = position_dodge(preserve = "single")) + +geom_point(pch = 19, +position = position_jitterdodge(0.5), +size=1) + +scale_fill_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), +values=c("forestgreen","orange", "purple")) + +#scale_colour_manual(values=c("forestgreen","orange")) + +theme_classic() + +ggtitle("F3 Dhinge") + +theme(axis.title.y=element_text(size=7), +axis.title.x=element_text(size=7), +axis.text.x=element_text(size=7), +legend.position = "none") + +#ylim(0, 0.2) + +stat_summary(fun.y=mean, +geom = "errorbar", +aes(ymax = ..y.., ymin = ..y..), +width = 0.6, +size=0.4, +linetype = "dashed", +position = position_dodge(preserve = "single")) + +facet_wrap(~pCO2_parents) +F3_Dinge_boxpots +F3_day7_boxpots <- F3_DF %>% dplyr::filter(Date %in%"4/12/2023") %>% +ggplot(aes(x = factor(pCO2_offspring), +y = resp_umol_hr_indiv, +fill = pCO2_offspring)) + +geom_boxplot(alpha = 0.5, # color hue +width=0.6, # boxplot width +outlier.size=0, # make outliers small +position = position_dodge(preserve = "single")) + +geom_point(pch = 19, +position = position_jitterdodge(0.5), +size=1) + +scale_fill_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), +values=c("forestgreen","orange", "purple")) + +#scale_colour_manual(values=c("forestgreen","orange")) + +theme_classic() + +ggtitle("F3 day 7") + +theme(axis.title.y=element_text(size=7), +axis.title.x=element_text(size=7), +axis.text.x=element_text(size=7), +legend.position = "none") + +#ylim(0, 0.2) + +stat_summary(fun.y=mean, +geom = "errorbar", +aes(ymax = ..y.., ymin = ..y..), +width = 0.6, +size=0.4, +linetype = "dashed", +position = position_dodge(preserve = "single")) + +facet_wrap(~pCO2_parents) +F3_day16_boxpots <- F3_DF %>% dplyr::filter(Date %in%"4/21/2023") %>% +ggplot(aes(x = factor(pCO2_offspring), +y = resp_umol_hr_indiv, +fill = pCO2_offspring)) + +geom_boxplot(alpha = 0.5, # color hue +width=0.6, # boxplot width +outlier.size=0, # make outliers small +position = position_dodge(preserve = "single")) + +geom_point(pch = 19, +position = position_jitterdodge(0.5), +size=1) + +scale_fill_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), +values=c("forestgreen","orange", "purple")) + +#scale_colour_manual(values=c("forestgreen","orange")) + +theme_classic() + +ggtitle("F3 day 16") + +theme(axis.title.y=element_text(size=7), +axis.title.x=element_text(size=7), +axis.text.x=element_text(size=7), +legend.position = "none") + +#ylim(0, 0.2) + +stat_summary(fun.y=mean, +geom = "errorbar", +aes(ymax = ..y.., ymin = ..y..), +width = 0.6, +size=0.4, +linetype = "dashed", +position = position_dodge(preserve = "single")) + +facet_wrap(~pCO2_parents) +print(ggarrange(F3_Dinge_boxpots,F3_day7_boxpots,F3_day16_boxpots, nrow =3)) +pdf(paste0("C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Respiration/plots/F3_RR_larvae_boxplots.pdf"), width=8, height=12) +print(ggarrange(F3_Dinge_boxpots,F3_day7_boxpots,F3_day16_boxpots, nrow =3)) +dev.off() +# LOAD PACKAGES ::::::::::::::::::::::::::::::::::::::::::::::::::::::: library(dplyr) -library(gplots) -library(RColorBrewer) -library(ggpubr) +library(ggplot2) +library(forcats) +library(lme4) +library(lmerTest) +library(see) +library(performance) +library(car) +library(kableExtra) +library(pander) +library(data.table) library(stringr) -library(lubridate) -library(readr) -path_out = 'C:/Users/samjg/Documents/Github_repositories/EAD-ASEB_EPA_LISS_Disease_Surveillance/Sonde_Data/output' -path.p <- "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB_EPA_LISS_Disease_Surveillance/Sonde_Data" -file.names.table <- data.frame(txt.files = -(basename(list.files(path = paste(path.p,'/',sep=''), -pattern = "csv$", -recursive = FALSE)))) -Master_Table <- data.frame() -# Notes on formatting datas -# ymd_hms for 072023_GOLD_Sonde.csv, 072023_LAUR_Sonde.csv, 082023_FENC_Sonde.csv, 082023_GOLD_Sonde.csv, 092023_FENC_Sonde.csv -# mdy_hms for 072023_ASHC_Sonde.csv -# mdy_hm for 082023_ASHC_Sonde.csv, 082023_LAUR_Sonde.csv, 092023_ASHC_Sonde.csv, 092023_GOLD_Sonde.csv, 102023_ASHC_Sonde.csv,102023_FENC_Sonde.csv, 102023_GOLD_Sonde.csv, 102023_LAUR_Sonde.csv, 112023_ASHC_Sonde.csv, 112023_FENC_Sonde.csv -# for loop to merge all csv files to one dataframe -for (m in 1:nrow(file.names.table)) { -# call the dir of the csv file in loop -raw_Sonde_rootdir <- paste(path.p,'/',file.names.table[m,1], sep='') #reads in the data files -# read it and fix dates and whatnot yay -D <- readLines(raw_Sonde_rootdir) -ind <- grep('Date Time',D) # the line where the data stars (theres a lot of junk beforehand) -# raw <- read.csv(paste(path.p,'/',file.names.table[m,1], sep=''),skip=ind-1) -raw <- read_csv(raw_Sonde_rootdir,skip=ind-1,col_types = cols()) -columns <- names(raw) # what are the column names? - youll see why -raw_df <- as.data.frame(raw) # ocnvert the dataframe now that you skipped the crappola -raw_df[2:(ncol(raw_df))] <- lapply(raw_df[2:(ncol(raw_df))],as.numeric) # convert all data to numeric -# ommit all numeric data, also the format of the units does not pair between datasets, cannot rbind later -names(raw_df) <- sapply(strsplit(names(raw_df), '\\s*[()]'), `[`, 1) -raw_df[,1] -# now for the dates - for some reason the sondes have diff format - fix it up! -if(gsub(".*/","", (gsub("\\","/",raw_Sonde_rootdir, fixed=T))) %in% -c('082023_ASHC_Sonde.csv', -'082023_LAUR_Sonde.csv', -'092023_ASHC_Sonde.csv', -'092023_GOLD_Sonde.csv', -'102023_ASHC_Sonde.csv', -'102023_FENC_Sonde.csv', -'102023_GOLD_Sonde.csv', -'102023_LAUR_Sonde.csv', -'112023_ASHC_Sonde.csv', -'112023_FENC_Sonde.csv')) { # date formatted as mdy_hm -raw_df[,1] <- mdy_hm(raw_df[,1]) # reformat as date using lubridate -} else if (gsub(".*/","", (gsub("\\","/",raw_Sonde_rootdir, fixed=T))) %in% -'072023_ASHC_Sonde.csv') { -raw_df[,1] <- mdy_hms(raw_df[,1]) # reformat as date using lubridate -} else (raw_df[,1] <- ymd_hms(raw_df[,1]) # reformat as date using lubridate -) -# add a column for the site -filename <- file.names.table[m,1] -raw_df <- raw_df %>% -dplyr::mutate(Site = sub(".*_(.*)_.*","\\1",filename), -Date = gsub("_.*","",filename)) -Master_Table <- rbind(raw_df,Master_Table) #bind to a cumulative list dataframe -} -path_out = 'C:/Users/samjg/Documents/Github_repositories/EAD-ASEB_EPA_LISS_Disease_Surveillance/Sonde_Data/output' -path.p <- "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB_EPA_LISS_Disease_Surveillance/Sonde_Data" -file.names.table <- data.frame(txt.files = -(basename(list.files(path = paste(path.p,'/',sep=''), -pattern = "csv$", -recursive = FALSE)))) -Master_Table <- data.frame() -for (m in 1:nrow(file.names.table)) { -# call the dir of the csv file in loop -raw_Sonde_rootdir <- paste(path.p,'/',file.names.table[m,1], sep='') #reads in the data files -# read it and fix dates and whatnot yay -D <- readLines(raw_Sonde_rootdir) -ind <- grep('Date Time',D) # the line where the data stars (theres a lot of junk beforehand) -# raw <- read.csv(paste(path.p,'/',file.names.table[m,1], sep=''),skip=ind-1) -raw <- read_csv(raw_Sonde_rootdir,skip=ind-1,col_types = cols()) -columns <- names(raw) # what are the column names? - youll see why -raw_df <- as.data.frame(raw) # ocnvert the dataframe now that you skipped the crappola -raw_df[2:(ncol(raw_df))] <- lapply(raw_df[2:(ncol(raw_df))],as.numeric) # convert all data to numeric -# ommit all numeric data, also the format of the units does not pair between datasets, cannot rbind later -names(raw_df) <- sapply(strsplit(names(raw_df), '\\s*[()]'), `[`, 1) -raw_df[,1] -# now for the dates - for some reason the sondes have diff format - fix it up! -if(gsub(".*/","", (gsub("\\","/",raw_Sonde_rootdir, fixed=T))) %in% -c('082023_ASHC_Sonde.csv', -'082023_LAUR_Sonde.csv', -'092023_ASHC_Sonde.csv', -'092023_GOLD_Sonde.csv', -'102023_ASHC_Sonde.csv', -'102023_FENC_Sonde.csv', -'102023_GOLD_Sonde.csv', -'102023_LAUR_Sonde.csv', -'112023_ASHC_Sonde.csv', -'112023_FENC_Sonde.csv')) { # date formatted as mdy_hm -raw_df[,1] <- mdy_hm(raw_df[,1]) # reformat as date using lubridate -} else if (gsub(".*/","", (gsub("\\","/",raw_Sonde_rootdir, fixed=T))) %in% -'072023_ASHC_Sonde.csv') { -raw_df[,1] <- mdy_hms(raw_df[,1]) # reformat as date using lubridate -} else (raw_df[,1] <- ymd_hms(raw_df[,1]) # reformat as date using lubridate -) -# add a column for the site -filename <- file.names.table[m,1] -raw_df <- raw_df %>% -dplyr::mutate(Site = sub(".*_(.*)_.*","\\1",filename), -Date = gsub("_.*","",filename)) -Master_Table <- rbind(raw_df,Master_Table) #bind to a cumulative list dataframe -} -file.names.table -path_out = 'C:/Users/samjg/Documents/Github_repositories/EAD-ASEB_EPA_LISS_Disease_Surveillance/Sonde_Data/output' -path.p <- "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB_EPA_LISS_Disease_Surveillance/Sonde_Data" -file.names.table <- data.frame(txt.files = -(basename(list.files(path = paste(path.p,'/',sep=''), -pattern = "csv$", -recursive = FALSE)))) -file.names.table -path_out -path.p <- "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB_EPA_LISS_Disease_Surveillance/Sonde_Data" -file.names.table <- data.frame(txt.files = -(basename(list.files(path = paste(path.p,'/',sep=''), -pattern = "csv$", -recursive = FALSE)))) -file.names.table -file.names.table <- data.frame(txt.files = -(basename(list.files(path = paste(path.p,'/',sep=''), -pattern = "csv$", -recursive = FALSE)))) -file.names.table -path.p -path.p <- "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB_EPA_LISS_Disease_Surveillance/Sonde_Data" -file.names.table <- data.frame(txt.files = -(basename(list.files(path = paste(path.p,'/',sep=''), -pattern = "csv$", -recursive = FALSE)))) -file.names.table -knitr::opts_chunk$set(echo = TRUE, -warning = FALSE, -message = FALSE, -cache = TRUE) -knitr::opts_knit$set(root.dir = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB_EPA_LISS_Disease_Surveillance/Sonde_Data") -# load libraries - notes show the install command needed to install (pre installed) -# Plotting +library(latex2exp) +library(Rmisc) +library(devtools) +library(ggpubr) +library(hrbrthemes) +# SET WORKING DIRECTORY ::::::::::::::::::::::::::::::::::::::::::::::: +knitr::opts_knit$set(root.dir = 'C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis') +# getwd() +resp.ref <- read.csv(file="Data/Respiration/metadata/Reference_resp_ID.csv", header=T) +resp.ref +resp.ref +summarySE(resp.ref, measurevar="Number", groupvars=c("Date", "Chamber_tank")) +knitr::opts_chunk$set(echo = TRUE) +# SET WORKING DIRECTORY +# knitr::opts_knit$set(root.dir = "C:/Users/katherine.mcfarland/Documents/GitHub/EAD-ASEB-Airradians_multigen_OA/larvae") # Katie's +knitr::opts_knit$set(root.dir = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis") # Sam's +?Surv library(ggplot2) +library(tidyr) library(dplyr) -library(gplots) -library(RColorBrewer) +library(rcompanion) +library(FSA) +library(car) +library(forcats) +library(kableExtra) # nice Rmd tables +library(emmeans) library(ggpubr) -library(stringr) -library(lubridate) -library(readr) -path_out = 'C:/Users/samjg/Documents/Github_repositories/EAD-ASEB_EPA_LISS_Disease_Surveillance/Sonde_Data/output' -path.p <- "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB_EPA_LISS_Disease_Surveillance/Sonde_Data" -file.names.table <- data.frame(txt.files = -(basename(list.files(path = paste(path.p,'/',sep=''), -pattern = "csv$", -recursive = FALSE)))) -Master_Table <- data.frame() -# Notes on formatting datas -# ymd_hms for 072023_GOLD_Sonde.csv, 072023_LAUR_Sonde.csv, 082023_FENC_Sonde.csv, 082023_GOLD_Sonde.csv, 092023_FENC_Sonde.csv -# mdy_hms for 072023_ASHC_Sonde.csv -# mdy_hm for 082023_ASHC_Sonde.csv, 082023_LAUR_Sonde.csv, 092023_ASHC_Sonde.csv, 092023_GOLD_Sonde.csv, 102023_ASHC_Sonde.csv,102023_FENC_Sonde.csv, 102023_GOLD_Sonde.csv, 102023_LAUR_Sonde.csv, 112023_ASHC_Sonde.csv, 112023_FENC_Sonde.csv -# for loop to merge all csv files to one dataframe -for (m in 1:nrow(file.names.table)) { -# call the dir of the csv file in loop -raw_Sonde_rootdir <- paste(path.p,'/',file.names.table[m,1], sep='') #reads in the data files -# read it and fix dates and whatnot yay -D <- readLines(raw_Sonde_rootdir) -ind <- grep('Date Time',D) # the line where the data stars (theres a lot of junk beforehand) -# raw <- read.csv(paste(path.p,'/',file.names.table[m,1], sep=''),skip=ind-1) -raw <- read_csv(raw_Sonde_rootdir,skip=ind-1,col_types = cols()) -columns <- names(raw) # what are the column names? - youll see why -raw_df <- as.data.frame(raw) # ocnvert the dataframe now that you skipped the crappola -raw_df[2:(ncol(raw_df))] <- lapply(raw_df[2:(ncol(raw_df))],as.numeric) # convert all data to numeric -# ommit all numeric data, also the format of the units does not pair between datasets, cannot rbind later -names(raw_df) <- sapply(strsplit(names(raw_df), '\\s*[()]'), `[`, 1) -raw_df[,1] -# now for the dates - for some reason the sondes have diff format - fix it up! -if(gsub(".*/","", (gsub("\\","/",raw_Sonde_rootdir, fixed=T))) %in% -c('082023_ASHC_Sonde.csv', -'082023_LAUR_Sonde.csv', -'092023_ASHC_Sonde.csv', -'092023_GOLD_Sonde.csv', -'102023_ASHC_Sonde.csv', -'102023_FENC_Sonde.csv', -'102023_GOLD_Sonde.csv', -'102023_LAUR_Sonde.csv', -'112023_ASHC_Sonde.csv', -'112023_FENC_Sonde.csv')) { # date formatted as mdy_hm -raw_df[,1] <- mdy_hm(raw_df[,1]) # reformat as date using lubridate -} else if (gsub(".*/","", (gsub("\\","/",raw_Sonde_rootdir, fixed=T))) %in% -'072023_ASHC_Sonde.csv') { -raw_df[,1] <- mdy_hms(raw_df[,1]) # reformat as date using lubridate -} else (raw_df[,1] <- ymd_hms(raw_df[,1]) # reformat as date using lubridate -) -# add a column for the site -filename <- file.names.table[m,1] -raw_df <- raw_df %>% -dplyr::mutate(Site = sub(".*_(.*)_.*","\\1",filename), -Date = gsub("_.*","",filename)) -Master_Table <- rbind(raw_df,Master_Table) #bind to a cumulative list dataframe -} -# call the dir of the csv file in loop -raw_Sonde_rootdir <- paste(path.p,'/',file.names.table[m,1], sep='') #reads in the data files -raw_Sonde_rootdir -file.names.table -file.names.table -(basename(list.files(path = paste(path.p,'/',sep='') -paste(path.p,'/',sep='') +library(survival) +library(Rmisc) +library(coxme) +library(survminer) +?Surv +?coxme +?glht diff --git a/RAnalysis/Scripts/RespirationRates_analysis_larvae.Rmd b/RAnalysis/Scripts/RespirationRates_analysis_larvae.Rmd index c062406..a42a1df 100644 --- a/RAnalysis/Scripts/RespirationRates_analysis_larvae.Rmd +++ b/RAnalysis/Scripts/RespirationRates_analysis_larvae.Rmd @@ -206,16 +206,41 @@ F2_Dinge_boxpots <- RR_master_F2_Dhinge %>% # omits a single dramatic outlier facet_wrap(~pCO2_parents) # OUTPUT THE PLOT -pdf(paste0(filename = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Respiration/plots/F2_RR_larvae.pdf"), width = 5, height = 5) +pdf(paste0(filename = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Respiration/plots/F2_RR_larvae.pdf"), width = 6, height = 4) print(F2_Dinge_MeanSE) dev.off() -pdf(paste0(filename = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Respiration/plots/F2_RR_larvae_boxplots.pdf"), width = 5, height = 5) +pdf(paste0(filename = "C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Respiration/plots/F2_RR_larvae_boxplots.pdf"), width = 6, height = 4) print(F2_Dinge_boxpots) dev.off() ``` +## F2 stats + +```{r F2 stats} + +F2_DF_mean_tank <- RR_master_F2_Dhinge %>% + dplyr::filter(!resp_umol_hr_indiv > 0.0002) %>% + dplyr::select(Date, pCO2_parents, pCO2_offspring,Chamber_tank, resp_umol_hr_indiv) %>% + dplyr::group_by(Date, pCO2_parents,pCO2_offspring, Chamber_tank) %>% + dplyr::summarise(mean_RR = mean(resp_umol_hr_indiv), + n = n(), + sd_RR = sd(mean_RR), + se_RR = sd_RR/(sqrt(n))) + +F2_TwoWayANOVA <- lm(mean_RR~pCO2_parents*pCO2_offspring, data=F2_DF_mean_tank) +shapiro.test(resid(F2_TwoWayANOVA)) # 0.3896 +leveneTest(F2_TwoWayANOVA) # 0.8975 +summary(aov(F2_TwoWayANOVA)) +# Df Sum Sq Mean Sq F value Pr(>F) +# pCO2_parents 1 2.057e-09 2.057e-09 3.100 0.112 +# pCO2_offspring 2 5.980e-10 2.990e-10 0.451 0.651 +# pCO2_parents:pCO2_offspring 2 6.520e-10 3.259e-10 0.491 0.627 +# Residuals 9 5.973e-09 6.637e-10 +``` + + # F3 figures all merged! @@ -495,6 +520,156 @@ print(ggarrange(Gen3_MEANSE_Dhinge,Gen3_MEANSE_day7,Gen3_MEANSE_day16, nrow =3)) dev.off() + +F3_Dinge_boxpots <- F3_DF %>% dplyr::filter(Date %in%"4/7/2023") %>% + ggplot(aes(x = factor(pCO2_offspring), + y = resp_umol_hr_indiv, + fill = pCO2_offspring)) + + geom_boxplot(alpha = 0.5, # color hue + width=0.6, # boxplot width + outlier.size=0, # make outliers small + position = position_dodge(preserve = "single")) + + geom_point(pch = 19, + position = position_jitterdodge(0.5), + size=1) + + scale_fill_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), + values=c("forestgreen","orange", "purple")) + + #scale_colour_manual(values=c("forestgreen","orange")) + + theme_classic() + + ggtitle("F3 Dhinge") + + theme(axis.title.y=element_text(size=7), + axis.title.x=element_text(size=7), + axis.text.x=element_text(size=7), + legend.position = "none") + + #ylim(0, 0.2) + + stat_summary(fun.y=mean, + geom = "errorbar", + aes(ymax = ..y.., ymin = ..y..), + width = 0.6, + size=0.4, + linetype = "dashed", + position = position_dodge(preserve = "single")) + + facet_wrap(~pCO2_parents) + + +F3_day7_boxpots <- F3_DF %>% dplyr::filter(Date %in%"4/12/2023") %>% + ggplot(aes(x = factor(pCO2_offspring), + y = resp_umol_hr_indiv, + fill = pCO2_offspring)) + + geom_boxplot(alpha = 0.5, # color hue + width=0.6, # boxplot width + outlier.size=0, # make outliers small + position = position_dodge(preserve = "single")) + + geom_point(pch = 19, + position = position_jitterdodge(0.5), + size=1) + + scale_fill_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), + values=c("forestgreen","orange", "purple")) + + #scale_colour_manual(values=c("forestgreen","orange")) + + theme_classic() + + ggtitle("F3 day 7") + + theme(axis.title.y=element_text(size=7), + axis.title.x=element_text(size=7), + axis.text.x=element_text(size=7), + legend.position = "none") + + #ylim(0, 0.2) + + stat_summary(fun.y=mean, + geom = "errorbar", + aes(ymax = ..y.., ymin = ..y..), + width = 0.6, + size=0.4, + linetype = "dashed", + position = position_dodge(preserve = "single")) + + facet_wrap(~pCO2_parents) + +F3_day16_boxpots <- F3_DF %>% dplyr::filter(Date %in%"4/21/2023") %>% + ggplot(aes(x = factor(pCO2_offspring), + y = resp_umol_hr_indiv, + fill = pCO2_offspring)) + + geom_boxplot(alpha = 0.5, # color hue + width=0.6, # boxplot width + outlier.size=0, # make outliers small + position = position_dodge(preserve = "single")) + + geom_point(pch = 19, + position = position_jitterdodge(0.5), + size=1) + + scale_fill_manual(breaks=c("500 uatm", "800 uatm", "1200 uatm"), + values=c("forestgreen","orange", "purple")) + + #scale_colour_manual(values=c("forestgreen","orange")) + + theme_classic() + + ggtitle("F3 day 16") + + theme(axis.title.y=element_text(size=7), + axis.title.x=element_text(size=7), + axis.text.x=element_text(size=7), + legend.position = "none") + + #ylim(0, 0.2) + + stat_summary(fun.y=mean, + geom = "errorbar", + aes(ymax = ..y.., ymin = ..y..), + width = 0.6, + size=0.4, + linetype = "dashed", + position = position_dodge(preserve = "single")) + + facet_wrap(~pCO2_parents) + +pdf(paste0("C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Respiration/plots/F3_RR_larvae_boxplots.pdf"), width=8, height=12) +print(ggarrange(F3_Dinge_boxpots,F3_day7_boxpots,F3_day16_boxpots, nrow =3)) +dev.off() + +``` + + + +## F3 stats + +```{r F3 stats} + +F3_DF_mean_tank <- F3_DF %>% + dplyr::select(Date, pCO2_parents, pCO2_offspring, Chamber_tank, resp_umol_hr_indiv) %>% + na.omit() %>% + dplyr::group_by(Date, pCO2_parents,pCO2_offspring, Chamber_tank) %>% + dplyr::summarise(mean_RR = mean(resp_umol_hr_indiv), + n = n(), + sd_RR = sd(resp_umol_hr_indiv), + se_RR = sd_RR/(sqrt(n))) + + +F3_DF_dhinge <- F3_DF_mean_tank %>% dplyr::filter(Date %in%"4/7/2023") +F3_DF_day7 <- F3_DF_mean_tank %>% dplyr::filter(Date %in%"4/12/2023") +F3_DF_day16 <- F3_DF_mean_tank %>% dplyr::filter(Date %in%"4/21/2023") + +F3_TwoWayANOVA_dhinge <- lm(mean_RR~pCO2_parents*pCO2_offspring, data=F3_DF_dhinge) +shapiro.test(resid(F3_TwoWayANOVA_dhinge)) # 0.1789 +leveneTest(F3_TwoWayANOVA_dhinge) # 0.7962 +summary(aov(F3_TwoWayANOVA_dhinge)) +# Df Sum Sq Mean Sq F value Pr(>F) +# pCO2_parents 2 3.886e-08 1.943e-08 2.137 0.147 +# pCO2_offspring 2 1.100e-09 5.520e-10 0.061 0.941 +# pCO2_parents:pCO2_offspring 4 7.604e-08 1.901e-08 2.091 0.124 +# Residuals 18 1.636e-07 9.089e-09 + +F3_TwoWayANOVA_day7 <- lm(mean_RR~pCO2_parents*pCO2_offspring, data=F3_DF_day7) +shapiro.test(resid(F3_TwoWayANOVA_day7)) # 6.619e-05 +leveneTest(F3_TwoWayANOVA_day7) # 0.319 +F3_SRH_day7 <- scheirerRayHare(mean_RR~pCO2_parents*pCO2_offspring, data=F3_DF_day7) +# Df Sum Sq H p.value +# pCO2_parents 2 72.67 1.15344 0.56174 +# pCO2_offspring 2 44.67 0.70899 0.70153 +# pCO2_parents:pCO2_offspring 4 69.33 1.10053 0.89419 +# Residuals 18 1451.33 + + + + +F3_TwoWayANOVA_day16 <- lm(mean_RR~pCO2_parents*pCO2_offspring, data=F3_DF_day16) +shapiro.test(resid(F3_TwoWayANOVA_day16)) # 0.6191 +leveneTest(F3_TwoWayANOVA_day16) # 0.8498 +summary(aov(F3_TwoWayANOVA_day16)) +# Df Sum Sq Mean Sq F value Pr(>F) +# pCO2_parents 2 6.314e-08 3.157e-08 1.241 0.334 +# pCO2_offspring 2 3.560e-09 1.780e-09 0.070 0.933 +# pCO2_parents:pCO2_offspring 3 9.866e-08 3.289e-08 1.293 0.335 +# Residuals 9 2.289e-07 2.544e-08 ``` diff --git a/RAnalysis/Scripts/RespirationRates_calc.Rmd b/RAnalysis/Scripts/RespirationRates_calc.Rmd index 2c35dca..4b483f0 100644 --- a/RAnalysis/Scripts/RespirationRates_calc.Rmd +++ b/RAnalysis/Scripts/RespirationRates_calc.Rmd @@ -54,6 +54,8 @@ nrow(resp.data) # 720 resp.ref <- read.csv(file="Data/Respiration/metadata/Reference_resp_ID.csv", header=T) nrow(resp.ref) # 720 + + ``` diff --git a/RAnalysis/Scripts/Survival.Rmd b/RAnalysis/Scripts/Survival.Rmd index 7e431c9..b9a2420 100644 --- a/RAnalysis/Scripts/Survival.Rmd +++ b/RAnalysis/Scripts/Survival.Rmd @@ -37,7 +37,7 @@ library(survminer) #Larval Survival ```{r, echo=FALSE} -df<-read.csv("Data/larval_survival.csv", header = T) %>% +df <- read.csv("Data/larval_survival.csv", header = T) %>% dplyr::mutate(Percent_Survival = as.numeric(Percent_Survival)) %>% dplyr::select(-treatment) %>% # redundant column with Exposure_OA @@ -51,8 +51,10 @@ df<-read.csv("Data/larval_survival.csv", header = T) %>% Exposure_OA == "HIGH" ~ "High pCO2")) %>% # change the levels for order when plotting - call as factors too! - dplyr::mutate(Parental_OA = as.factor(forcats::fct_relevel(Parental_OA, c("Low pCO2", "Moderate pCO2", "High pCO2")))) %>% - dplyr::mutate(Exposure_OA = as.factor(forcats::fct_relevel(Exposure_OA, c("Low pCO2", "Moderate pCO2", "High pCO2")))) %>% + dplyr::mutate(Parental_OA = as.factor(forcats::fct_relevel(Parental_OA, + c("Low pCO2", "Moderate pCO2", "High pCO2")))) %>% + dplyr::mutate(Exposure_OA = as.factor(forcats::fct_relevel(Exposure_OA, + c("Low pCO2", "Moderate pCO2", "High pCO2")))) %>% # rename them dplyr::rename(Parent_pCO2 = Parental_OA) %>% @@ -169,8 +171,11 @@ binary_df_F1_mod <- data.frame() # start dataframe - only mod history binary_df_F1_low <- data.frame() # start dataframe - only low history # run it for (i in 1:nrow(df_F1)) { - count_alive <- round(df_F1[i,]$Count_alive/1000) # divide all counts by 1000, round to nearest - count_dead <- round(df_F1[i,]$Count_dead/1000) # divide all counts by 1000, round to nearest + #count_alive <- round(df_F1[i,]$Count_alive/1000) # divide all counts by 1000, round to nearest + #count_dead <- round(df_F1[i,]$Count_dead/1000) # divide all counts by 1000, round to nearest + + count_alive <- round(df_F1[i,]$Count_alive) # divide all counts by 1000, round to nearest + count_dead <- round(df_F1[i,]$Count_dead) # divide all counts by 1000, round to nearest loopDF <- data.frame(matrix(0,ncol = 6, nrow = (count_alive + count_dead))) colnames(loopDF) <- (c('Parent_pCO2','Larvae_pCO2','Generation','Tank','Age','Count_dead')) @@ -193,7 +198,7 @@ for (i in 1:nrow(df_F1)) { -# (2) Mean counts alive and dead by tank - better presentative for plots! +# (2) Mean counts alive and dead by tank - better representated for plots! # for loop for plotting!! # difference here is we take a mean for counts alive and dead within treatment, no tank included and a singluar line for plotting @@ -233,7 +238,7 @@ head(binary_df_F1_all) %>% kbl() %>% kable_classic(full_width = F, html_font = " ```{r F1 build surv object and plot em} - +?Surv # (1) Including tank - to use as surv object in coxme statistics (random factor included!) not surtable for plotting # build the surv object! # with tank included in the data frame (first for loop above) @@ -276,7 +281,7 @@ dev.off() ``` ```{r F1 run coxme on surv stats} - +?coxme head(binary_df_F1_all) # use this data from the first fr loop - has live and dead counts with reoslution of treatment AND replicate tank!!! @@ -295,7 +300,7 @@ summary(F1_survFIT_cox) #pairwise comparisons library(multcomp) - +?glht glht(F1_survFIT_cox, linfct = mcp(Larvae_pCO2 = "Dunnett"),alternative = "greater") # Linear Hypotheses: # Estimate @@ -708,8 +713,12 @@ binary_df_F3_mod <- data.frame() # start dataframe binary_df_F3_low <- data.frame() # start dataframe # run it for (i in 1:nrow(df_F3)) { - count_alive <- round(df_F3[i,]$Count_alive/1000) # divide all counts by 1000, round to nearest - count_dead <- round(df_F3[i,]$Count_dead/1000) # divide all counts by 1000, round to nearest + #count_alive <- round(df_F3[i,]$Count_alive/1000) # divide all counts by 1000, round to nearest + #count_dead <- round(df_F3[i,]$Count_dead/1000) # divide all counts by 1000, round to nearest + + count_alive <- round(df_F3[i,]$Count_alive) # divide all counts by 1000, round to nearest + count_dead <- round(df_F3[i,]$Count_dead) # divide all counts by 1000, round to nearest + loopDF <- data.frame(matrix(0,ncol = 6, nrow = (count_alive + count_dead))) colnames(loopDF) <- (c('Parent_pCO2','Larvae_pCO2','Generation','Tank','Age','Count_dead')) @@ -760,8 +769,12 @@ binary_df_F3_MEANS_low <- data.frame() # start dataframe # run it for (i in 1:nrow(df_F3_MEANS)) { - count_alive <- round(df_F3_MEANS[i,]$mean_Count_alive/1000) # divide all counts by 1000, round to nearest - count_dead <- round(df_F3_MEANS[i,]$mean_Count_dead/1000) # divide all counts by 1000, round to nearest + #count_alive <- round(df_F3_MEANS[i,]$mean_Count_alive/1000) # divide all counts by 1000, round to nearest + #count_dead <- round(df_F3_MEANS[i,]$mean_Count_dead/1000) # divide all counts by 1000, round to nearest + + count_alive <- round(df_F3_MEANS[i,]$mean_Count_alive) # divide all counts by 1000, round to nearest + count_dead <- round(df_F3_MEANS[i,]$mean_Count_dead) # divide all counts by 1000, round to nearest + loopDF <- data.frame(matrix(0,ncol = 4, nrow = (count_alive + count_dead))) colnames(loopDF) <- (c('Parent_pCO2','Larvae_pCO2','Age','Count_dead')) @@ -789,7 +802,7 @@ for (i in 1:nrow(df_F3_MEANS)) { ``` ```{r F3 build surv object} - +nrow(binary_df_F3_all) # create surv object # with tank included in the data frame (first for loop above) surv_F3_all <- Surv(time = as.numeric(binary_df_F3_all$Age), @@ -804,6 +817,32 @@ surv_F3_mod <- Surv(time = as.numeric(binary_df_F3_mod$Age), surv_F3_low <- Surv(time = as.numeric(binary_df_F3_low$Age), event = as.numeric(binary_df_F3_low$Count_dead), type = "right") + +F3_low_survplot <- ggsurvplot(survfit(surv_F3_low ~ Larvae_pCO2, binary_df_F3_low), + risk.table = F, pval = F , conf.int = TRUE, + palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8, 0.3)), + ggtheme = theme_survminer(), + break.x.by = 2 ,xlab = "Age (days post fertilization)", + legend.title = "treatment") + +F3_mod_survplot <- ggsurvplot(survfit(surv_F3_mod ~ Larvae_pCO2, binary_df_F3_mod), + risk.table = F, pval = F , conf.int = TRUE, + palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8, 0.3)), + ggtheme = theme_survminer(), + break.x.by = 2 ,xlab = "Age (days post fertilization)", + legend.title = "treatment") + +F3_high_survplot <- ggsurvplot(survfit(surv_F3_high ~ Larvae_pCO2, binary_df_F3_high), + risk.table = F, pval = F , conf.int = TRUE, + palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8, 0.3)), + ggtheme = theme_survminer(), + break.x.by = 2 ,xlab = "Age (days post fertilization)", + legend.title = "treatment") + +F3_full_plot <- ggarrange(F3_low_survplot, F3_mod_survplot, F3_high_survplot, ncol = 3) +print(F3_full_plot) + + # create surv object # without tank included - means wihtin treatment SEPCIFICALLY FOR PLOTTING second for loop above surv_F3_MEANS <- Surv(time = as.numeric(binary_df_F3_MEANS$Age), @@ -836,6 +875,15 @@ surv_F3_MEANS_high <- Surv(time = as.numeric(binary_df_F3_MEANS_high$Age), ```{r F3 plot em} +pdf("C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Survival/F3_survival.pdf", width = 10) +ggsurvplot(survfit(surv_F3_all ~ Larvae_pCO2, binary_df_F3_all), + risk.table = F, pval = F , conf.int = TRUE, + palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.2)), + break.x.by = 2 ,xlab = "Age (days post fertilization)", + facet.by = "Parent_pCO2", + legend.title = "treatment") +dev.off() + pdf("C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Survival/F3_survival.pdf", width = 10) ggsurvplot(survfit(surv_F3_MEANS ~ Larvae_pCO2, binary_df_F3_MEANS), risk.table = F, pval = F , conf.int = TRUE, @@ -846,19 +894,19 @@ ggsurvplot(survfit(surv_F3_MEANS ~ Larvae_pCO2, binary_df_F3_MEANS), dev.off() # # pdf("C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Survival/F3_surival_low_history.pdf") -# ggsurvplot(survfit(surv_F3_MEANS_low ~ Larvae_pCO2, binary_df_F3_MEANS_low), -# risk.table = F, pval = F , conf.int = TRUE, -# palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.3)), -# break.x.by = 2 ,xlab = "Age (days post fertilization)", -# legend.title = "treatment") + ggsurvplot(survfit(surv_F3_MEANS_low ~ Larvae_pCO2, binary_df_F3_MEANS_low), + risk.table = F, pval = F , conf.int = TRUE, + palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.3)), + break.x.by = 2 ,xlab = "Age (days post fertilization)", + legend.title = "treatment") # dev.off() # # pdf("C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Survival/F3_moderate_history.pdf") -# ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), -# risk.table = F, pval = F , conf.int = TRUE, -# palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.2)), -# break.x.by = 2 ,xlab = "Age (days post fertilization)", -# legend.title = "treatment") +ggsurvplot(survfit(surv_F3_MEANS_mod ~ Larvae_pCO2, binary_df_F3_MEANS_mod), + risk.table = F, pval = F , conf.int = TRUE, + palette = alpha(c("forestgreen","orange", "purple"), c(1,0.8,0.2)), + break.x.by = 2 ,xlab = "Age (days post fertilization)", + legend.title = "treatment") # dev.off() # # pdf("C:/Users/samjg/Documents/Github_repositories/EAD-ASEB-Airradians_Larvae_crossgen_OA/RAnalysis/Output/Survival/F3_high_history.pdf") @@ -893,15 +941,15 @@ summary(survFitME_all_tank) # Larvae_pCO2High pCO2:Parent_pCO2Moderate pCO2 0.008771872 1.0088105 0.02630283 0.33 # Larvae_pCO2Moderate pCO2:Parent_pCO2High pCO2 0.035694517 1.0363392 0.02605141 1.37 # Larvae_pCO2High pCO2:Parent_pCO2High pCO2 0.006232918 1.0062524 0.02611643 0.24 -p - 1.7e-01 - 6.0e-01 - 4.5e-02 - 3.8e-05 - 2.5e-01 - 7.4e-01 - 1.7e-01 - 8.1e-01 +# p +# 1.7e-01 +# 6.0e-01 +# 4.5e-02 +# 3.8e-05 +# 2.5e-01 +# 7.4e-01 +# 1.7e-01 +# 8.1e-01 #The standard deviation of the tank random effect is very small (0.015) which suggests that tank does not really affect the outcome. Also the liklihood ratio test (Chi?sq) value is the same (527597), meaning that including tank as a random effect does not improve the model. We likely do not need to include it as a random effect.