Set up

read in relevant files

I took the sheet Diana sent me and saved a csv with each tab of that sheet then read those in here.

remove dead and lost

We only had 52 removed. This seems a bit low to me and I wonder if the scallops that were removed weren’t marked in this? or if something else strange happened. But rolling with this for now.

## [1] 425  16
## [1] 373  16

look at data

## Looking for trends

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

add a pH where minimum pH is 0 (to avoid weird intercepts later on)

Not totally sure if this is necessary for temperature too, but it makes sense to me.

## [1] 7.36
## [1] 6.14

# Look at all species together Not sure this really matters much, or makes much sense but I thought I would look for overall trends. ## LMER Best fit model included random effect of tank (intercept) and average pH but the residuals etc. were wacky. You can see in the plot that some species (juv arctica, scallop) do not seem to follow the trend.

## boundary (singular) fit: see help('isSingular')
## Backward reduced random-effect table:
## 
##            Eliminated npar  logLik    AIC         LRT Df Pr(>Chisq)
## <none>                   6 -1772.8 3557.6                          
## (1 | Tank)          1    5 -1772.8 3555.6 -5.9572e-10  1          1
## 
## Backward reduced fixed-effect table:
##                         Eliminated Df Sum of Sq    RSS    AIC F value    Pr(>F)
## temp.average:pH.average          1  1     386.3 304800 2507.3  0.4683    0.4942
## temp.average                     2  1     158.8 304959 2505.5  0.1928    0.6608
## pH.average                       0  1   23858.4 328818 2531.6 29.0251 1.272e-07
##                            
## temp.average:pH.average    
## temp.average               
## pH.average              ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## mean ~ pH.average
## 
## Call:
## lm(formula = mean ~ pH.average, data = color.all)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -77.019  -9.765   8.429  20.873  48.838 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  380.526     49.652   7.664 1.59e-13 ***
## pH.average   -34.679      6.437  -5.387 1.27e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28.67 on 371 degrees of freedom
## Multiple R-squared:  0.07256,    Adjusted R-squared:  0.07006 
## F-statistic: 29.03 on 1 and 371 DF,  p-value: 1.272e-07

## [1] 3565.999
  mean
Predictors Estimates CI p
(Intercept) 380.53 282.89 – 478.16 <0.001
pH average -34.68 -47.34 – -22.02 <0.001
Observations 373
R2 / R2 adjusted 0.073 / 0.070

LMER again but NORMALIZED

In the above model the intercepts weren’t making much sense to me, which when you think about it makes sense. In the models above the intercept is at X=0 but a pH or temperature of 0 in our study system has no meaning. I therefore modeled here with the pH and temperature averages, normalized to the lowest occuring average.

You can see these generally have the same output, just with shifted intercepts.

## boundary (singular) fit: see help('isSingular')
## Backward reduced random-effect table:
## 
##            Eliminated npar  logLik    AIC LRT Df Pr(>Chisq)
## <none>                   6 -1772.8 3557.6                  
## (1 | Tank)          1    5 -1772.8 3555.6   0  1          1
## 
## Backward reduced fixed-effect table:
##                               Eliminated Df Sum of Sq    RSS    AIC F value
## temp.normalized:pH.normalized          1  1     386.3 304800 2507.3  0.4683
## temp.normalized                        2  1     158.8 304959 2505.5  0.1928
## pH.normalized                          0  1   23858.4 328818 2531.6 29.0251
##                                  Pr(>F)    
## temp.normalized:pH.normalized    0.4942    
## temp.normalized                  0.6608    
## pH.normalized                 1.272e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## mean ~ pH.normalized
## 
## Call:
## lm(formula = mean ~ pH.normalized, data = color.all)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -77.019  -9.765   8.429  20.873  48.838 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    125.287      2.699  46.419  < 2e-16 ***
## pH.normalized  -34.679      6.437  -5.387 1.27e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28.67 on 371 degrees of freedom
## Multiple R-squared:  0.07256,    Adjusted R-squared:  0.07006 
## F-statistic: 29.03 on 1 and 371 DF,  p-value: 1.272e-07

## [1] 3565.999
  mean mean
Predictors Estimates CI p Estimates CI p
(Intercept) 380.53 282.89 – 478.16 <0.001 125.29 119.98 – 130.59 <0.001
pH average -34.68 -47.34 – -22.02 <0.001
pH normalized -34.68 -47.34 – -22.02 <0.001
Observations 373 373
R2 / R2 adjusted 0.073 / 0.070 0.073 / 0.070

ANOVA to see if same signal is retained

ANOVA is really not appropriate here because our population/distribution is so wacky. But it does show the same trend (pH significant), although type III ANOVA which accounts for the unbalanced design of the expiriment, brings the significance of pH down to 0.052.

## Call:
##    aov(formula = mean ~ pH * temp, data = color.all)
## 
## Terms:
##                        pH      temp   pH:temp Residuals
## Sum of Squares   27296.93    402.66   2690.29 298427.76
## Deg. of Freedom         3         2         6       361
## 
## Residual standard error: 28.75186
## Estimated effects may be unbalanced
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## pH            3  27297    9099  11.007 6.23e-07 ***
## temp          2    403     201   0.244    0.784    
## pH:temp       6   2690     448   0.542    0.776    
## Residuals   361 298428     827                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Design is unbalanced - use se.contrast() for se's
## Tables of means
## Grand mean
##          
## 113.1428 
## 
##  pH 
##       7.4   7.6   7.8     8
##     126.3 115.3 104.1 106.9
## rep  92.0  95.0  88.0  98.0
## 
##  temp 
##         6     9  12
##     111.5 113.7 114
## rep 103.0 185.0  85
## 
##  pH:temp 
##      temp
## pH    6      9      12    
##   7.4 124.31 128.18 123.99
##   rep  25.00  48.00  19.00
##   7.6 111.38 114.15 122.00
##   rep  27.00  45.00  23.00
##   7.8 102.63 107.32  99.59
##   rep  24.00  42.00  22.00
##   8   107.50 105.38 109.65
##   rep  27.00  50.00  21.00
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mean ~ pH * temp, data = color.all)
## 
## $pH
##               diff        lwr         upr     p adj
## 7.6-7.4 -10.998892 -21.853779  -0.1440043 0.0456879
## 7.8-7.4 -22.157394 -33.222645 -11.0921438 0.0000023
## 8-7.4   -19.384413 -30.157250  -8.6115757 0.0000283
## 7.8-7.6 -11.158502 -22.138005  -0.1790002 0.0447465
## 8-7.6    -8.385521 -19.070263   2.2992215 0.1804570
## 8-7.8     2.772981  -8.125408  13.6713714 0.9131445

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group  11   0.308  0.984
##       361

