Comparing factors that influence oogonia shedding vs retention.

For each run: I take the average unit conditions for salinity, pH, temperature, and dissolved oxygen. For all parameters except pH I take the average of Hach reading taken in unit during an experimental run. For pH, I use a corrected pH value (correcting the pH value read by the Hach by the actual pH value analized in lab). I average the oogonia on the discs for each unit. I average the oogonia left in each experimental thalli.

I have the linear model that I ran preivously and then AIC ranking.

To do: put data on GitHub

Read in data

exp<-read.csv("C:/Users/chels/Box Sync/Thesis/Data/Working data/Final data/wegener_experiment.csv", header = TRUE, sep=",", fileEncoding="UTF-8-BOM", stringsAsFactors = FALSE)

Convert to number

exp$oo_disc_1<-as.numeric(exp$oo_disc_1)
exp$oo_disc_2<-as.numeric(exp$oo_disc_2)
exp$oo_disc_3<-as.numeric(exp$oo_disc_3)
exp$oo_disc_4<-as.numeric(exp$oo_disc_4)

exp$oo_in_1<-as.numeric(exp$oo_in_1)
exp$oo_in_2<-as.numeric(exp$oo_in_2)
exp$oo_in_3<-as.numeric(exp$oo_in_3)

exp$unit_temp_hach<-as.numeric(exp$unit_temp_hach)
exp$unit_do_hach<-as.numeric(exp$unit_do_hach)
exp$treatment<-as.character(exp$treatment)

Average

#average oogonia released per thalli/unit
exp$avg_oo_thalli<-rowMeans(exp[,c('oo_disc_1', 'oo_disc_2', 'oo_disc_3','oo_disc_4')], na.rm=TRUE)

#average oogonia remaining in experimental thalli
exp$avg_oo_in<-rowMeans(exp[,c('oo_in_1', 'oo_in_2', 'oo_in_3')], na.rm=TRUE)

Averaging conditions per unit per run

#subset per run
r1<-subset(exp, run=='1')
r2<-subset(exp, run=='2')
r3<-subset(exp,run=='3')

#average of conditions per thalli/unit
#temperature
r1tempav<-aggregate( unit_temp_hach ~ thalli,r1, mean)
r2tempav<-aggregate( unit_temp_hach ~ thalli,r2, mean)
r3tempav<-aggregate( unit_temp_hach ~ thalli,r3, mean)

#dissolved oxygen
r1doav<-aggregate( unit_do_hach ~ thalli,r1, mean)
r2doav<-aggregate( unit_do_hach ~ thalli,r2, mean)
r3doav<-aggregate( unit_do_hach ~ thalli,r3, mean)

#salinity
r1salav<-aggregate( unit_sal_hach ~ thalli,r1, mean)
r2salav<-aggregate( unit_sal_hach ~ thalli,r2, mean)
r3salav<-aggregate( unit_sal_hach ~ thalli,r3, mean)

#pH
r1phav<-aggregate( unit_ph_c ~ thalli,r1, mean)
r2phav<-aggregate( unit_ph_c ~ thalli,r2, mean)
r3phav<-aggregate( unit_ph_c ~ thalli,r3, mean)

Clean up data for modeling and merging

#Remove rows with NA and make new dataframe
oo<-exp[!is.na(exp$avg_oo_thalli), ]

#select columns to keep
oo<-data.frame(oo$thalli, oo$avg_oo_thalli, oo$treatment, oo$avg_oo_in)

#add run number back
r1tempav$run<-'1'
r2tempav$run<-'2'
r3tempav$run<-'3'

#combine dataframes by run 
r1complete<-cbind(r1tempav, r1doav, r1salav, r1phav)
r2complete<-cbind(r2tempav, r2doav, r2salav, r2phav)
r3complete<-cbind(r3tempav, r3doav, r3salav, r3phav)

#combine
runs<-rbind(r1complete, r2complete, r3complete)

#remove redunant columns
runs<-runs[-c(4, 6, 8)]

#combine with oogonia datafram
final<-cbind(oo,runs)

#remove redundant columns
final<-final[-c(5)]

