Goal of this markdown: Using a AIC model ranking to analyze my evironmental and field data. Looking at models with environmental data fitlered between field surveys, one month (30 days) before field surveys, and two months (60 days) before field surveys.

I am addressing research question #3:
Is variation in salinity, air temperature, or pH associated with changes in abundances or reproductive effort of F. distichus populations in SFE? –> compare models and see which ones have the most evidence (based on fit of model and evidence)
- assumptions: environmental variables matter (lit evidence and prior knowledge)
- Environmental variables have been filtered for tides

AIC link:
https://www.statology.org/aic-in-r/

Models
I’ll try the same model list with a one month lag and then a two month lag. Focusing on 1) daily salinity range and low salinity events (number of days with a daily maximum less than 10), 2) high air temperature events (number of days with a maximum daily temperature over 26C) , 3) daily pH range and high pH events, and 4) adult cover (large thalli density). Previously I tested both raw and square root transformed for juvenile abundance and the untransformed data fit better and had a higher R^2.

  1. Juv ~ adult cover + daily salinity range + high air temp event + daily pH range + adult cover: daily air temp event
  2. Juv ~ adult cover + low salinity event + high air temp event + high pH event + adult cover : high air temp event
  3. Juv ~ adult cover + daily salinity range + high air temp event + adult cover : high air temp event 4. Juv ~ adult cover + low salinity event + high air temp event + adult cover : high air temp event
  4. Juv ~ adult cover + high air temp event + adult cover : high air temp event
  5. Juv ~ adult cover + daily salinity range + high air temp event
  6. Juv ~ low salinity event + high pH event + high air temp event
  7. Juv ~ daily salinity range + daily pH range + high air temp event
  8. Juv ~ adult cover + high air temp event
  9. Juv ~ daily salinity range + daily pH range
  10. Juv ~ low salinity event + high pH event
  11. Juv ~ low salinity event + high air temp event

outliers 3: highest mean juvenile denisty 4: third highest mean juvenile density 32: fourth highest mean juvenile density 37: CODING ERROR with max daily air temp greater than 26 –> result of a coding error I found, fixed, and checked for at all sites. No long an outlier

Next steps - try median instead of mean juvenile abundance

Set up
Loading packages, reading in and prepping data

rm(list=ls())

library(tidyverse)
library(ggpubr)
library(scales)
library(chron)
library(plotly)
library(taRifx)
library(aweek)
library(easypackages)
library(renv)
library(here)
library(ggthemes)
library(gridExtra)
library(patchwork)
library(tidyquant)
library(recipes) 
library(cranlogs)
library(knitr)
library(openair)
library(xts)
library(dplyr)

#AIC
library(AICcmodavg)
library(plm) #for random factors
library(lme4)
library(car)

#Table 
library(sjPlot) 
library(sjmisc)
library(sjlabelled)
library(insight)
library(tibble)

Envirnomental data fitlered 30 days before field survey

Envirnomental data 30 days leading up to survey day (excluding day of survey) and adult thalli density from the field survey one month before. Juvenile density is the mean.

Read in data

#read in data
one.month<-read.csv(
      "https://raw.githubusercontent.com/Cmwegener/thesis/master/data/environment_field/envi.field.all.one.month.lag.csv",
    header = TRUE
  )
one.month$date<-as.Date(one.month$date, format = c("%Y-%m-%d"))

AIC

#models
m1 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + daily.sal.range + max.daily.air.gt26 + daily.ph.range + lag.no.large.fuc.q:max.daily.air.gt26, data=one.month)

m2 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.sal.lt10 + max.daily.air.gt26 + max.daily.ph.gt8 + lag.no.large.fuc.q:max.daily.air.gt26, data=one.month)

m3 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + daily.sal.range + max.daily.air.gt26 + lag.no.large.fuc.q:max.daily.air.gt26, data=one.month)

m4 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.sal.lt10 + max.daily.air.gt26 + lag.no.large.fuc.q:max.daily.air.gt26, data=one.month)

m5 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.air.gt26 + lag.no.large.fuc.q:max.daily.air.gt26, data=one.month)