## Anova Table (Type III tests)
## 
## Response: mean
##             Sum Sq  Df  F value  Pr(>F)    
## (Intercept) 386338   1 467.3427 < 2e-16 ***
## pH            6423   3   2.5898 0.05269 .  
## temp           370   2   0.2238 0.79960    
## pH:temp       2690   6   0.5424 0.77588    
## Residuals   298428 361                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Shapiro-Wilk normality test
## 
## data:  aov_residuals
## W = 0.87162, p-value < 2.2e-16

Plot

## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following object is masked from 'package:purrr':
## 
##     compact
## The following objects are masked from 'package:plotly':
## 
##     arrange, mutate, rename, summarise
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following object is masked from 'package:ggpubr':
## 
##     mutate
##    temp  pH  N      mean       sd       se        ci
## 1     6 7.4 25 124.31220 30.77729 6.155459 10.531264
## 2     6 7.6 27 111.38259 27.62450 5.316337  9.067639
## 3     6 7.8 24 102.63350 25.12788 5.129208  8.790803
## 4     6   8 27 107.49989 24.26398 4.669605  7.964562
## 5     9 7.4 48 128.18169 33.99215 4.906344  8.232486
## 6     9 7.6 45 114.15244 29.07736 4.334596  7.283119
## 7     9 7.8 42 107.31557 26.61527 4.106825  6.911286
## 8     9   8 50 105.38400 26.71152 3.777579  6.333303
## 9    12 7.4 19 123.99026 29.04887 6.664268 11.556265
## 10   12 7.6 23 122.00217 34.74127 7.244055 12.439089
## 11   12 7.8 22  99.58973 25.71482 5.482418  9.433832
## 12   12   8 21 109.64567 27.27758 5.952455 10.266308

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

so there are probably some issues with this because it is not a normal distribution, but earlier I saw that juv arctica were the darkest ~ 60 and I suspect if we took juv arctica out it would be normal… or at least more normal. I will try modeling by each species now.

Individual Species Groups

mercenaria

check structure

Looks pretty normal to me!

## Classes 'data.table' and 'data.frame':   122 obs. of  18 variables:
##  $ Tank           : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ File path      : chr  NA "C:\\Users\\thatc\\OneDrive - Iowa State University\\Desktop\\jpg\\mercenaria\\b38.jpeg" NA "C:\\Users\\thatc\\OneDrive - Iowa State University\\Desktop\\jpg\\mercenaria\\2nd round\\b40.jpeg" ...
##  $ file name      : chr  "b38.jpeg" "b39" "b40.jpeg" "b41.jpeg" ...
##  $ shell          : chr  "b38" "b39" "b40" "b41" ...
##  $ N              : num  456633 316471 207098 195915 274746 ...
##  $ mean           : num  146 142 134 133 140 ...
##  $ stdev          : num  24.1 22.7 20.8 23 16.8 ...
##  $ mode           : num  139 128 126 137 144 130 133 134 145 122 ...
##  $ pH             : Factor w/ 4 levels "7.4","7.6","7.8",..: 2 2 2 2 2 2 2 2 2 4 ...
##  $ temp           : Factor w/ 3 levels "6","9","12": 3 3 3 3 3 3 3 3 3 2 ...
##  $ dead?          : chr  NA NA NA NA ...
##  $ species        : chr  "mercenaria" "mercenaria" "mercenaria" "mercenaria" ...
##  $ temp.average   : num  12.1 12.1 12.1 12.1 12.1 ...
##  $ temp.stdev     : num  0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.63 ...
##  $ pH.average     : num  7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 8 ...
##  $ pH.stdev       : num  0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.07 ...
##  $ pH.normalized  : num  0.21 0.21 0.21 0.21 0.21 0.21 0.21 0.21 0.21 0.64 ...
##  $ temp.normalized: num  5.92 5.92 5.92 5.92 5.92 5.92 5.92 5.92 5.92 3.05 ...
##  - attr(*, ".internal.selfref")=<externalptr>

## 
##  Shapiro-Wilk normality test
## 
## data:  mercenaria$mean
## W = 0.9925, p-value = 0.758

linear model

The best fit linear model included the random effect of tank (intercept) and both temperature and pH.

