Goal of this markdown: Using a mixed model approach to analyze my field data. In this markdown I’m just focusing on the field data (ie. abundance, reproduction) and will address the environmental data (salinity, ect.) in a separate document.

What are the next steps?
- Should I add interactions terms?
- Confirm interpretation of results
- We talked about changing the family of distribution but do we need to do that if the transformations helped the fit?
- Should I start doing this with the environmental data?

What’s new:
- Percent cover as number instead of percentage (ie. 50% = 50) - Attempted to compare models to a null, not sure this is neccessary. Performs a Likelihood Ratio Test and output AIC and BIC values. (https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q1/015590.html)
- Data transformations for model 2.2 and 3.6. Having a hard time finding and understanding how to do data transformations within lme4 so I’m just going to do data transformations in the df for now. Achieves the same goal I was hoping for a different coding approach so I’ll look into it more - In the process of reading link from Karina (https://cran.r-project.org/web/packages/lme4/vignettes/lmer.pdf )
- Still need to look into some of the warnings and messages I’m getting but I want to finish reading the paper

Questions:
- Not sure how to interpret slope with the transformation terms. For example, model 2.2 I squared the week number so now I’m getting two slopes (for week number and week numer ^2). Neither of the slopes pass through zero so it’s significant but I’m not sure how to interpret it. - Do I need interaction terms? Not sure how to do those with mixed models so this is something I need to look into.
- For research question 3 I just looked at week number instead of season. Season was a character variable so I’m not sure if it’s as insightful as week number. Maybe I should look into assigning a range of week numbers to a particular season? So that it somehow remains a continuous variable? I’m not sure if I’ll be able to do this

Looking into data transformations with lme4. I’d like to be able to do these in the equation instead of transforming my df.I’m not seeing a clear answer so I may have to.

Mixed model is the best approach because my data is “clustered/grouped” by site but the sites are random representatives of Central SFE (grouped = random). Code for random term: (1 | term).

For “time” variable I will be using number of weeks since the first survey (with the first survey week being 0) which will allow me to convert it to a numneric variable and use quadratic terms as needed. This also allows the model to interpolate back to a defined “0”.

For now I’m not looking at interactions terms. The third link below has some insight on how to approach interaction terms in mixed model. I think this will be more important when I’m looking at the environmental data.

Helpful resources: https://m-clark.github.io/mixed-models-with-R/introduction.html https://ourcodingclub.github.io/tutorials/mixed-models/ http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#model-definition
https://cran.r-project.org/web/packages/lme4/vignettes/lmer.pdf https://ademos.people.uic.edu/Chapter18.html understanding the output: https://campus.datacamp.com/courses/hierarchical-and-mixed-effects-models-in-r/linear-mixed-effect-models?ex=7 https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q1/subject.html#12692

Notes:
- you can compare models by using anova(model1, model2) and see which model accounts for more of the variance. This can also be used to compare a null model

Questions for analysis:
ABUNDANCE (site only, no time) —-> Not doing these. Check with Karina but I think we’ve concluded that is this not insightful so I will be skipping this for now. I think this was more a question when we were trying to understand what the data looked like more. If I do answer this, I don’t think it’d be a mixed model. I could probably to a simple lm with cover ~ site, look at anova
1. Is there a difference in Fucus abundance between sites?
1.1 Is there a difference in Fucus percent cover amoung sites?
1.2 Is there a difference in Fucus total density amoung sites?

ABUNDANCE (site and time)
2. Does Fucus abundance over time? (site:time)
2.1 Does Fucus percent cover change over time?
2.2 Does Fucus total density change over time? 2.3 Does large thalli density change over time?
2.4 Does small thalli density change over time? notes: Find best fit (Gaussain, Poisson (density), normal). Seasonality time by quadrat. Percent cover as a proporation (i.e. 50% = 0.5)

REPRODUCTIVE PHENOLOGY (with seasonality) –> for now I’m doing week number
3. Does reproduction change over time?
3.1 Does oogonia/receptacle change over time?
3.2 Does oogonia/thalli change over time?
3.3 Does conceptical/thalli change over time?
3.4 Does percent reproductive tissue change over time? 3.5 Does percent reproductive apices change over time?○ 3.6 Does number of conceptacles per receptacles change over time?

####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)

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

Read in data

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

Convert date to number of weeks of surveys. Starting at week 0 being the first survey. This will make it a contiuous variable and then I can add a quadratic term to it. This is different that “week number of the year” so I have to use a customized approach. I was able to find a similar question online: https://stackoverflow.com/questions/54406803/assigning-custom-week-number-in-r/54407091, https://stackoverflow.com/questions/19780415/lubridate-week-to-find-consecutive-week-number-for-multi-year-periods

# subtracts minimun week number, +0 sets what week you want to start it on
all$week_nr <- (interval(min(all$date), all$date) %/% weeks(1)) + 0

2. Is there a difference in Fucus abundance over time?

2.1 Is there a difference in Fucus percent cover over time?
GOOD FIT (confirmed)- no transformations needed
Cover ~ site(random) + time

mlm2.1<-lmer(cover~ week_nr + (1|site), data=all)
mlm2.1
## Linear mixed model fit by REML ['lmerMod']
## Formula: cover ~ week_nr + (1 | site)
##    Data: all
## REML criterion at convergence: 472.2036
## Random effects:
##  Groups   Name        Std.Dev.
##  site     (Intercept)  4.609  
##  Residual             11.360  
## Number of obs: 61, groups:  site, 4
## Fixed Effects:
## (Intercept)      week_nr  
##      43.281       -0.331
summary(mlm2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: cover ~ week_nr + (1 | site)
##    Data: all
## 
## REML criterion at convergence: 472.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.95103 -0.60564  0.04063  0.82081  1.85494 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept)  21.24    4.609  
##  Residual             129.06   11.360  
## Number of obs: 61, groups:  site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 43.28091    3.58459  12.074
## week_nr     -0.33101    0.07136  -4.639
## 
## Correlation of Fixed Effects:
##         (Intr)
## week_nr -0.649
anova(mlm2.1)
## Analysis of Variance Table
##         npar Sum Sq Mean Sq F value
## week_nr    1 2776.8  2776.8  21.516
plot(mlm2.1)

confint(mlm2.1)
## Computing profile confidence intervals ...
##                  2.5 %    97.5 %
## .sig01       0.0000000 11.199785
## .sigma       9.4712927 13.689329
## (Intercept) 36.1784747 50.418449
## week_nr     -0.4716595 -0.189689
ranef(mlm2.1)$site
##        (Intercept)
## by.rb     3.766621
## hs.fp    -4.092564
## nd.eos   -2.533065
## pc.cc     2.859007
#comparing to a null? cover ~ site
null<-lmer(cover~ (1|site), data=all)
#performs a Likelihood Ratio Test (LRT) and prints values of AIC and BIC
anova(mlm2.1, null)
## refitting model(s) with ML (instead of REML)
## Data: all
## Models:
## null: cover ~ (1 | site)
## mlm2.1: cover ~ week_nr + (1 | site)
##        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## null      3 496.92 503.26 -245.46   490.92                         
## mlm2.1    4 480.45 488.89 -236.22   472.45 18.477  1   1.72e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#another way of confirming week number significantly affects percent cover?

Interpretation of output: There is a difference in Fucus percent cover over time. The intercept is significant with a value of 43% (95% confidence interval 36.18 to 50.42). Percent cover is slighting decreasing over time with a slope of -.33 meaning the decrease was about .33% per week.

2.2 Is there a difference in Fucus total density over time?
Good fit- with transformations density ~ site (random) + time

mlm2.2<-lmer(no.fuc.q ~ week_nr + (1|site), data=all)
plot(mlm2.2)

Transformation
sqrt(no.fuc.q ~ week_nr + week_nr^2 +((1|site)

#square root transform total density
all$sqrt.no.fuc.q<-sqrt(all$no.fuc.q)

#square week number
all$sqr.week_nr<-all$week_nr^2

#model
mlm2.2<-lmer(sqrt.no.fuc.q ~ week_nr + sqr.week_nr + (1|site), data=all)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
plot(mlm2.2) #looks much better!

mlm2.2
## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt.no.fuc.q ~ week_nr + sqr.week_nr + (1 | site)
##    Data: all
## REML criterion at convergence: 264.6063
## Random effects:
##  Groups   Name        Std.Dev.
##  site     (Intercept) 0.9207  
##  Residual             1.7505  
## Number of obs: 61, groups:  site, 4
## Fixed Effects:
## (Intercept)      week_nr  sqr.week_nr  
##    8.576820    -0.163591     0.001625  
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
summary(mlm2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt.no.fuc.q ~ week_nr + sqr.week_nr + (1 | site)
##    Data: all
## 
## REML criterion at convergence: 264.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9615 -0.6973 -0.1418  0.6245  2.8434 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.8477   0.9207  
##  Residual             3.0642   1.7505  
## Number of obs: 61, groups:  site, 4
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)  8.5768195  0.7265103  11.806
## week_nr     -0.1635911  0.0414342  -3.948
## sqr.week_nr  0.0016248  0.0006303   2.578
## 
## Correlation of Fixed Effects:
##             (Intr) wek_nr
## week_nr     -0.622       
## sqr.week_nr  0.509 -0.964
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
anova(mlm2.2)
## Analysis of Variance Table
##             npar Sum Sq Mean Sq F value
## week_nr        1 93.103  93.103 30.3842
## sqr.week_nr    1 20.360  20.360  6.6446
confint(mlm2.2)
## Computing profile confidence intervals ...
##                     2.5 %       97.5 %
## .sig01       0.2387882508  2.125503456
## .sigma       1.4463085387  2.090183145
## (Intercept)  7.1614415982  9.993792120
## week_nr     -0.2448236486 -0.082579515
## sqr.week_nr  0.0003930787  0.002861275
ranef(mlm2.2)$site
##        (Intercept)
## by.rb    0.7133954
## hs.fp   -1.1931716
## nd.eos   0.1816081
## pc.cc    0.2981681

Interpretation of output: Total Fucus density does change over time. Intercept 8.58 is within the 95% confidence interval (7.16 to 9.99) and the slope does not pass through zero. Need help interpretting with the transformation (losing 0.16 fucus a week?).

2.3 Does adult density change over time?
Good fit
no.large.fuc.q ~ time + site(random)

mlm2.3<-lmer(no.large.fuc.q~ week_nr + (1|site), data=all)
## boundary (singular) fit: see ?isSingular
plot(mlm2.3)

mlm2.3
## Linear mixed model fit by REML ['lmerMod']
## Formula: no.large.fuc.q ~ week_nr + (1 | site)
##    Data: all
## REML criterion at convergence: 322.7336
## Random effects:
##  Groups   Name        Std.Dev.
##  site     (Intercept) 0.000   
##  Residual             4.006   
## Number of obs: 57, groups:  site, 4
## Fixed Effects:
## (Intercept)      week_nr  
##     10.6130      -0.1026  
## optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
summary(mlm2.3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: no.large.fuc.q ~ week_nr + (1 | site)
##    Data: all
## 
## REML criterion at convergence: 322.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7479 -0.7366 -0.0400  0.5788  3.2462 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept)  0.00    0.000   
##  Residual             16.05    4.006   
## Number of obs: 57, groups:  site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 10.61298    1.10413   9.612
## week_nr     -0.10259    0.02776  -3.695
## 
## Correlation of Fixed Effects:
##         (Intr)
## week_nr -0.877
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
anova(mlm2.3)
## Analysis of Variance Table
##         npar Sum Sq Mean Sq F value
## week_nr    1 219.15  219.15  13.654
confint(mlm2.3)
## Computing profile confidence intervals ...
##                  2.5 %      97.5 %
## .sig01       0.0000000  1.71489927
## .sigma       3.3101735  4.78512466
## (Intercept)  8.4510800 12.77462914
## week_nr     -0.1569484 -0.04822257
ranef(mlm2.3)$site
##        (Intercept)
## by.rb            0
## hs.fp            0
## nd.eos           0
## pc.cc            0

Interpretation of output The intercept 10.61 is between the 95% confidence interval (8.45 to 12.77) and the slope -0.10 does not pass through zero, meaning large Fucus density slightly decreases over time. A loss of .10 large thalli means 1 is lost ~10 weeks, a very slow rate. This supports the idea that large fucus thalli provide a farily constant and stable habitat for invertebrates (vs small thalli that decrease). This also supports that large thalli more significantly contribute to percent cover while small thalli contribute to total density.

2.4 Is juvenile density affected by the season?
Good fit- with transformation
no.small.fuc.q ~ time + site(random)

mlm2.4<-lmer(no.small.fuc.q~ week_nr + (1|site), data=all)
plot(mlm2.4)

Doesn’t look like a great fit. Pretty much looks identical to total density which isn’t surprising since small thalli largely influcence the total density. Going to try the same transformation that worked for total density

Transformation
sqrt(no.small.fuc.q ~ week_nr + week_nr^2 +((1|site)

#square root transform small thalli density
all$sqrt.small.fuc.q<-sqrt(all$no.small.fuc.q)

#square week number
all$sqr.week_nr<-all$week_nr^2

#model
mlm2.4<-lmer(sqrt.small.fuc.q ~ week_nr + sqr.week_nr + (1|site), data=all)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
plot(mlm2.4) #looks much better!

mlm2.4
## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt.small.fuc.q ~ week_nr + sqr.week_nr + (1 | site)
##    Data: all
## REML criterion at convergence: 239.0961
## Random effects:
##  Groups   Name        Std.Dev.
##  site     (Intercept) 1.106   
##  Residual             1.590   
## Number of obs: 57, groups:  site, 4
## Fixed Effects:
## (Intercept)      week_nr  sqr.week_nr  
##   10.247571    -0.325272     0.003775  
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
summary(mlm2.4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt.small.fuc.q ~ week_nr + sqr.week_nr + (1 | site)
##    Data: all
## 
## REML criterion at convergence: 239.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.32751 -0.72528 -0.06033  0.63549  2.35548 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 1.224    1.106   
##  Residual             2.528    1.590   
## Number of obs: 57, groups:  site, 4
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) 10.2475713  0.8657281  11.837
## week_nr     -0.3252723  0.0453725  -7.169
## sqr.week_nr  0.0037751  0.0006538   5.774
## 
## Correlation of Fixed Effects:
##             (Intr) wek_nr
## week_nr     -0.669       
## sqr.week_nr  0.578 -0.970
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
anova(mlm2.4)
## Analysis of Variance Table
##             npar  Sum Sq Mean Sq F value
## week_nr        1 105.141 105.141  41.593
## sqr.week_nr    1  84.285  84.285  33.342
confint(mlm2.4)
## Computing profile confidence intervals ...
##                    2.5 %      97.5 %
## .sig01       0.411659277  2.47614374
## .sigma       1.303977627  1.91067909
## (Intercept)  8.569383190 11.92719657
## week_nr     -0.414146987 -0.23649205
## sqr.week_nr  0.002496403  0.00505627
ranef(mlm2.4)$site
##        (Intercept)
## by.rb    0.9490921
## hs.fp   -1.4722460
## nd.eos   0.2542086
## pc.cc    0.2689453

Interpretation of output: Small thalli density does change over time. Intercept 10.25 is between 95% confidence interval (8.57 to 11.93) and the slope does not go through zero. Small thalli do not contribute to a constant habitat for invertebrates even though they contribute largely to total density.

3. Is reproduction affectd by time?

3.1 Is oogonia/conceptacle affected by the season?
GOOD FIT (confirmed)- no transformations needed
avg.oog ~ time + site(random)

mlm3.1<-lmer(oog.per.con~ week_nr + (1|site), data=all)
mlm3.1
## Linear mixed model fit by REML ['lmerMod']
## Formula: oog.per.con ~ week_nr + (1 | site)
##    Data: all
## REML criterion at convergence: 473.46
## Random effects:
##  Groups   Name        Std.Dev.
##  site     (Intercept)  6.526  
##  Residual             15.227  
## Number of obs: 57, groups:  site, 4
## Fixed Effects:
## (Intercept)      week_nr  
##     30.1167      -0.1847
summary(mlm3.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: oog.per.con ~ week_nr + (1 | site)
##    Data: all
## 
## REML criterion at convergence: 473.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3848 -0.8534 -0.1623  0.5926  2.2798 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept)  42.59    6.526  
##  Residual             231.85   15.227  
## Number of obs: 57, groups:  site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  30.1167     4.9963   6.028
## week_nr      -0.1847     0.1054  -1.751
## 
## Correlation of Fixed Effects:
##         (Intr)
## week_nr -0.640
anova(mlm3.1)
## Analysis of Variance Table
##         npar Sum Sq Mean Sq F value
## week_nr    1 711.05  711.05  3.0668
plot(mlm3.1)

confint(mlm3.1)
## Computing profile confidence intervals ...
##                 2.5 %      97.5 %
## .sig01       0.000000 15.78872110
## .sigma      12.610099 18.47927351
## (Intercept) 20.125452 40.02032213
## week_nr     -0.393969  0.02299488
ranef(mlm3.1)$site
##        (Intercept)
## by.rb    -3.625301
## hs.fp     7.187345
## nd.eos    1.475826
## pc.cc    -5.037871

Interpretation of output: The number of oogonia per conceptacle is not affected by time/ does not significantlly change over time. Although the intercept value of 30.12 is within the 95% confidence interval (20.13 to 40.02) the slope crosses zero making the slope effectively horizantal meaning it is not changing over time. This is different from what is found in the literature where populations have clear reproductive peaks and seasonality. This suggests that Fucus distichus populations in Centeral SFE are dribble spawners and reproduce throughtout the year.

3.2 Is oogonia/thalli affected by the season?
Good fit- look into pulling outlier and look into reason for split in data

oog.thalli ~ time + site(random)

mlm3.2<-lmer(oog.thalli~ week_nr + (1|site), data=all)
plot(mlm3.2)

mlm3.2
## Linear mixed model fit by REML ['lmerMod']
## Formula: oog.thalli ~ week_nr + (1 | site)
##    Data: all
## REML criterion at convergence: 1570.186
## Random effects:
##  Groups   Name        Std.Dev.
##  site     (Intercept) 142500  
##  Residual             325272  
## Number of obs: 57, groups:  site, 4
## Fixed Effects:
## (Intercept)      week_nr  
##      253623        -1181
summary(mlm3.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: oog.thalli ~ week_nr + (1 | site)
##    Data: all
## 
## REML criterion at convergence: 1570.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1956 -0.4123 -0.2753  0.0754  5.1008 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  site     (Intercept) 2.031e+10 142500  
##  Residual             1.058e+11 325272  
## Number of obs: 57, groups:  site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   253623     107750   2.354
## week_nr        -1180       2252  -0.524
## 
## Correlation of Fixed Effects:
##         (Intr)
## week_nr -0.634
anova(mlm3.2)
## Analysis of Variance Table
##         npar    Sum Sq   Mean Sq F value
## week_nr    1 2.906e+10 2.906e+10  0.2747
confint(mlm3.2)
## Computing profile confidence intervals ...
##                  2.5 %     97.5 %
## .sig01           0.000 342616.962
## .sigma      269368.769 394669.288
## (Intercept)  37897.736 467591.042
## week_nr      -5649.412   3256.497
ranef(mlm3.2)$site
##        (Intercept)
## by.rb    -74813.99
## hs.fp    178270.77
## nd.eos   -21403.14
## pc.cc    -82053.65

Interpretation of output: The intercept 253623 is outside the 95% confidence interval (37897 to 394669) and the slope -1181 is effectively zero (-5649 to 3256), meaning there is no significant change in the number of oogonia per thalli over time. This again supports the idea that Fucus in SFE are dribble spawners without clear reproductive seasonality or peaks.

Why is the data split? I’m going to subset the data by site. Not sure if this is the best way to do this.

#subset by site
HS<-subset(all, site == "hs.fp")
BY<-subset(all, site == "by.rb")
ND<-subset(all, site == "nd.eos")
PC<-subset(all, site == "pc.cc")

#can't use lmer since there's no random factor. Looking at residual vs fitted plot
hs<-lm(oog.thalli~ week_nr, data=HS)
plot(hs) #300,000 - 600,000

by<-lm(oog.thalli~ week_nr, data=BY)
plot(by) #90,000 - 150,000

nd<-lm(oog.thalli~ week_nr, data=ND)
plot(nd) #50,000 - 350,000

pc<-lm(oog.thalli~ week_nr, data=PC)
plot(pc) # 0 - 200,000

Horseshoe Bay contributing to most of the points on the right of the graph. Brickyard Park and Paradise Cay contributing more to the left side. North Dock spanning both sides? I’m really not sure if this is the best way to look at this

3.3 Is conceptical/thalli affected by the season?
Good fit- look into removing outlier, look into error message

con.thalli ~ time + site(random)

mlm3.3<-lmer(con.thalli~ week_nr + (1|site), data=all)
plot(mlm3.3)

mlm3.3
## Linear mixed model fit by REML ['lmerMod']
## Formula: con.thalli ~ week_nr + (1 | site)
##    Data: all
## REML criterion at convergence: 1209.187
## Random effects:
##  Groups   Name        Std.Dev.
##  site     (Intercept)  4556   
##  Residual             12295   
## Number of obs: 57, groups:  site, 4
## Fixed Effects:
## (Intercept)      week_nr  
##    11494.49        52.72
summary(mlm3.3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: con.thalli ~ week_nr + (1 | site)
##    Data: all
## 
## REML criterion at convergence: 1209.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4127 -0.5385 -0.2877  0.3413  4.3648 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  site     (Intercept)  20754818  4556   
##  Residual             151173666 12295   
## Number of obs: 57, groups:  site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 11494.49    3810.67   3.016
## week_nr        52.72      85.14   0.619
## 
## Correlation of Fixed Effects:
##         (Intr)
## week_nr -0.678
anova(mlm3.3)
## Analysis of Variance Table
##         npar   Sum Sq  Mean Sq F value
## week_nr    1 57966650 57966650  0.3834
confint(mlm3.3) 
## Computing profile confidence intervals ...
##                  2.5 %     97.5 %
## .sig01          0.0000 11436.0272
## .sigma      10182.3117 14920.2082
## (Intercept)  3937.4548 18965.3398
## week_nr      -116.5206   220.1806
ranef(mlm3.3)$site
##        (Intercept)
## by.rb    -935.4477
## hs.fp    5390.3473
## nd.eos  -1419.5503
## pc.cc   -3035.3493

Interpretation of output: Intercept with 95% ci 11494.49 (95%ci 3937.45 to 18965.34), slope passes zero so conceptacle per thalli does not signficantly change over time –> dribble spawner.

3.4 Is percent reproductive tissue affected by the season?
Good fit- look into split in data

perc.rdw ~ time + site(random)

mlm3.4<-lmer(perc.rdw~ week_nr + (1|site), data=all)
plot(mlm3.4)

mlm3.4
## Linear mixed model fit by REML ['lmerMod']
## Formula: perc.rdw ~ week_nr + (1 | site)
##    Data: all
## REML criterion at convergence: 451.767
## Random effects:
##  Groups   Name        Std.Dev.
##  site     (Intercept)  8.243  
##  Residual             12.258  
## Number of obs: 57, groups:  site, 4
## Fixed Effects:
## (Intercept)      week_nr  
##    18.83654      0.02305
summary(mlm3.4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: perc.rdw ~ week_nr + (1 | site)
##    Data: all
## 
## REML criterion at convergence: 451.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9107 -0.6444 -0.1481  0.5326  2.6233 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept)  67.94    8.243  
##  Residual             150.26   12.258  
## Number of obs: 57, groups:  site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 18.83654    5.12528   3.675
## week_nr      0.02305    0.08489   0.271
## 
## Correlation of Fixed Effects:
##         (Intr)
## week_nr -0.503
anova(mlm3.4)
## Analysis of Variance Table
##         npar Sum Sq Mean Sq F value
## week_nr    1 11.075  11.075  0.0737
confint(mlm3.4)
## Computing profile confidence intervals ...
##                  2.5 %     97.5 %
## .sig01       2.8769867 18.5365830
## .sigma      10.1516123 14.8778297
## (Intercept)  8.2502663 29.4114459
## week_nr     -0.1448988  0.1907738
ranef(mlm3.4)$site
##        (Intercept)
## by.rb    -7.941243
## hs.fp     4.187305
## nd.eos    8.536717
## pc.cc    -4.782779

Interpretation of output: Although the intercept value of 18.84 is within the 95% confidence interval (8.25 to 29.41) the slope passes through zero (-0.14 to 0.19) meaning percent reproductive tissue does not change over time. This again supports the idea that SFE Fucus populations are dribble spawners. This contrasts from what is published in the literature with other populations having distinct reproductive seasons or peaks.

3.5 Does percent reproductive apices change over time? Good fit- no transformations needed perc.ra ~ time + site(random)

mlm3.5<-lmer(perc.ra~ week_nr + (1|site), data=all)
plot(mlm3.5) #good fit

mlm3.5
## Linear mixed model fit by REML ['lmerMod']
## Formula: perc.ra ~ week_nr + (1 | site)
##    Data: all
## REML criterion at convergence: 478.3566
## Random effects:
##  Groups   Name        Std.Dev.
##  site     (Intercept) 10.38   
##  Residual             15.62   
## Number of obs: 57, groups:  site, 4
## Fixed Effects:
## (Intercept)      week_nr  
##    26.34445      0.05378
summary(mlm3.5)
## Linear mixed model fit by REML ['lmerMod']
## Formula: perc.ra ~ week_nr + (1 | site)
##    Data: all
## 
## REML criterion at convergence: 478.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4864 -0.7186 -0.1143  0.4345  2.3079 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 107.8    10.38   
##  Residual             243.9    15.62   
## Number of obs: 57, groups:  site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 26.34445    6.48264   4.064
## week_nr      0.05378    0.10817   0.497
## 
## Correlation of Fixed Effects:
##         (Intr)
## week_nr -0.506
anova(mlm3.5)
## Analysis of Variance Table
##         npar Sum Sq Mean Sq F value
## week_nr    1 60.295  60.295  0.2472
confint(mlm3.5)
## Computing profile confidence intervals ...
##                  2.5 %     97.5 %
## .sig01       3.5905486 23.3765899
## .sigma      12.9343695 18.9556625
## (Intercept) 12.9665303 39.7066407
## week_nr     -0.1602124  0.2674688
ranef(mlm3.5)$site
##        (Intercept)
## by.rb   -10.318643
## hs.fp     6.099902
## nd.eos   10.059337
## pc.cc    -5.840597

Interpretation of output: Although the intercept value of 26.34 is within the 95% confidence interval (12.97 to 39.71) the slope passes through zero (-0.16 to 0.26) meaning percent reproductive apices does not change over time. This again supports the idea that SFE Fucus populations are dribble spawners. ng distinct reproductive seasons or peaks.

3.6 Does number of conceptacles per receptacles change over time? Good fit- no transformations needed no.concept.recp ~ time + site(random)

mlm3.6<-lmer(no.concept.recp~ week_nr + (1|site), data=all)
plot(mlm3.6) #good fit, split by four sites?

mlm3.6
## Linear mixed model fit by REML ['lmerMod']
## Formula: no.concept.recp ~ week_nr + (1 | site)
##    Data: all
## REML criterion at convergence: 679.7148
## Random effects:
##  Groups   Name        Std.Dev.
##  site     (Intercept) 62.00   
##  Residual             97.62   
## Number of obs: 57, groups:  site, 4
## Fixed Effects:
## (Intercept)      week_nr  
##   292.29384      0.03721
summary(mlm3.6)
## Linear mixed model fit by REML ['lmerMod']
## Formula: no.concept.recp ~ week_nr + (1 | site)
##    Data: all
## 
## REML criterion at convergence: 679.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.09555 -0.68771  0.05725  0.68473  2.73033 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 3844     62.00   
##  Residual             9530     97.62   
## Number of obs: 57, groups:  site, 4
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 292.29384   39.36878   7.425
## week_nr       0.03721    0.67609   0.055
## 
## Correlation of Fixed Effects:
##         (Intr)
## week_nr -0.521
anova(mlm3.6)
## Analysis of Variance Table
##         npar Sum Sq Mean Sq F value
## week_nr    1 28.861  28.861   0.003
confint(mlm3.6)
## Computing profile confidence intervals ...
##                  2.5 %     97.5 %
## .sig01       20.396311 140.291854
## .sigma       80.846027 118.476905
## (Intercept) 211.115166 372.900895
## week_nr      -1.303964   1.369501
ranef(mlm3.6)$site
##        (Intercept)
## by.rb   -25.306821
## hs.fp    82.270654
## nd.eos   -8.635105
## pc.cc   -48.328727

Interpretation of output: Although the intercept value of 292.29 is within the 95% confidence interval (211.12 to 372.90) the slope passes through zero (-1.30 to 1.37) meaning the number of conceptacles per receptacle does not change over time. This again supports the idea that SFE Fucus populations are dribble spawners. ng distinct reproductive seasons or peaks.