m6 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + daily.sal.range + max.daily.air.gt26, data=one.month)

m7 <- lm(no.small.fuc.q ~ max.daily.sal.lt10 + max.daily.ph.gt8 + max.daily.air.gt26, data=one.month)

m8 <- lm(no.small.fuc.q ~ daily.sal.range + daily.ph.range + max.daily.air.gt26, data=one.month)

m9 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.air.gt26, data=one.month)

m10 <- lm(no.small.fuc.q ~ daily.sal.range + daily.ph.range, data=one.month)

m11 <- lm(no.small.fuc.q ~ max.daily.sal.lt10 + max.daily.ph.gt8, data=one.month)

m12 <- lm(no.small.fuc.q ~ max.daily.sal.lt10 + max.daily.air.gt26, data=one.month)

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

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

#calculate AIC of each model
aictab(cand.set = models, modnames = mod.names)
## 
## Model selection based on AICc:
## 
##     K   AICc Delta_AICc AICcWt Cum.Wt      LL
## m1  7 395.48       0.00   0.96   0.96 -188.93
## m2  7 401.78       6.30   0.04   1.00 -192.09
## m8  5 424.88      29.40   0.00   1.00 -206.61
## m10 4 425.32      29.84   0.00   1.00 -208.12
## m11 4 430.60      35.12   0.00   1.00 -210.76
## m7  5 433.06      37.58   0.00   1.00 -210.70
## m6  5 515.41     119.92   0.00   1.00 -252.06
## m3  6 517.94     122.46   0.00   1.00 -252.06
## m9  4 521.93     126.45   0.00   1.00 -256.55
## m5  5 524.24     128.75   0.00   1.00 -256.48
## m4  6 526.78     131.30   0.00   1.00 -256.48
## m12 4 569.63     174.15   0.00   1.00 -280.43

Top ranked
m1 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + daily.sal.range + max.daily.air.gt26 + daily.ph.range + lag.no.large.fuc.q:max.daily.air.gt26, data=one.month)
Fitting and analyzing the top ranked model.

m1 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + daily.sal.range + max.daily.air.gt26 + daily.ph.range + lag.no.large.fuc.q:max.daily.air.gt26, data=one.month)
plot(m1) 

m1
## 
## Call:
## lm(formula = no.small.fuc.q ~ lag.no.large.fuc.q + daily.sal.range + 
##     max.daily.air.gt26 + daily.ph.range + lag.no.large.fuc.q:max.daily.air.gt26, 
##     data = one.month)
## 
## Coefficients:
##                           (Intercept)                     lag.no.large.fuc.q  
##                             58.753199                               1.000281  
##                       daily.sal.range                     max.daily.air.gt26  
##                             -4.504252                              -0.790638  
##                        daily.ph.range  lag.no.large.fuc.q:max.daily.air.gt26  
##                             -1.891923                               0.002399
summary(m1)
## 
## Call:
## lm(formula = no.small.fuc.q ~ lag.no.large.fuc.q + daily.sal.range + 
##     max.daily.air.gt26 + daily.ph.range + lag.no.large.fuc.q:max.daily.air.gt26, 
##     data = one.month)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.061 -24.344  -7.122  12.649 104.767 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                           58.753199  20.286896   2.896  0.00666 **
## lag.no.large.fuc.q                     1.000281   1.342969   0.745  0.46165   
## daily.sal.range                       -4.504252   1.774704  -2.538  0.01605 * 
## max.daily.air.gt26                    -0.790638   0.882614  -0.896  0.37685   
## daily.ph.range                        -1.891923   3.429397  -0.552  0.58489   
## lag.no.large.fuc.q:max.daily.air.gt26  0.002399   0.094643   0.025  0.97993   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.42 on 33 degrees of freedom
##   (22 observations deleted due to missingness)
## Multiple R-squared:  0.2047, Adjusted R-squared:  0.0842 
## F-statistic: 1.699 on 5 and 33 DF,  p-value: 0.1624
anova(m1)
## Analysis of Variance Table
## 
## Response: no.small.fuc.q
##                                       Df Sum Sq Mean Sq F value  Pr(>F)  
## lag.no.large.fuc.q                     1   1580  1580.1  1.4150 0.24271  
## daily.sal.range                        1   5706  5706.1  5.1099 0.03051 *
## max.daily.air.gt26                     1   1857  1856.8  1.6628 0.20620  
## daily.ph.range                         1    341   340.9  0.3053 0.58429  
## lag.no.large.fuc.q:max.daily.air.gt26  1      1     0.7  0.0006 0.97993  
## Residuals                             33  36850  1116.7                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(m1)
##                                            2.5 %      97.5 %
## (Intercept)                           17.4791983 100.0271992
## lag.no.large.fuc.q                    -1.7320093   3.7325716
## daily.sal.range                       -8.1149154  -0.8935885
## max.daily.air.gt26                    -2.5863298   1.0050534
## daily.ph.range                        -8.8690826   5.0852374
## lag.no.large.fuc.q:max.daily.air.gt26 -0.1901534   0.1949505