## Backward reduced random-effect table:
## 
##            Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                   6 -418.07 848.15                         
## (1 | Tank)          0    5 -428.99 867.98 21.838  1  2.967e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                               Eliminated Sum Sq Mean Sq NumDF  DenDF F value
## temp.normalized:pH.normalized          1   1.77    1.77     1 10.983  0.0342
## temp.normalized                        0 721.68  721.68     1 11.859 13.9285
## pH.normalized                          0 672.77  672.77     1 12.061 12.9845
##                                 Pr(>F)   
## temp.normalized:pH.normalized 0.856749   
## temp.normalized               0.002920 **
## pH.normalized                 0.003595 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## mean ~ temp.normalized + pH.normalized + (1 | Tank)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean ~ temp.normalized + pH.normalized + (1 | Tank)
##    Data: mercenaria
## 
## REML criterion at convergence: 840.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0419 -0.6219 -0.0391  0.5781  3.6860 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 25.99    5.098   
##  Residual             51.81    7.198   
## Number of obs: 122, groups:  Tank, 16
## 
## Fixed effects:
##                 Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)     121.6203     3.3600  11.7627  36.197 1.99e-13 ***
## temp.normalized   2.6750     0.7168  11.8589   3.732  0.00292 ** 
## pH.normalized   -22.5901     6.2691  12.0610  -3.603  0.00360 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.nr
## temp.nrmlzd -0.632       
## pH.normalzd -0.656  0.017

## [1] 850.3746
  mean
Predictors Estimates CI p
(Intercept) 121.62 114.97 – 128.27 <0.001
temp normalized 2.68 1.26 – 4.09 <0.001
pH normalized -22.59 -35.01 – -10.17 <0.001
Random Effects
σ2 51.81
Ï„00 Tank 25.99
ICC 0.33
N Tank 16
Observations 122
Marginal R2 / Conditional R2 0.430 / 0.621

ANOVA

The ANOVA results show similar pattern to the linear model but with interaction term on the edge (0.045). Tukey showed no significant difference between 7.6-7.4 & 8-7.8. All pH groups were different, with the smallest difference between 12 & 9. Since the design is unbalanced I also did a type III ANOVA which showed significant pH and pH temp interaction but no significant temperature term. I trust this more as it is accounting for the balance (or lack there of) in our design. Passed normality of residuals and homogenity of variance.

## Call:
##    aov(formula = mean ~ pH * temp, data = mercenaria)
## 
## Terms:
##                       pH     temp  pH:temp Residuals
## Sum of Squares  3817.842 4255.010  842.373  6939.814
## Deg. of Freedom        3        2        6       110
## 
## Residual standard error: 7.942872
## Estimated effects may be unbalanced
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## pH            3   3818  1272.6  20.172 1.73e-10 ***
## temp          2   4255  2127.5  33.722 3.79e-12 ***
## pH:temp       6    842   140.4   2.225   0.0459 *  
## Residuals   110   6940    63.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Design is unbalanced - use se.contrast() for se's
## Tables of means
## Grand mean
##          
## 121.9273 
## 
##  pH 
##       7.4   7.6   7.8     8
##     129.1 124.8 114.1 118.4
## rep  31.0  32.0  25.0  34.0
## 
##  temp 
##         6     9    12
##     113.1 123.8 128.5
## rep  35.0  56.0  31.0
## 
##  pH:temp 
##      temp
## pH    6      9      12    
##   7.4 120.31 134.08 129.04
##   rep   9.00  16.00   6.00
##   7.6 117.79 122.66 135.27
##   rep   9.00  14.00   9.00
##   7.8 104.65 115.97 122.33
##   rep   8.00  10.00   7.00
##   8   108.42 120.66 124.28
##   rep   9.00  16.00   9.00
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mean ~ pH * temp, data = mercenaria)
## 
## $pH
##               diff        lwr        upr     p adj
## 7.6-7.4  -4.265891  -9.487724  0.9559414 0.1497220
## 7.8-7.4 -14.979716 -20.549671 -9.4097614 0.0000000
## 8-7.4   -10.724193 -15.869897 -5.5784877 0.0000020
## 7.8-7.6 -10.713825 -16.244790 -5.1828597 0.0000103
## 8-7.6    -6.458301 -11.561777 -1.3548261 0.0070115
## 8-7.8     4.255524  -1.203627  9.7146737 0.1820655
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mean ~ pH * temp, data = mercenaria)
## 
## $temp
##           diff        lwr       upr     p adj
## 9-6  10.659263  6.5930486 14.725478 0.0000000
## 12-6 15.380816 10.7265116 20.035120 0.0000000
## 12-9  4.721552  0.4969832  8.946121 0.0244987

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group  11   0.308  0.984
##       361

## Anova Table (Type III tests)
## 
## Response: mean
##             Sum Sq  Df   F value    Pr(>F)    
## (Intercept) 130269   1 2064.8306 < 2.2e-16 ***
## pH            1434   3    7.5782 0.0001182 ***
## temp          1092   2    8.6540 0.0003232 ***
## pH:temp        842   6    2.2253 0.0458596 *  
## Residuals     6940 110                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Shapiro-Wilk normality test
## 
## data:  aov_residuals
## W = 0.9885, p-value = 0.3978

Plot

This plot illustrates the trend in both pH and temperature for Mercenaria shell coloration.

##    temp  pH  N     mean        sd       se       ci
## 1     6 7.4  9 120.3091  5.836794 1.945598 3.617933
## 2     6 7.6  9 117.7914  5.050849 1.683616 3.130766
## 3     6 7.8  8 104.6461  4.402510 1.556522 2.948954
## 4     6   8  9 108.4222  7.452191 2.484064 4.619235
## 5     9 7.4 16 134.0775  9.955011 2.488753 4.362909
## 6     9 7.6 14 122.6619  8.737814 2.335279 4.135624
## 7     9 7.8 10 115.9668  5.505253 1.740914 3.191292
## 8     9   8 16 120.6646 10.680070 2.670018 4.680675
## 9    12 7.4  6 129.0363  3.197103 1.305212 2.630065
## 10   12 7.6  9 135.2718  6.861502 2.287167 4.253098
## 11   12 7.8  7 122.3261  8.823586 3.335002 6.480510
## 12   12   8  9 124.2776  8.105061 2.701687 5.023917

*** ## mya ### check structure Appears pretty normal, one outlier.