#rename columns
names(final)[1] <- "thalli"
names(final)[2] <- "avg.oo.disc"
names(final)[3] <- "treatment"
names(final)[4] <- "avg.oo.in"

#clean up environment
rm(list=setdiff(ls(), "final"))

Model

#oogonia released
m1<- lmer(avg.oo.disc ~ unit_sal_hach + unit_ph_c + unit_do_hach + unit_temp_hach + (1|run), data=final)
## boundary (singular) fit: see ?isSingular
plot(m1)

m1
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## avg.oo.disc ~ unit_sal_hach + unit_ph_c + unit_do_hach + unit_temp_hach +  
##     (1 | run)
##    Data: final
## REML criterion at convergence: 805.9728
## Random effects:
##  Groups   Name        Std.Dev.
##  run      (Intercept)   0.0   
##  Residual             314.2   
## Number of obs: 60, groups:  run, 3
## Fixed Effects:
##    (Intercept)   unit_sal_hach       unit_ph_c    unit_do_hach  unit_temp_hach  
##        4949.32          -50.08          -54.51          -61.45         -140.46  
## optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## avg.oo.disc ~ unit_sal_hach + unit_ph_c + unit_do_hach + unit_temp_hach +  
##     (1 | run)
##    Data: final
## 
## REML criterion at convergence: 806
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5738 -0.5391 -0.1981  0.3392  3.9545 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  run      (Intercept)     0      0.0   
##  Residual             98744    314.2   
## Number of obs: 60, groups:  run, 3
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)     4949.32    2883.81   1.716
## unit_sal_hach    -50.08      21.47  -2.333
## unit_ph_c        -54.51      97.44  -0.559
## unit_do_hach     -61.45     127.49  -0.482
## unit_temp_hach  -140.46     107.63  -1.305
## 
## Correlation of Fixed Effects:
##             (Intr) unt_s_ unt_p_ unt_d_
## unit_sl_hch -0.028                     
## unit_ph_c   -0.130 -0.405              
## unit_do_hch -0.895 -0.020 -0.172       
## unt_tmp_hch -0.953 -0.122  0.059  0.843
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
anova(m1)
## Analysis of Variance Table
##                npar Sum Sq Mean Sq F value
## unit_sal_hach     1 854078  854078  8.6494
## unit_ph_c         1      6       6  0.0001
## unit_do_hach      1 130649  130649  1.3231
## unit_temp_hach    1 168181  168181  1.7032
confint(m1)
## Computing profile confidence intervals ...
##                     2.5 %       97.5 %
## .sig01            0.00000   129.766903
## .sigma          254.11386   363.905484
## (Intercept)    -550.02569 10553.071124
## unit_sal_hach   -91.01707    -9.146142
## unit_ph_c      -240.32627   131.314420
## unit_do_hach   -306.97551   181.677930
## unit_temp_hach -355.34858    64.783901
tab_model(m1)
  avg.oo.disc
Predictors Estimates CI p
(Intercept) 4949.32 -702.85 – 10601.49 0.086
unit_sal_hach -50.08 -92.15 – -8.01 0.020
unit_ph_c -54.51 -245.49 – 136.48 0.576
unit_do_hach -61.45 -311.33 – 188.43 0.630
unit_temp_hach -140.46 -351.41 – 70.49 0.192
Random Effects
σ2 98743.93
τ00 run 0.00
N run 3
Observations 60
Marginal R2 / Conditional R2 0.165 / NA 
#oogonia retained
m2<- lmer(avg.oo.in ~ unit_sal_hach + unit_ph_c + unit_do_hach + unit_temp_hach + (1|run), data=final)
## boundary (singular) fit: see ?isSingular
plot(m2)