Table

tab_model(m1)
  no.small.fuc.q
Predictors Estimates CI p
(Intercept) 58.75 17.48 – 100.03 0.007
lag.no.large.fuc.q 1.00 -1.73 – 3.73 0.462
daily.sal.range -4.50 -8.11 – -0.89 0.016
max.daily.air.gt26 -0.79 -2.59 – 1.01 0.377
daily.ph.range -1.89 -8.87 – 5.09 0.585
lag.no.large.fuc.q *
max.daily.air.gt26 		0.00 -0.19 – 0.19 0.980
Observations 39
R2 / R2 adjusted 0.205 / 0.084

Second ranked
m2 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.sal.lt10 + max.daily.air.gt26 + max.daily.ph.gt8 + lag.no.large.fuc.q:max.daily.air.gt26, data=one.month) Fitting and analyzing the top ranked model.

m2 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.sal.lt10 + max.daily.air.gt26 + max.daily.ph.gt8 + lag.no.large.fuc.q:max.daily.air.gt26, data=one.month)
plot(m2) 

m2
## 
## Call:
## lm(formula = no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.sal.lt10 + 
##     max.daily.air.gt26 + max.daily.ph.gt8 + lag.no.large.fuc.q:max.daily.air.gt26, 
##     data = one.month)
## 
## Coefficients:
##                           (Intercept)                     lag.no.large.fuc.q  
##                             15.677698                               1.537695  
##                    max.daily.sal.lt10                     max.daily.air.gt26  
##                             -0.384003                               0.051113  
##                      max.daily.ph.gt8  lag.no.large.fuc.q:max.daily.air.gt26  
##                              0.676607                               0.006729
summary(m2)
## 
## Call:
## lm(formula = no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.sal.lt10 + 
##     max.daily.air.gt26 + max.daily.ph.gt8 + lag.no.large.fuc.q:max.daily.air.gt26, 
##     data = one.month)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -44.297 -18.246 -11.544   8.979 117.353 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)                           15.677698  13.880438   1.129    0.267
## lag.no.large.fuc.q                     1.537695   1.412855   1.088    0.284
## max.daily.sal.lt10                    -0.384003   2.730995  -0.141    0.889
## max.daily.air.gt26                     0.051113   0.890520   0.057    0.955
## max.daily.ph.gt8                       0.676607   0.899435   0.752    0.457
## lag.no.large.fuc.q:max.daily.air.gt26  0.006729   0.102540   0.066    0.948
## 
## Residual standard error: 36.23 on 33 degrees of freedom
##   (22 observations deleted due to missingness)
## Multiple R-squared:  0.06521,    Adjusted R-squared:  -0.07642 
## F-statistic: 0.4604 on 5 and 33 DF,  p-value: 0.8027
anova(m2)
## Analysis of Variance Table
## 
## Response: no.small.fuc.q
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## lag.no.large.fuc.q                     1   1580 1580.11  1.2039 0.2805
## max.daily.sal.lt10                     1     88   87.93  0.0670 0.7974
## max.daily.air.gt26                     1    611  610.56  0.4652 0.5000
## max.daily.ph.gt8                       1    737  737.30  0.5617 0.4589
## lag.no.large.fuc.q:max.daily.air.gt26  1      6    5.65  0.0043 0.9481
## Residuals                             33  43313 1312.53
confint(m2)
##                                             2.5 %     97.5 %
## (Intercept)                           -12.5622666 43.9176619
## lag.no.large.fuc.q                     -1.3367796  4.4121702
## max.daily.sal.lt10                     -5.9402538  5.1722480
## max.daily.air.gt26                     -1.7606630  1.8628891
## max.daily.ph.gt8                       -1.1533069  2.5065202
## lag.no.large.fuc.q:max.daily.air.gt26  -0.2018911  0.2153492