## Classes 'data.table' and 'data.frame':   150 obs. of  18 variables:
##  $ Tank           : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ File path      : chr  "C:\\Users\\thatc\\OneDrive - Iowa State University\\Tank Experiment_Boron Isotopes Paper\\Mya arenaria pics\\jpg\\r37.jpeg" "C:\\Users\\thatc\\OneDrive - Iowa State University\\Tank Experiment_Boron Isotopes Paper\\Mya arenaria pics\\jpg\\r38.jpeg" "C:\\Users\\thatc\\OneDrive - Iowa State University\\Tank Experiment_Boron Isotopes Paper\\Mya arenaria pics\\jpg\\r39.jpeg" "C:\\Users\\thatc\\OneDrive - Iowa State University\\Tank Experiment_Boron Isotopes Paper\\Mya arenaria pics\\jpg\\r40.jpeg" ...
##  $ file name      : chr  "r37" "r38" "r39" "r40" ...
##  $ shell          : chr  "r37" "r38" "r39" "r40" ...
##  $ N              : num  417483 267263 707170 587852 661697 ...
##  $ mean           : num  158 149 140 140 130 ...
##  $ stdev          : num  39.6 40.6 41.2 46.1 45.9 ...
##  $ mode           : num  174 181 117 171 127 96 91 180 172 170 ...
##  $ pH             : Factor w/ 4 levels "7.4","7.6","7.8",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ temp           : Factor w/ 3 levels "6","9","12": 3 3 3 3 3 3 3 3 3 3 ...
##  $ dead?          : chr  NA NA NA NA ...
##  $ species        : chr  "mya" "mya" "mya" "mya" ...
##  $ temp.average   : num  12.1 12.1 12.1 12.1 12.1 ...
##  $ temp.stdev     : num  0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 ...
##  $ pH.average     : num  7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 ...
##  $ pH.stdev       : num  0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 ...
##  $ pH.normalized  : num  0.21 0.21 0.21 0.21 0.21 0.21 0.21 0.21 0.21 0.21 ...
##  $ temp.normalized: num  5.92 5.92 5.92 5.92 5.92 5.92 5.92 5.92 5.92 5.92 ...
##  - attr(*, ".internal.selfref")=<externalptr>

## 
##  Shapiro-Wilk normality test
## 
## data:  mya$mean
## W = 0.98651, p-value = 0.1533

linear model

Best fit model included only pH and the random effect of temp (shown by change in intercepts) the estimate was -48.57 for pH.

## Backward reduced random-effect table:
## 
##            Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                   6 -545.54 1103.1                         
## (1 | Tank)          0    5 -554.14 1118.3 17.207  1  3.353e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                               Eliminated Sum Sq Mean Sq NumDF  DenDF F value
## temp.normalized:pH.normalized          1  132.2   132.2     1 11.919  1.6289
## temp.normalized                        2  166.9   166.9     1 12.903  2.0580
## pH.normalized                          0 3528.6  3528.6     1 13.935 43.5486
##                                  Pr(>F)    
## temp.normalized:pH.normalized    0.2262    
## temp.normalized                  0.1752    
## pH.normalized                 1.219e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## mean ~ pH.normalized + (1 | Tank)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean ~ pH.normalized + (1 | Tank)
##    Data: mya
## 
## REML criterion at convergence: 1100.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.90190 -0.59486  0.05656  0.76082  2.64882 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 36.27    6.023   
##  Residual             81.03    9.002   
## Number of obs: 150, groups:  Tank, 16
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)    149.361      3.055  13.872  48.887  < 2e-16 ***
## pH.normalized  -48.568      7.360  13.935  -6.599 1.22e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## pH.normalzd -0.836

## [1] 1108.496
  mean
Predictors Estimates CI p
(Intercept) 149.36 143.32 – 155.40 <0.001
pH normalized -48.57 -63.11 – -34.02 <0.001
Random Effects
σ2 81.03
Ï„00 Tank 36.27
ICC 0.31
N Tank 16
Observations 150
Marginal R2 / Conditional R2 0.513 / 0.663

ANOVA

ANOVA results show similar patterns to the linear model but with temp term on the edge (0.027) and interaction term very significant. Tukey showed no significant difference for 8-7.8 and for 12-6 and 12-9. Passes homogenous variances and normal residuals. Type III ANOVA showed significant pH and temperature terms but no siginificant interaction.

## Call:
##    aov(formula = mean ~ pH * temp, data = mya)
## 
## Terms:
##                        pH      temp   pH:temp Residuals
## Sum of Squares  20422.167   616.664  2217.976 11551.518
## Deg. of Freedom         3         2         6       138
## 
## Residual standard error: 9.149134
## Estimated effects may be unbalanced
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## pH            3  20422    6807  81.324  < 2e-16 ***
## temp          2    617     308   3.683 0.027639 *  
## pH:temp       6   2218     370   4.416 0.000408 ***
## Residuals   138  11552      84                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Design is unbalanced - use se.contrast() for se's
## Tables of means
## Grand mean
##          
## 132.6371 
## 
##  pH 
##       7.4 7.6   7.8     8
##     150.3 136 121.1 122.9
## rep  37.0  39  38.0  36.0
## 
##  temp 
##       6     9    12
##     134 133.5 128.9
## rep  40  77.0  33.0
## 
##  pH:temp 
##      temp
## pH    6      9      12    
##   7.4 150.52 150.84 148.48
##   rep  10.00  20.00   7.00
##   7.6 134.93 135.39 138.40
##   rep  10.00  19.00  10.00
##   7.8 125.84 125.76 106.13
##   rep   9.00  20.00   9.00
##   8   125.01 121.70 122.84
##   rep  11.00  18.00   7.00
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mean ~ pH * temp, data = mya)
## 
## $pH
##               diff        lwr        upr    p adj
## 7.6-7.4 -14.261505 -19.721963  -8.801046 0.000000
## 7.8-7.4 -29.173417 -34.668743 -23.678091 0.000000
## 8-7.4   -27.371255 -32.941378 -21.801131 0.000000
## 7.8-7.6 -14.911912 -20.335376  -9.488449 0.000000
## 8-7.6   -13.109750 -18.608989  -7.610511 0.000000
## 8-7.8     1.802162  -3.731699   7.336024 0.831887
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mean ~ pH * temp, data = mya)
## 
## $temp
##           diff        lwr         upr     p adj
## 9-6  -0.512583  -4.737445  3.71227911 0.9554883
## 12-6 -5.192278 -10.289921 -0.09463448 0.0448506
## 12-9 -4.679695  -9.189824 -0.16956547 0.0400484

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group  11   0.308  0.984
##       361