m2
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## avg.oo.in ~ unit_sal_hach + unit_ph_c + unit_do_hach + unit_temp_hach +  
##     (1 | run)
##    Data: final
## REML criterion at convergence: 434.3046
## Random effects:
##  Groups   Name        Std.Dev.
##  run      (Intercept)  0.00   
##  Residual             10.71   
## Number of obs: 60, groups:  run, 3
## Fixed Effects:
##    (Intercept)   unit_sal_hach       unit_ph_c    unit_do_hach  unit_temp_hach  
##        225.809           1.249          -6.329          -7.711          -6.288  
## optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
summary(m2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## avg.oo.in ~ unit_sal_hach + unit_ph_c + unit_do_hach + unit_temp_hach +  
##     (1 | run)
##    Data: final
## 
## REML criterion at convergence: 434.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.22097 -0.65458  0.07002  0.64055  2.07414 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  run      (Intercept)   0.0     0.00   
##  Residual             114.7    10.71   
## Number of obs: 60, groups:  run, 3
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)    225.8091    98.3042   2.297
## unit_sal_hach    1.2485     0.7318   1.706
## unit_ph_c       -6.3290     3.3217  -1.905
## unit_do_hach    -7.7113     4.3461  -1.774
## unit_temp_hach  -6.2880     3.6689  -1.714
## 
## Correlation of Fixed Effects:
##             (Intr) unt_s_ unt_p_ unt_d_
## unit_sl_hch -0.028                     
## unit_ph_c   -0.130 -0.405              
## unit_do_hch -0.895 -0.020 -0.172       
## unt_tmp_hch -0.953 -0.122  0.059  0.843
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
anova(m2)
## Analysis of Variance Table
##                npar Sum Sq Mean Sq F value
## unit_sal_hach     1  80.25   80.25  0.6994
## unit_ph_c         1 587.76  587.76  5.1225
## unit_do_hach      1  43.02   43.02  0.3750
## unit_temp_hach    1 337.04  337.04  2.9374
confint(m2)
## Computing profile confidence intervals ...
## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep
## Warning in FUN(X[[i]], ...): non-monotonic profile for .sig01
## Warning in confint.thpr(pp, level = level, zeta = zeta): bad spline fit
## for .sig01: falling back to linear interpolation
##                      2.5 %       97.5 %
## .sig01           0.0000000 3.797932e+00
## .sigma           8.6623059 1.240491e+01
## (Intercept)     38.3468551 4.132714e+02
## unit_sal_hach   -0.1468886 2.643951e+00
## unit_ph_c      -12.6633370 5.257105e-03
## unit_do_hach   -15.9991122 5.764360e-01
## unit_temp_hach -13.2843854 7.083806e-01
tab_model(m2)
  avg.oo.in
Predictors Estimates CI p
(Intercept) 225.81 33.14 – 418.48 0.022
unit_sal_hach 1.25 -0.19 – 2.68 0.088
unit_ph_c -6.33 -12.84 – 0.18 0.057
unit_do_hach -7.71 -16.23 – 0.81 0.076
unit_temp_hach -6.29 -13.48 – 0.90 0.087
Random Effects
σ2 114.74
τ00 run 0.00
N run 3
Observations 60
Marginal R2 / Conditional R2 0.134 / NA

AIC Ranking

Oogonia released grouped by run

#models
m1 <- lmer(avg.oo.disc ~ unit_sal_hach + unit_ph_c + unit_do_hach + unit_temp_hach + unit_sal_hach:unit_ph_c + (1|run), data=final)
## boundary (singular) fit: see ?isSingular
m2 <- lmer(avg.oo.disc ~ unit_sal_hach + unit_ph_c + unit_temp_hach + unit_sal_hach:unit_ph_c + (1|run), data=final)
## boundary (singular) fit: see ?isSingular
m3 <- lmer(avg.oo.disc ~ unit_sal_hach + unit_ph_c + unit_do_hach + unit_sal_hach:unit_ph_c + (1|run), data=final)
## boundary (singular) fit: see ?isSingular
m4 <- lmer(avg.oo.disc ~ unit_sal_hach + unit_ph_c + unit_sal_hach:unit_ph_c + (1|run), data=final)
## boundary (singular) fit: see ?isSingular
m5 <- lmer(avg.oo.disc ~ unit_sal_hach + unit_ph_c + (1|run), data=final)
## boundary (singular) fit: see ?isSingular
m6 <- lmer(avg.oo.disc ~ unit_sal_hach + (1|run), data=final)
## boundary (singular) fit: see ?isSingular
m7 <- lmer(avg.oo.disc ~ unit_ph_c + (1|run), data=final)