Table

tab_model(m2)
  no.small.fuc.q
Predictors Estimates CI p
(Intercept) 15.68 -12.56 – 43.92 0.267
lag.no.large.fuc.q 1.54 -1.34 – 4.41 0.284
max.daily.sal.lt10 -0.38 -5.94 – 5.17 0.889
max.daily.air.gt26 0.05 -1.76 – 1.86 0.955
max.daily.ph.gt8 0.68 -1.15 – 2.51 0.457
lag.no.large.fuc.q *
max.daily.air.gt26 		0.01 -0.20 – 0.22 0.948
Observations 39
R2 / R2 adjusted 0.065 / -0.076

Envirnomental data fitlered 60 days before field survey

Envirnomental data for 30 days starting 60 days from survey date (excludes the 30 days leading up to the survey date). Adult thalli density from two surveys/months before. Juvenile density is the mean.

Read in data

#read in data
two.month<-read.csv(
      "https://raw.githubusercontent.com/Cmwegener/thesis/master/data/environment_field/envi.field.all.two.month.lag.csv",
    header = TRUE
  )
two.month$date<-as.Date(two.month$date, format = c("%Y-%m-%d"))

AIC

#models
m1 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + daily.sal.range + max.daily.air.gt26 + daily.ph.range + lag.no.large.fuc.q:max.daily.air.gt26, data=two.month)

m2 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.sal.lt10 + max.daily.air.gt26 + max.daily.ph.gt8 + lag.no.large.fuc.q:max.daily.air.gt26, data=two.month)

m3 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + daily.sal.range + max.daily.air.gt26 + lag.no.large.fuc.q:max.daily.air.gt26, data=two.month)

m4 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.sal.lt10 + max.daily.air.gt26 + lag.no.large.fuc.q:max.daily.air.gt26, data=two.month)

m5 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.air.gt26 + lag.no.large.fuc.q:max.daily.air.gt26, data=two.month)

m6 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + daily.sal.range + max.daily.air.gt26, data=two.month)

m7 <- lm(no.small.fuc.q ~ max.daily.sal.lt10 + max.daily.ph.gt8 + max.daily.air.gt26, data=two.month)

m8 <- lm(no.small.fuc.q ~ daily.sal.range + daily.ph.range + max.daily.air.gt26, data=two.month)

m9 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.air.gt26, data=two.month)

m10 <- lm(no.small.fuc.q ~ daily.sal.range + daily.ph.range, data=two.month)

m11 <- lm(no.small.fuc.q ~ max.daily.sal.lt10 + max.daily.ph.gt8, data=two.month)

m12 <- lm(no.small.fuc.q ~ max.daily.sal.lt10 + max.daily.air.gt26, data=two.month)

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

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

#calculate AIC of each model
aictab(cand.set = models, modnames = mod.names)
## 
## Model selection based on AICc:
## 
##     K   AICc Delta_AICc AICcWt Cum.Wt      LL
## m1  7 347.63       0.00   0.69   0.69 -164.82
## m2  7 349.20       1.57   0.31   1.00 -165.60
## m10 4 428.06      80.42   0.00   1.00 -209.49
## m11 4 429.37      81.73   0.00   1.00 -210.14
## m8  5 429.76      82.12   0.00   1.00 -209.05
## m7  5 431.93      84.29   0.00   1.00 -210.13
## m9  4 449.02     101.39   0.00   1.00 -220.06
## m6  5 451.25     103.61   0.00   1.00 -219.93
## m5  5 451.42     103.78   0.00   1.00 -220.01
## m3  6 453.76     106.13   0.00   1.00 -219.88
## m4  6 453.91     106.27   0.00   1.00 -219.95
## m12 4 567.12     219.48   0.00   1.00 -279.17