## Anova Table (Type III tests)
## 
## Response: mean
##             Sum Sq  Df   F value    Pr(>F)    
## (Intercept) 226568   1 2706.6955 < 2.2e-16 ***
## pH            4228   3   16.8354 2.242e-09 ***
## temp            30   2    0.1766 0.8383140    
## pH:temp       2218   6    4.4162 0.0004079 ***
## Residuals    11552 138                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Shapiro-Wilk normality test
## 
## data:  aov_residuals
## W = 0.98342, p-value = 0.06818

Plot

This plot depicts the strong relationship between pH and shell color in Mya, seemingly regardless of temperature.

##    temp  pH  N     mean        sd       se       ci
## 1     6 7.4 10 150.5219  5.412784 1.711673 3.137689
## 2     6 7.6 10 134.9325  5.662082 1.790508 3.282203
## 3     6 7.8  9 125.8381  6.966905 2.322302 4.318432
## 4     6   8 11 125.0058  5.317460 1.603274 2.905873
## 5     9 7.4 20 150.8367 14.428803 3.226378 5.578837
## 6     9 7.6 19 135.3885  9.583922 2.198703 3.812690
## 7     9 7.8 20 125.7649  7.166237 1.602419 2.770796
## 8     9   8 18 121.7049  7.966948 1.877828 3.266681
## 9    12 7.4  7 148.4751  6.220035 2.350952 4.568324
## 10   12 7.6 10 138.3984 13.599189 4.300441 7.883194
## 11   12 7.8  9 106.1281  8.849082 2.949694 5.485098
## 12   12   8  7 122.8366  4.856334 1.835522 3.566750

scallop

check structure

Seemed weirder than other groups, skewed right with four outliers. When these were removed (r36 r76 r80 r94) the distribution was normal. The rest was done with this subset.

str(scallop)
## Classes 'data.table' and 'data.frame':   35 obs. of  18 variables:
##  $ Tank           : Factor w/ 16 levels "H1","H10","H11",..: 2 2 2 3 3 3 4 4 5 5 ...
##  $ File path      : chr  "C:\\Users\\thatc\\OneDrive - Iowa State University\\Desktop\\jpg\\scallops\\r55.jpg" NA NA "C:\\Users\\thatc\\OneDrive - Iowa State University\\Desktop\\jpg\\scallops\\r36.jpg" ...
##  $ file name      : chr  "r55.jpg" NA NA "r36.jpg" ...
##  $ shell          : chr  "r55" "r34" "r67" "r36" ...
##  $ N              : num  2733522 2087642 1239933 964342 1839238 ...
##  $ mean           : num  104.9 83.1 96.6 143.2 80.8 ...
##  $ stdev          : num  41.1 22.8 32.2 27.2 23.1 ...
##  $ mode           : num  117 77 91 143 74 73 101 78 95 100 ...
##  $ pH             : Factor w/ 4 levels "7.4","7.6","7.8",..: 4 4 4 4 4 4 4 4 1 1 ...
##  $ temp           : Factor w/ 3 levels "6","9","12": 2 2 2 1 1 1 2 2 3 3 ...
##  $ dead?          : chr  NA NA NA NA ...
##  $ species        : chr  "scallop" "scallop" "scallop" "scallop" ...
##  $ temp.average   : num  9.19 9.19 9.19 6.26 6.26 ...
##  $ temp.stdev     : num  0.63 0.63 0.63 0.78 0.78 0.78 0.33 0.33 0.75 0.75 ...
##  $ pH.average     : num  8 8 8 8.02 8.02 8.02 8 8 7.44 7.44 ...
##  $ pH.stdev       : num  0.07 0.07 0.07 0.06 0.06 0.06 0.06 0.06 0.07 0.07 ...
##  $ pH.normalized  : num  0.64 0.64 0.64 0.66 0.66 ...
##  $ temp.normalized: num  3.05 3.05 3.05 0.12 0.12 ...
##  - attr(*, ".internal.selfref")=<externalptr>
simple.eda.ggplot(scallop$mean) #skewed right, could remove outliers? or transform?

shapiro.test(scallop$mean) #again not normal
## 
##  Shapiro-Wilk normality test
## 
## data:  scallop$mean
## W = 0.87181, p-value = 0.000751
Q <- quantile(scallop$mean, probs=c(.25, .75), na.rm = FALSE)
iqr <- IQR(scallop$mean)
up <-  Q[2]+1.5*iqr # Upper Range  
low<- Q[1]-1.5*iqr # Lower Range
scallop.no.out<- subset(scallop, scallop$mean > (Q[1] - 1.5*iqr) & scallop$mean < (Q[2]+1.5*iqr))
dim(scallop)-dim(scallop.no.out) # 4 removed
## [1] 4 0
simple.eda.ggplot(scallop.no.out$mean) #slightly better