m8 <- lmer(avg.oo.disc ~ unit_temp_hach + (1|run), data=final)
## boundary (singular) fit: see ?isSingular
m9 <- lmer(avg.oo.disc ~ unit_do_hach + (1|run), data=final)
## boundary (singular) fit: see ?isSingular
#define list of models
models <- list(m1,m2,m3,m4,m5,m6,m7,m8,m9)

#specify model names- just keeping them the same
mod.names <- c('m1','m2','m3','m4', 'm5','m6','m7','m8', 'm9')

#calculate AIC of each model
aictab(cand.set = models, modnames = mod.names)
## Warning in aictab.AIClmerMod(cand.set = models, modnames = mod.names): 
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects
## 
## Model selection based on AICc:
## 
##    K   AICc Delta_AICc AICcWt Cum.Wt  Res.LL
## m1 8 814.05       0.00   0.98   0.98 -397.61
## m2 7 823.14       9.09   0.01   0.99 -403.49
## m3 7 824.39      10.34   0.01   1.00 -404.12
## m4 6 833.55      19.49   0.00   1.00 -409.98
## m5 5 841.59      27.54   0.00   1.00 -415.24
## m6 4 850.02      35.97   0.00   1.00 -420.65
## m8 4 853.61      39.55   0.00   1.00 -422.44
## m7 4 854.05      40.00   0.00   1.00 -422.66
## m9 4 854.86      40.81   0.00   1.00 -423.07

Table of top ranked model (avg.oo.disc ~ unit_sal_hach + unit_ph_c + unit_do_hach + unit_temp_hach + unit_sal_hach:unit_ph_c + (1|run), data=final)

tab_model(m1)
  avg.oo.disc
Predictors Estimates CI p
(Intercept) -1769.73 -35418.20 – 31878.74 0.918
unit_sal_hach 199.35 -1032.47 – 1431.17 0.751
unit_ph_c 806.81 -3448.71 – 5062.33 0.710
unit_do_hach -59.30 -311.35 – 192.74 0.645
unit_temp_hach -146.37 -360.93 – 68.20 0.181
unit_sal_hach * unit_ph_c -31.63 -187.72 – 124.47 0.691
Random Effects
σ2 100279.69
τ00 run 0.00
N run 3
Observations 60
Marginal R2 / Conditional R2 0.165 / NA 
tab_model(m1,
  show.ci = 0.95,
  show.se = TRUE,
  show.stat = FALSE,
  emph.p = TRUE,
  df.method = "kr",
  p.style = "numeric_stars",
  title = "AIC top ranked models for experiment",
  linebreak = TRUE,
  dv.labels = c(
    "oogonia released"
      ),
  pred.labels = c("(Intercept)", "salinity", "pH", "dissolved oxygen", "water temperature", "salinity:pH"),
  CSS = css_theme("cells")
  )
AIC top ranked models for experiment
  oogonia released
Predictors Estimates std. Error CI p
(Intercept) -1769.73 17818.70 -37494.11 – 33954.64 0.921
salinity 199.35 628.89 -1062.06 – 1460.76 0.753
pH 806.81 2171.90 -3549.53 – 5163.16 0.712
dissolved oxygen -59.30 260.62 -1621.33 – 1502.72 0.847
water temperature -146.37 220.04 -2048.75 – 1756.02 0.611
salinity:pH -31.63 79.67 -191.43 – 128.17 0.693
Random Effects
σ2 100279.69
τ00 run 0.00
N run 3
Observations 60
Marginal R2 / Conditional R2 0.165 / NA
  • p<0.05   ** p<0.01   *** p<0.001

Oogonia retained grouped by run

#models
m1 <- lmer(avg.oo.in ~ unit_sal_hach + unit_ph_c + unit_do_hach + unit_temp_hach + unit_sal_hach:unit_ph_c + (1|run), data=final)
## boundary (singular) fit: see ?isSingular
m2 <- lmer(avg.oo.in ~ unit_sal_hach + unit_ph_c + unit_temp_hach + unit_sal_hach:unit_ph_c + (1|run), data=final)