The top two models are the same. The main difference is that the AICcWt for the second model is much higher than the previous method (30 days leading to survey date)

Top ranked
m1 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + daily.sal.range + max.daily.air.gt26 + daily.ph.range + lag.no.large.fuc.q:max.daily.air.gt26, data=two.month)
Fitting and analyzing the top ranked model.

m1 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + daily.sal.range + max.daily.air.gt26 + daily.ph.range + lag.no.large.fuc.q:max.daily.air.gt26, data=two.month)
plot(m1) 

m1
## 
## Call:
## lm(formula = no.small.fuc.q ~ lag.no.large.fuc.q + daily.sal.range + 
##     max.daily.air.gt26 + daily.ph.range + lag.no.large.fuc.q:max.daily.air.gt26, 
##     data = two.month)
## 
## Coefficients:
##                           (Intercept)                     lag.no.large.fuc.q  
##                              43.36859                               -0.56716  
##                       daily.sal.range                     max.daily.air.gt26  
##                              -1.98120                               -0.39524  
##                        daily.ph.range  lag.no.large.fuc.q:max.daily.air.gt26  
##                              -0.97481                                0.03345
summary(m1)
## 
## Call:
## lm(formula = no.small.fuc.q ~ lag.no.large.fuc.q + daily.sal.range + 
##     max.daily.air.gt26 + daily.ph.range + lag.no.large.fuc.q:max.daily.air.gt26, 
##     data = two.month)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.507 -16.291  -8.085   6.955  82.769 
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                           43.36859   17.13167   2.531   0.0168 *
## lag.no.large.fuc.q                    -0.56716    1.08098  -0.525   0.6037  
## daily.sal.range                       -1.98120    1.40581  -1.409   0.1690  
## max.daily.air.gt26                    -0.39524    0.78034  -0.506   0.6162  
## daily.ph.range                        -0.97481    2.58670  -0.377   0.7089  
## lag.no.large.fuc.q:max.daily.air.gt26  0.03345    0.08203   0.408   0.6863  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25.8 on 30 degrees of freedom
##   (25 observations deleted due to missingness)
## Multiple R-squared:  0.1027, Adjusted R-squared:  -0.04683 
## F-statistic: 0.6868 on 5 and 30 DF,  p-value: 0.6371
anova(m1)
## Analysis of Variance Table
## 
## Response: no.small.fuc.q
##                                       Df  Sum Sq Mean Sq F value  Pr(>F)  
## lag.no.large.fuc.q                     1    18.3   18.33  0.0275 0.86933  
## daily.sal.range                        1  2031.7 2031.73  3.0516 0.09089 .
## max.daily.air.gt26                     1    46.9   46.91  0.0705 0.79249  
## daily.ph.range                         1    78.8   78.78  0.1183 0.73326  
## lag.no.large.fuc.q:max.daily.air.gt26  1   110.7  110.70  0.1663 0.68635  
## Residuals                             30 19973.5  665.78                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(m1)
##                                            2.5 %     97.5 %
## (Intercept)                            8.3810510 78.3561243
## lag.no.large.fuc.q                    -2.7748120  1.6404978
## daily.sal.range                       -4.8522385  0.8898369
## max.daily.air.gt26                    -1.9889056  1.1984253
## daily.ph.range                        -6.2575546  4.3079373
## lag.no.large.fuc.q:max.daily.air.gt26 -0.1340773  0.2009735

Table

tab_model(m1)
  no.small.fuc.q