shapiro.test(scallop.no.out$mean) #better, so will continue without outliers, but should make note that they were 
## 
##  Shapiro-Wilk normality test
## 
## data:  scallop.no.out$mean
## W = 0.96715, p-value = 0.4444
#summary(comparedf(scallop, scallop.no.out, by = "shell")) #r36 r76 r80 r94

linear model

The only sigificant predictor was Intercept (Tank) and the best fit model only included this random term.

## Backward reduced random-effect table:
## 
##            Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)  
## <none>                   6 -110.27 232.53                       
## (1 | Tank)          0    5 -111.66 233.32 2.7865  1    0.09506 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                               Eliminated Sum Sq Mean Sq NumDF  DenDF F value
## temp.normalized:pH.normalized          1   4.21    4.21     1 12.351  0.0424
## temp.normalized                        2  24.75   24.75     1 11.725  0.2519
## pH.normalized                          3 371.01  371.01     1 15.608  3.8232
##                                Pr(>F)  
## temp.normalized:pH.normalized 0.84026  
## temp.normalized               0.62501  
## pH.normalized                 0.06869 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## mean ~ (1 | Tank)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean ~ (1 | Tank)
##    Data: scallop.no.out
## 
## REML criterion at convergence: 239.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.37454 -0.52271  0.00075  0.49317  2.22303 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 78.86    8.881   
##  Residual             99.05    9.952   
## Number of obs: 31, groups:  Tank, 15
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)   93.227      2.969 14.067    31.4 1.98e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean ~ temp.normalized * pH.normalized + (1 | Tank)
##    Data: scallop.no.out
## 
## REML criterion at convergence: 220.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.31627 -0.54105 -0.01677  0.37491  1.94213 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 71.84    8.476   
##  Residual             99.40    9.970   
## Number of obs: 31, groups:  Tank, 15
## 
## Fixed effects:
##                               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)                     98.776      9.637  13.695  10.249 8.42e-08 ***
## temp.normalized                  1.232      2.895  11.107   0.425    0.679    
## pH.normalized                  -20.174     21.917  13.495  -0.920    0.373    
## temp.normalized:pH.normalized   -1.391      6.757  12.351  -0.206    0.840    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.nr pH.nrm
## temp.nrmlzd -0.812              
## pH.normalzd -0.855  0.698       
## tmp.nrml:H.  0.679 -0.850 -0.805

## [1] 245.5285
  mean
Predictors Estimates CI p
(Intercept) 93.23 87.15 – 99.31 <0.001
Random Effects
σ2 99.05
Ï„00 Tank 78.86
ICC 0.44
N Tank 15
Observations 31
Marginal R2 / Conditional R2 0.000 / 0.443

ANOVA

The ANOVA results show a similar to model but with pH term significant (0.00467). As in other cases the design is unbalanced but here we have 0 scallops t 7.6 and 12 so we cant do an interaction term I believe (think you need a minimum of two reps in each interaction). I therefore did an ANOVA that didnt include interaction (mean~pH+temp). This also indicated significant pH term in the normal ANOVA (0.00397) and type III (0.007861). Tukey showed no significant difference between 8 and all other pHs and between 7.8-7.6. Residuals appeared normal.

## Call:
##    aov(formula = mean ~ pH * temp, data = scallop.no.out)
## 
## Terms:
##                        pH      temp   pH:temp Residuals
## Sum of Squares  2139.2838   16.7010  695.6561 2413.8567
## Deg. of Freedom         3         2         5        20
## 
## Residual standard error: 10.98603
## 1 out of 12 effects not estimable
## Estimated effects may be unbalanced
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## pH           3 2139.3   713.1   5.908 0.00467 **
## temp         2   16.7     8.4   0.069 0.93337   
## pH:temp      5  695.7   139.1   1.153 0.36609   
## Residuals   20 2413.9   120.7                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Design is unbalanced - use se.contrast() for se's
## Tables of means
## Grand mean
##          
## 93.01642 
## 
##  pH 
##       7.4   7.6   7.8    8
##     109.1 89.02 85.93 92.9
## rep   6.0  8.00  9.00  8.0
## 
##  temp 
##         6    9    12
##     91.99 93.6 93.26
## rep 10.00 15.0  6.00
## 
##  pH:temp 
##      temp
## pH    6      9      12    
##   7.4 110.10 101.46 113.93
##   rep   1.00   2.00   3.00
##   7.6  91.25  86.79       
##   rep   4.00   4.00   0.00
##   7.8  86.88  88.53  79.29
##   rep   3.00   4.00   2.00
##   8    81.42  97.43  93.24
##   rep   2.00   5.00   1.00
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mean ~ pH * temp, data = scallop.no.out)
## 
## $pH
##               diff        lwr        upr     p adj
## 7.6-7.4 -20.118333 -36.724796 -3.5118705 0.0141694
## 7.8-7.4 -23.208056 -39.414304 -7.0018073 0.0035487
## 8-7.4   -16.231458 -32.837921  0.3750045 0.0568515
## 7.8-7.6  -3.089722 -18.031145 11.8517004 0.9372936
## 8-7.6     3.886875 -11.487722 19.2614721 0.8928948
## 8-7.8     6.976597  -7.964825 21.9180199 0.5693393
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mean ~ pH * temp, data = scallop.no.out)
## 
## $temp
##            diff        lwr      upr     p adj
## 9-6   1.6070815  -9.739957 12.95412 0.9319116
## 12-6  1.2656866 -13.087308 15.61868 0.9729709
## 12-9 -0.3413949 -13.767392 13.08460 0.9977211

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group  11   0.308  0.984
##       361