m3 <- lmer(avg.oo.in ~ unit_sal_hach + unit_ph_c + unit_do_hach + unit_sal_hach:unit_ph_c + (1|run), data=final)

m4 <- lmer(avg.oo.in ~ unit_sal_hach + unit_ph_c + unit_sal_hach:unit_ph_c + (1|run), data=final)

m5 <- lmer(avg.oo.in ~ unit_sal_hach + unit_ph_c + (1|run), data=final)
## boundary (singular) fit: see ?isSingular
m6 <- lmer(avg.oo.in ~ unit_sal_hach + (1|run), data=final)

m7 <- lmer(avg.oo.in ~ unit_ph_c + (1|run), data=final)

m8 <- lmer(avg.oo.in ~ unit_temp_hach + (1|run), data=final)

m9 <- lmer(avg.oo.in ~ unit_do_hach + (1|run), data=final)

#define list of models
models <- list(m1,m2,m3,m4,m5,m6,m7,m8,m9)

#specify model names- just keeping them the same
mod.names <- c('m1','m2','m3','m4', 'm5','m6','m7','m8', 'm9')

#calculate AIC of each model
aictab(cand.set = models, modnames = mod.names)
## Warning in aictab.AIClmerMod(cand.set = models, modnames = mod.names): 
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects
## 
## Model selection based on AICc:
## 
##    K   AICc Delta_AICc AICcWt Cum.Wt  Res.LL
## m1 8 449.00       0.00   0.79   0.79 -215.09
## m3 7 453.66       4.66   0.08   0.87 -218.75
## m2 7 454.01       5.01   0.06   0.93 -218.93
## m4 6 455.19       6.19   0.04   0.97 -220.80
## m5 5 456.69       7.69   0.02   0.98 -222.79
## m7 4 458.05       9.05   0.01   0.99 -224.66
## m9 4 459.43      10.43   0.00   1.00 -225.35
## m8 4 461.05      12.05   0.00   1.00 -226.16
## m6 4 463.10      14.10   0.00   1.00 -227.19

Table of top ranked model (avg.oo.disc ~ unit_sal_hach + unit_ph_c + unit_do_hach + unit_temp_hach + unit_sal_hach:unit_ph_c + (1|run), data=final). Same model as oogonia released

tab_model(m1)
  avg.oo.in
Predictors Estimates CI p
(Intercept) -90.61 -1236.11 – 1054.88 0.877
unit_sal_hach 12.99 -28.94 – 54.93 0.544
unit_ph_c 34.23 -110.64 – 179.10 0.643
unit_do_hach -7.61 -16.19 – 0.97 0.082
unit_temp_hach -6.57 -13.87 – 0.74 0.078
unit_sal_hach * unit_ph_c -1.49 -6.80 – 3.82 0.583
Random Effects
σ2 116.22
τ00 run 0.00
N run 3
Observations 60
Marginal R2 / Conditional R2 0.136 / NA 
tab_model(m1,
  show.ci = 0.95,
  show.se = TRUE,
  show.stat = FALSE,
  emph.p = TRUE,
  df.method = "kr",
  p.style = "numeric_stars",
  title = "AIC top ranked models for experiment",
  linebreak = TRUE,
  dv.labels = c(
    "oogonia retained"
      ),
  pred.labels = c("(Intercept)", "salinity", "pH", "dissolved oxygen", "water temperature", "salinity:pH"),
  CSS = css_theme("cells")
  )
AIC top ranked models for experiment
  oogonia retained
Predictors Estimates std. Error CI p
(Intercept) -90.61 606.60 -1306.78 – 1125.55 0.882
salinity 12.99 21.41 -29.95 – 55.94 0.546
pH 34.23 73.94 -114.07 – 182.54 0.645
dissolved oxygen -7.61 8.87 -60.79 – 45.57 0.506
water temperature -6.57 7.49 -71.33 – 58.20 0.521
salinity:pH -1.49 2.71 -6.93 – 3.95 0.585
Random Effects
σ2 116.22
τ00 run 0.00
N run 3
Observations 60
Marginal R2 / Conditional R2 0.136 / NA
  • p<0.05   ** p<0.01   *** p<0.001