Predictors Estimates CI p
(Intercept) 43.37 8.38 – 78.36 0.017
lag.no.large.fuc.q -0.57 -2.77 – 1.64 0.604
daily.sal.range -1.98 -4.85 – 0.89 0.169
max.daily.air.gt26 -0.40 -1.99 – 1.20 0.616
daily.ph.range -0.97 -6.26 – 4.31 0.709
lag.no.large.fuc.q *
max.daily.air.gt26 		0.03 -0.13 – 0.20 0.686
Observations 36
R2 / R2 adjusted 0.103 / -0.047 
tab_model(m1,
  show.ci = 0.95,
  show.se = TRUE,
  show.stat = FALSE,
  emph.p = TRUE,
  df.method = "kr",
  p.style = "numeric_stars",
  title = "Linear model for small thalli density",
  linebreak = TRUE,
  dv.labels = c("small thalli density"
      ),
  pred.labels = c(
    "(Intercept)", 
    "large thalli density",
    "daily salinity range",
    "number of days with a maximum daily air temperature greater than 26C",
    "daily pH range",
    "large thalli density : number of days with a maximum daily air temperature greater than 26C"
      ),
  CSS = css_theme("cells")
  )
Linear model for small thalli density
  small thalli density
Predictors Estimates std. Error CI p
(Intercept) 43.37 * 17.13 8.38 – 78.36 0.017
large thalli density -0.57 1.08 -2.77 – 1.64 0.604
daily salinity range -1.98 1.41 -4.85 – 0.89 0.169
number of days with a maximum daily air temperature greater than 26C -0.40 0.78 -1.99 – 1.20 0.616
daily pH range -0.97 2.59 -6.26 – 4.31 0.709
large thalli density : number of days with a maximum daily air temperature greater than 26C 0.03 0.08 -0.13 – 0.20 0.686
Observations 36
R2 / R2 adjusted 0.103 / -0.047
  • p<0.05   ** p<0.01   *** p<0.001

Second ranked
m2 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.sal.lt10 + max.daily.air.gt26 + max.daily.ph.gt8 + lag.no.large.fuc.q:max.daily.air.gt26, data=one.month) Fitting and analyzing the top ranked model.

m2 <- lm(no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.sal.lt10 + max.daily.air.gt26 + max.daily.ph.gt8 + lag.no.large.fuc.q:max.daily.air.gt26, data=one.month)
plot(m2) 

m2
## 
## Call:
## lm(formula = no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.sal.lt10 + 
##     max.daily.air.gt26 + max.daily.ph.gt8 + lag.no.large.fuc.q:max.daily.air.gt26, 
##     data = one.month)
## 
## Coefficients:
##                           (Intercept)                     lag.no.large.fuc.q  
##                             15.677698                               1.537695  
##                    max.daily.sal.lt10                     max.daily.air.gt26  
##                             -0.384003                               0.051113  
##                      max.daily.ph.gt8  lag.no.large.fuc.q:max.daily.air.gt26  
##                              0.676607                               0.006729
summary(m2)
## 
## Call:
## lm(formula = no.small.fuc.q ~ lag.no.large.fuc.q + max.daily.sal.lt10 + 
##     max.daily.air.gt26 + max.daily.ph.gt8 + lag.no.large.fuc.q:max.daily.air.gt26, 
##     data = one.month)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -44.297 -18.246 -11.544   8.979 117.353 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)                           15.677698  13.880438   1.129    0.267
## lag.no.large.fuc.q                     1.537695   1.412855   1.088    0.284
## max.daily.sal.lt10                    -0.384003   2.730995  -0.141    0.889
## max.daily.air.gt26                     0.051113   0.890520   0.057    0.955
## max.daily.ph.gt8                       0.676607   0.899435   0.752    0.457
## lag.no.large.fuc.q:max.daily.air.gt26  0.006729   0.102540   0.066    0.948
## 
## Residual standard error: 36.23 on 33 degrees of freedom
##   (22 observations deleted due to missingness)
## Multiple R-squared:  0.06521,    Adjusted R-squared:  -0.07642 
## F-statistic: 0.4604 on 5 and 33 DF,  p-value: 0.8027
anova(m2)
## Analysis of Variance Table
## 
## Response: no.small.fuc.q
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## lag.no.large.fuc.q                     1   1580 1580.11  1.2039 0.2805
## max.daily.sal.lt10                     1     88   87.93  0.0670 0.7974
## max.daily.air.gt26                     1    611  610.56  0.4652 0.5000
## max.daily.ph.gt8                       1    737  737.30  0.5617 0.4589
## lag.no.large.fuc.q:max.daily.air.gt26  1      6    5.65  0.0043 0.9481
## Residuals                             33  43313 1312.53
confint(m2)
##                                             2.5 %     97.5 %
## (Intercept)                           -12.5622666 43.9176619
## lag.no.large.fuc.q                     -1.3367796  4.4121702
## max.daily.sal.lt10                     -5.9402538  5.1722480
## max.daily.air.gt26                     -1.7606630  1.8628891
## max.daily.ph.gt8                       -1.1533069  2.5065202
## lag.no.large.fuc.q:max.daily.air.gt26  -0.2018911  0.2153492