## 
##  Shapiro-Wilk normality test
## 
## data:  aov_residuals
## W = 0.98236, p-value = 0.8746
## Call:
##    aov(formula = mean ~ pH + temp, data = scallop.no.out)
## 
## Terms:
##                       pH     temp Residuals
## Sum of Squares  2139.284   16.701  3109.513
## Deg. of Freedom        3        2        25
## 
## Residual standard error: 11.1526
## Estimated effects may be unbalanced
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## pH           3 2139.3   713.1   5.733 0.00397 **
## temp         2   16.7     8.4   0.067 0.93524   
## Residuals   25 3109.5   124.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Design is unbalanced - use se.contrast() for se's
## Tables of means
## Grand mean
##          
## 93.01642 
## 
##  pH 
##       7.4   7.6   7.8    8
##     109.1 89.02 85.93 92.9
## rep   6.0  8.00  9.00  8.0
## 
##  temp 
##         6    9    12
##     91.99 93.6 93.26
## rep 10.00 15.0  6.00
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mean ~ pH + temp, data = scallop.no.out)
## 
## $pH
##               diff        lwr        upr     p adj
## 7.6-7.4 -20.118333 -36.685715 -3.5509522 0.0131236
## 7.8-7.4 -23.208056 -39.376164 -7.0399471 0.0029780
## 8-7.4   -16.231458 -32.798840  0.3359228 0.0563324
## 7.8-7.6  -3.089722 -17.995982 11.8165373 0.9400402
## 8-7.6     3.886875 -11.451540 19.2252895 0.8972455
## 8-7.8     6.976597  -7.929662 21.8828567 0.5791539
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mean ~ pH + temp, data = scallop.no.out)
## 
## $temp
##            diff        lwr      upr     p adj
## 9-6   1.6070815  -9.733744 12.94791 0.9338108
## 12-6  1.2656866 -13.079449 15.61082 0.9737536
## 12-9 -0.3413949 -13.760040 13.07725 0.9977886

## Anova Table (Type III tests)
## 
## Response: mean
##             Sum Sq Df  F value   Pr(>F)    
## (Intercept)  38329  1 308.1605 1.43e-15 ***
## pH            1845  3   4.9440 0.007861 ** 
## temp            17  2   0.0671 0.935235    
## Residuals     3110 25                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Shapiro-Wilk normality test
## 
## data:  aov_residuals
## W = 0.95551, p-value = 0.2211

Plot

These plots show that there may be something going on with pH but we have poor coverage and large error due to the small sample size. This is also only data with the four outliers removed.

##    temp  pH N      mean         sd       se        ci
## 1     6 7.4 1 110.09600         NA       NA       NaN
## 2     6 7.6 4  91.24725 12.0427456 6.021373 14.170479
## 3     6 7.8 3  86.88233 10.9553238 6.325059 18.469082
## 4     6   8 2  81.41500  0.8061017 0.570000  3.598838
## 5     9 7.4 2 101.45500  2.4720453 1.748000 11.036438
## 6     9 7.6 4  86.78575 10.8310821 5.415541 12.744736
## 7     9 7.8 4  88.52925 13.9115648 6.955782 16.369484
## 8     9   8 5  97.43140  8.7275949 3.903099  8.320809
## 9    12 7.4 3 113.93433 14.5771304 8.416110 24.574920
## 10   12 7.8 2  79.28850  8.3530524 5.906500 37.292173
## 11   12   8 1  93.24000         NA       NA        NA

***

juv arctica

check structure

This all looks pretty normal to me.

## Classes 'data.table' and 'data.frame':   66 obs. of  18 variables:
##  $ Tank           : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ File path      : chr  "C:\\Users\\thatc\\OneDrive - Iowa State University\\Desktop\\jpg\\juv Arctica\\g29.jpeg" NA NA NA ...
##  $ file name      : chr  "g29.jpeg" NA NA NA ...
##  $ shell          : chr  "g29" "g30" "g31" "g32" ...
##  $ N              : num  381123 339242 764948 602176 620811 ...
##  $ mean           : num  52.5 41 58.5 52.7 51.8 ...
##  $ stdev          : num  9.96 17.64 19.74 16.74 23.15 ...
##  $ mode           : num  47 28 46 51 46 50 68 47 59 51 ...
##  $ pH             : Factor w/ 4 levels "7.4","7.6","7.8",..: 2 2 2 2 4 4 4 4 4 4 ...
##  $ temp           : Factor w/ 3 levels "6","9","12": 3 3 3 3 2 2 2 2 1 1 ...
##  $ dead?          : chr  NA NA NA NA ...
##  $ species        : chr  "juv arctica" "juv arctica" "juv arctica" "juv arctica" ...
##  $ temp.average   : num  12.06 12.06 12.06 12.06 9.19 ...
##  $ temp.stdev     : num  0.42 0.42 0.42 0.42 0.63 0.63 0.63 0.63 0.78 0.78 ...
##  $ pH.average     : num  7.57 7.57 7.57 7.57 8 8 8 8 8.02 8.02 ...
##  $ pH.stdev       : num  0.04 0.04 0.04 0.04 0.07 0.07 0.07 0.07 0.06 0.06 ...
##  $ pH.normalized  : num  0.21 0.21 0.21 0.21 0.64 ...
##  $ temp.normalized: num  5.92 5.92 5.92 5.92 3.05 ...
##  - attr(*, ".internal.selfref")=<externalptr>

## 
##  Shapiro-Wilk normality test
## 
## data:  juv.arctica$mean
## W = 0.9902, p-value = 0.8847

linear model

The best fit model only includes the intercept, but unlike in other cases this doesnt seem to mean ‘Tank’/the random term as there is only one coefficent/estamite. I need to read more on what this means.