Table

tab_model(m2)
  no.small.fuc.q
Predictors Estimates CI p
(Intercept) 15.68 -12.56 – 43.92 0.267
lag.no.large.fuc.q 1.54 -1.34 – 4.41 0.284
max.daily.sal.lt10 -0.38 -5.94 – 5.17 0.889
max.daily.air.gt26 0.05 -1.76 – 1.86 0.955
max.daily.ph.gt8 0.68 -1.15 – 2.51 0.457
lag.no.large.fuc.q *
max.daily.air.gt26 		0.01 -0.20 – 0.22 0.948
Observations 39
R2 / R2 adjusted 0.065 / -0.076

not sure what a negative R^2 adjusted means

Trying similar approach thing with adults

I’m curious to see if a similar of models has an effect on adult/large thalli. I’m curious to see if the same or different conditions affect adult thalli and if it provides some insight about a bottle neck event that occurs for juveniles. Not sure if it makes sense to have the lag adult density so I’m going to remove it.

Envirnomental data for 30 days starting 60 days from survey date (excludes the 30 days leading up to the survey date). Adult density is the mean.

Read in data

#read in data
two.month<-read.csv(
      "https://raw.githubusercontent.com/Cmwegener/thesis/master/data/environment_field/envi.field.all.two.month.lag.csv",
    header = TRUE
  )
two.month$date<-as.Date(two.month$date, format = c("%Y-%m-%d"))

AIC

#models
m1 <- lm(no.large.fuc.q ~ daily.sal.range + max.daily.air.gt26 + daily.ph.range, data=two.month)

m2 <- lm(no.large.fuc.q ~ max.daily.sal.lt10 + max.daily.air.gt26 + max.daily.ph.gt8, data=two.month)

m3 <- lm(no.large.fuc.q ~ daily.sal.range + max.daily.air.gt26, data=two.month)

m4 <- lm(no.large.fuc.q ~ daily.sal.range + daily.ph.range, data=two.month)

m5 <- lm(no.large.fuc.q ~ max.daily.sal.lt10 + max.daily.ph.gt8, data=two.month)

m6 <- lm(no.large.fuc.q ~ max.daily.sal.lt10 + max.daily.air.gt26, data=two.month)

m7 <- lm(no.large.fuc.q ~ max.daily.air.gt26, data=two.month)


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

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

#calculate AIC of each model
aictab(cand.set = models, modnames = mod.names)
## 
## Model selection based on AICc:
## 
##    K   AICc Delta_AICc AICcWt Cum.Wt      LL
## m5 4 255.26       0.00   0.50   0.50 -123.09
## m1 5 257.21       1.95   0.19   0.68 -122.77
## m2 5 257.32       2.06   0.18   0.86 -122.83
## m4 4 257.82       2.57   0.14   1.00 -124.37
## m6 4 333.48      78.22   0.00   1.00 -162.36
## m7 3 336.47      81.21   0.00   1.00 -165.01
## m3 4 338.68      83.42   0.00   1.00 -164.96

Top ranked
m5 <- lm(no.large.fuc.q ~ max.daily.sal.lt10 + max.daily.ph.gt8, data=two.month)

Fitting and analyzing the top ranked model.

m5 <- lm(no.large.fuc.q ~ max.daily.sal.lt10 + max.daily.ph.gt8, data=two.month)

plot(m5) 

m5
## 
## Call:
## lm(formula = no.large.fuc.q ~ max.daily.sal.lt10 + max.daily.ph.gt8, 
##     data = two.month)
## 
## Coefficients:
##        (Intercept)  max.daily.sal.lt10    max.daily.ph.gt8  
##             7.2260              0.4215             -0.1091
summary(m5)
## 
## Call:
## lm(formula = no.large.fuc.q ~ max.daily.sal.lt10 + max.daily.ph.gt8, 
##     data = two.month)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0880 -2.9623 -0.5797  1.7961 15.6740 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          7.2260     0.9960   7.255  9.6e-09 ***
## max.daily.sal.lt10   0.4215     0.2052   2.054   0.0467 *  
## max.daily.ph.gt8    -0.1091     0.0934  -1.168   0.2497    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.706 on 39 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.1132, Adjusted R-squared:  0.06771 
## F-statistic: 2.489 on 2 and 39 DF,  p-value: 0.0961
anova(m5)
## Analysis of Variance Table
## 
## Response: no.large.fuc.q
##                    Df Sum Sq Mean Sq F value  Pr(>F)  
## max.daily.sal.lt10  1  80.00  79.997  3.6125 0.06476 .
## max.daily.ph.gt8    1  30.23  30.235  1.3653 0.24971  
## Residuals          39 863.63  22.144                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(m5)
##                          2.5 %     97.5 %
## (Intercept)         5.21142291 9.24066561
## max.daily.sal.lt10  0.00644685 0.83659773
## max.daily.ph.gt8   -0.29804480 0.07978105

Table

tab_model(m5)
  no.large.fuc.q
Predictors Estimates CI p
(Intercept) 7.23 5.21 – 9.24 <0.001
max.daily.sal.lt10 0.42 0.01 – 0.84 0.047
max.daily.ph.gt8 -0.11 -0.30 – 0.08 0.250
Observations 42
R2 / R2 adjusted 0.113 / 0.068

Trying this again with without the ph terms because there’s no compelling lit evidence to support it. More interested in salniity events and heat. I’d like to compare what effects adults vs juveniles so should I keep the pH terms? AIC

#models
m1 <- lm(no.large.fuc.q ~ daily.sal.range + max.daily.air.gt26, data=two.month)

m2 <- lm(no.large.fuc.q ~ max.daily.sal.lt10 + max.daily.air.gt26, data=two.month)

m3 <- lm(no.large.fuc.q ~ daily.sal.range, data=two.month)

m4 <- lm(no.large.fuc.q ~ max.daily.sal.lt10, data=two.month)

m5 <- lm(no.large.fuc.q ~ max.daily.air.gt26, data=two.month)



#define list of models
models <- list(m1,m2,m3,m4,m5)

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

#calculate AIC of each model
aictab(cand.set = models, modnames = mod.names)
## 
## Model selection based on AICc:
## 
##    K   AICc Delta_AICc AICcWt Cum.Wt      LL
## m4 3 332.63       0.00   0.51   0.51 -163.09
## m2 4 333.48       0.85   0.33   0.84 -162.36
## m5 3 336.47       3.84   0.07   0.92 -165.01
## m3 3 337.02       4.39   0.06   0.98 -165.28
## m1 4 338.68       6.05   0.02   1.00 -164.96
tab_model(m4)
  no.large.fuc.q
Predictors Estimates CI p
(Intercept) 6.58 5.36 – 7.80 <0.001
max.daily.sal.lt10 0.38 0.02 – 0.74 0.040
Observations 57
R2 / R2 adjusted 0.074 / 0.057

Positive effect of low salinity events for larger thalli density??? Not sure about my approach here. More of a curiousty quest but it’d be nice in the discussion to compare what events affects juveniles vs adults.