## Backward reduced random-effect table:
## 
##            Eliminated npar  logLik    AIC     LRT Df Pr(>Chisq)
## <none>                   6 -207.27 426.54                      
## (1 | Tank)          1    5 -207.43 424.86 0.32521  1     0.5685
## 
## Backward reduced fixed-effect table:
##                               Eliminated Df Sum of Sq    RSS    AIC F value
## temp.normalized:pH.normalized          1  1     0.040 2353.8 241.89  0.0011
## pH.normalized                          2  1    74.861 2428.7 241.96  2.0036
## temp.normalized                        3  1   112.739 2541.4 242.96  2.9708
##                                Pr(>F)  
## temp.normalized:pH.normalized 0.97417  
## pH.normalized                 0.16185  
## temp.normalized               0.08961 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## mean ~ 1
## 
## Call:
## lm(formula = mean ~ 1, data = juv.arctica)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.9847  -4.1685  -0.2952   4.0418  15.7703 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  59.9697     0.7697   77.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.253 on 65 degrees of freedom

## hat values (leverages) are all = 0.01515152
##  and there are no factor predictors; no plot no. 5

## [1] 432.2551
  mean
Predictors Estimates CI p
(Intercept) 59.97 58.43 – 61.51 <0.001
Observations 66
R2 / R2 adjusted 0.000 / 0.000

ANOVA

The summary for a normal ANOVA shows similar results to model but with pH term on the edge (0.06). Here we see that as with the rest our design is unbalanced… but all groups are covered. No pH or temperature groups are significantly different. No evidence to suggest that the variance across groups is statistically significantly different. Type III ANOVA is apparently more accurate in unbalanced designs, which we have here. Therefore these results are likely more accurate than the raw aov output and explain why this could be more close to the best fit model, the intercept is the only significant term. I trust this more than the straight up anova as it accounts for the unbalanced design. The intercept is the estimate of the dependent variable when all the independent variables are 0. So it just means mean is significantly different from 0.

## Call:
##    aov(formula = mean ~ pH * temp, data = juv.arctica)
## 
## Terms:
##                        pH      temp   pH:temp Residuals
## Sum of Squares   255.3089  139.0573  365.4394 1781.6465
## Deg. of Freedom         3         2         6        54
## 
## Residual standard error: 5.743993
## Estimated effects may be unbalanced
##             Df Sum Sq Mean Sq F value Pr(>F)  
## pH           3  255.3   85.10   2.579  0.063 .
## temp         2  139.1   69.53   2.107  0.131  
## pH:temp      6  365.4   60.91   1.846  0.107  
## Residuals   54 1781.6   32.99                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Design is unbalanced - use se.contrast() for se's
## Tables of means
## Grand mean
##          
## 59.96971 
## 
##  pH 
##       7.4  7.6   7.8     8
##     63.41 58.6 58.25 59.86
## rep 15.00 16.0 16.00 19.00
## 
##  temp 
##         6     9   12
##     60.63 60.81 57.3
## rep 16.00 35.00 15.0
## 
##  pH:temp 
##      temp
## pH    6     9     12   
##   7.4 64.93 61.37 66.82
##   rep  4.00  8.00  3.00
##   7.6 58.22 62.51 51.15
##   rep  4.00  8.00  4.00
##   7.8 58.21 59.77 55.24
##   rep  4.00  8.00  4.00
##   8   61.40 60.07 57.74
##   rep  4.00 11.00  4.00
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mean ~ pH * temp, data = juv.arctica)
## 
## $pH
##               diff        lwr       upr     p adj
## 7.6-7.4 -4.8117917 -10.284199 0.6606158 0.1036631
## 7.8-7.4 -5.1621042 -10.634512 0.3103033 0.0711566
## 8-7.4   -3.5537193  -8.812926 1.7054873 0.2886749
## 7.8-7.6 -0.3503125  -5.733732 5.0331068 0.9981528
## 8-7.6    1.2580724  -3.908475 6.4246200 0.9166657
## 8-7.8    1.6083849  -3.558163 6.7749325 0.8422669
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mean ~ pH * temp, data = juv.arctica)
## 
## $temp
##            diff       lwr       upr     p adj
## 9-6   0.1800824 -3.997443 4.3576079 0.9940681
## 12-6 -3.3304606 -8.305577 1.6446560 0.2488274
## 12-9 -3.5105430 -7.782564 0.7614781 0.1268255

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group  11   0.308  0.984
##       361

## Anova Table (Type III tests)
## 
## Response: mean
##              Sum Sq Df  F value Pr(>F)    
## (Intercept) 16862.2  1 511.0769 <2e-16 ***
## pH            122.7  3   1.2401 0.3042    
## temp           77.3  2   1.1721 0.3175    
## pH:temp       365.4  6   1.8460 0.1074    
## Residuals    1781.6 54                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Shapiro-Wilk normality test
## 
## data:  aov_residuals
## W = 0.98425, p-value = 0.5675

Plot

These plots show how there really doesn’t seem to be much relationship between the color of juvenile Arctica shells and the pH or temperature treatments.

##    temp  pH  N     mean       sd       se        ci
## 1     6 7.4  4 64.92725 8.478716 4.239358  9.976750
## 2     6 7.6  4 58.22325 4.905617 2.452809  5.772350
## 3     6 7.8  4 58.21125 8.986877 4.493438 10.574694
## 4     6   8  4 61.39900 4.614352 2.307176  5.429624
## 5     9 7.4  8 61.37288 3.151784 1.114324  2.111174
## 6     9 7.6  8 62.50863 6.428547 2.272835  4.306064
## 7     9 7.8  8 59.77137 4.823539 1.705379  3.230974
## 8     9   8 11 60.06555 6.548514 1.974451  3.578616
## 9    12 7.4  3 66.82267 3.444560 1.988717  5.807026
## 10   12 7.6  4 51.15500 7.338287 3.669144  8.634829
## 11   12 7.8  4 55.24025 3.043270 1.521635  3.580960
## 12   12   8  4 57.74125 2.830496 1.415248  3.330593