Field mlm visulatizations

Goal: To make visualizations for my mixed linear models for my field data. Mixed models here: https://cmwegener.github.io/thesis/mixed_models_field.html. I numbered/labeled the models the same in both markdowns to make it easier to read the two (hopefully).

Components:
- First half of markdown- Graph with predicted values for each week. X-axis is the week number, y-axis is the predicted value (ex. density) and each point with standard error shading. In the background will the the raw data points. I also have simplified equations calculed from model outputs on each graph.
- Second half- Table of results with fixed and random effects. Look in the literature for examples. Look into tidy table tools and the tables I’ve made in R previously.

What’s new:
- added percent reproductive apices and number of conceptacles per receptacle
- added equations taken from model output

Questions:
- I’m still unclear about our transformed models. Having some trouble with models 2.2 (total density) and model 2.4 (small thalli density) which are the transformed models. The top graph is supposed to graph the overall predicted values of the model but it’s linear while the graph by site is curved. Can we just show the graph of the model per site for these? The raw data points I have behind the lines for these are the raw untransformed data points, should I be graphing the transformed data points instead? Also I simplified the equation to remove the squareroot on the y value but the intercept now isn’t matching what’s graphed. Overall think I need help with these two.
- There is no difference between sites (confirmed by table output) for model 2.3 (large thalli density) so the graph of predicted values persite looks like ther’s one line. I had to jitter it a considerable amount to be able to see the four lines. Can we just show the overall predicted values model and state that there’s no difference between sites?

• Table outputting p-values, why? https://cran.r-project.org/web/packages/sjPlot/vignettes/tab_mixed.html
• Only output (small) difference between model and table is the confidence intervals, why?
  https://stats.stackexchange.com/questions/146988/getting-degrees-of-freedom-from-lmer
  https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html
  http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#why-doesnt-lme4-display-denominator-degrees-of-freedomp-values-what-other-options-do-i-have
  
• I found a way to specify the method to calculate df but I’m not sure the best approach. From ?tab_model df.method and p.val let’s you choose the computing method:
  “Method for computing degrees of freedom for p-values, standard errors and confidence intervals (CI). Only applies to mixed models. Use df.method =”wald" for a faster, but less precise computation. df.method = “kenward” (or df.method = “kr”) uses Kenward-Roger approximation for the degrees of freedom. df.method = “satterthwaite” uses Satterthwaite’s approximation and “ml1” uses a “m-l-1” heuristic see degrees_of_freedom for details). Use show.df = TRUE to show the approximated degrees of freedom for each coefficient."

Useful links:
https://lmudge13.github.io/sample_code/mixed_effects.html
https://stackoverflow.com/questions/9447329/how-to-plot-the-results-of-a-mixed-model
https://terpconnect.umd.edu/~egurarie/teaching/Biol709/Topic3/Lecture16_MixedEffectsModels.html#mixed_effects_models
https://stackoverflow.com/questions/31075407/plot-mixed-effects-model-in-ggplot https://cran.r-project.org/web/packages/ggeffects/vignettes/ggeffects.html
https://aosmith.rbind.io/2018/11/16/plot-fitted-lines/
https://strengejacke.github.io/sjPlot/articles/tab_model_estimates.html
https://cran.r-project.org/web/packages/sjPlot/vignettes/tab_mixed.html
Customizing ledgend http://www.cookbook-r.com/Graphs/Legends_(ggplot2)/
https://www.zoology.ubc.ca/~schluter/R/Model.html
Printing table in different formats: https://github.com/strengejacke/sjPlot/issues/712

Setting up the same way I did in my mixed model markdown (https://cmwegener.github.io/thesis/mixed_models_field.html)

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)
library(AICcmodavg)
library(plm) #for random factors
library(lme4)
library(car)
library(xts)
library(ggeffects)

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

####Graph of predicted values over time plotted over raw data points#### Figuring out the code

Starting with the simplest model (Cover ~ site(random) + time). Plotting the model’s predicted value with a grey band that is +/- the standard error from the predicted value with the raw data behind it. https://ourcodingclub.github.io/tutorials/mixed-models/#presenting

#model
mlm2.1<-lmer(cover~ week_nr + (1|site), data=all)


# Extract the prediction data frame
pred.mm <- ggpredict(mlm2.1, terms = c("week_nr"))  # this gives overall predictions for the model

# Plot the predictions 

(
  ggplot(pred.mm) +
    geom_line(aes(x = x, y = predicted)) +  # slope
    geom_ribbon(
      aes(                                # grey band is +/- the standard error from predicted value
        x = x,
        ymin = predicted - std.error,
        ymax = predicted + std.error
      ),
      fill = "lightgrey",
      alpha = 0.5
    ) +  # error band
    geom_point(data = all,                      # adding the raw data (scaled values)
               aes(
                 x = week_nr, y = cover, colour = site
               )) +
    labs(x = "weeks since first survey", y = "Percent Cover",
         title = "Percent cover declined ~0.33% per week",
         subtitle = "percent cover ~ week since first survey + site(random)") +
    theme_minimal()
)

Linear regression line per site. This is just a linear model per site, not actually what our model predicts –> Not what we want

ggplot(data = all,
       aes(x = week_nr, y = cover, colour = site)) + geom_point() +
    geom_smooth(se=FALSE, method=lm)

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

The lines here are fit based on our model. Using: https://stats.stackexchange.com/questions/98958/plots-to-illustrate-results-of-linear-mixed-effect-model/99049

#not what we want
ggplot(all, aes(x=week_nr, y=cover, colour=site)) +
  geom_point() +
  geom_line(aes(y=predict(mlm2.1), group=site)) +
  geom_line(data=all, aes(y=predict(mlm2.1, level=0, all=all), size="SFE")) +
  scale_size_manual(name="Predictions", values=c("Site"=0.5, "SFE"=1)) +
  theme_bw(base_size=22) 

## Warning in predict.merMod(mlm2.1, level = 0, all = all): unused arguments
## ignored

Interesting that Brickyard Park and Paradise Cay are closer together and then Point Chauncy and Horseshoe are closer. Visually that’s how Byron and I thought they compared/more similar to

Lines based off of the predicted values from our model. Using: https://aosmith.rbind.io/2018/11/16/plot-fitted-lines/

#add predicted values from the model to the dataset
all$predlm = predict(mlm2.1)

#then use geom_line to add fitted lines based on this new column
ggplot(all, aes(x = week_nr, y = cover, color = site)) +
  geom_point() +
  geom_line(aes(y = predlm), size = 1) 

Stacking and working on aesthetics of plots

#model
mlm2.1<-lmer(cover~ week_nr + (1|site), data=all)

# Extract the prediction data frame
pred.mm <- ggpredict(mlm2.1, terms = c("week_nr"))  # this gives overall predictions for the model


# Plot the predictions of model 
model.plot <-
  ggplot(pred.mm) +
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(
    aes(
      x = x,
      ymin = predicted - std.error,
      ymax = predicted + std.error
    ),
    fill = "lightgrey",
    alpha = 0.5
  ) +
  labs(x = "",
       y = "",
       title = "Predicted values of model") +
  geom_point(data = all,
             aes(x = week_nr, y = cover, colour = site)) +
  scale_colour_discrete(                                    #key modifiers
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  theme_minimal()



#add predicted values from the model to the dataset
all$predlm = predict(mlm2.1)

#then use geom_line to add fitted lines based on this new column
site <-
  ggplot(all, aes(x = week_nr, y = cover, colour = site)) +
  geom_point() +
  scale_colour_discrete(
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  geom_line(aes(y = predlm), size = 1) +
  labs(x = "",
       y = "",
       title = "Predicted values of model per site") +
  theme_minimal()


#arrange plots
plot<-ggarrange(model.plot, site, ncol=1, nrow=2, common.legend = TRUE, legend="right")

#making subtitle for top
title <-
  expression(atop(
    bold("percent cover ~ week since first survey + site(random)"),
    scriptstyle("Percent cover declined ~0.33% per week")
  ))

#annotating figure
annotate_figure(
  plot,
  top = text_grob(title),
  left = text_grob(
    "Percent Cover",
    color = "black",
    rot = 90
  ),
  bottom = text_grob(
    "Weeks since first survey",
    color = "black",
    rot = 0
      )
)

Apply to all models!

2. Does Fucus abundance change over time?

2.1 Does Fucus percent cover change over time? Cover ~ site(random) + time

#model
mlm2.1<-lmer(cover~ week_nr + (1|site), data=all)

# Extract the prediction data frame
pred.mm <- ggpredict(mlm2.1, terms = c("week_nr"))  # this gives overall predictions for the model

# Plot the predictions of model 
model.plot <-
  ggplot(pred.mm) +
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(
    aes(
      x = x,
      ymin = predicted - std.error,
      ymax = predicted + std.error
    ),
    fill = "lightgrey",
    alpha = 0.5
  ) +
  labs(x = "",
       y = "",
       title = "Predicted values of model") +
  geom_point(data = all,
             aes(x = week_nr, y = cover, colour = site)) +
  scale_colour_discrete(                                    #key modifiers
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  theme_minimal()



#add predicted values from the model to the dataset
all$predlm = predict(mlm2.1)

#then use geom_line to add fitted lines based on this new column
site <-
  ggplot(all, aes(x = week_nr, y = cover, colour = site)) +
  geom_point() +
  scale_colour_discrete(
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  geom_line(aes(y = predlm), size = 1) +
  labs(x = "",
       y = "",
       title = "Predicted values of model per site") +
  theme_minimal()


#arrange plots
plot<-ggarrange(model.plot, site, ncol=1, nrow=2, common.legend = TRUE, legend="right")

#making subtitle for top
title <-
  expression(atop(
    bold("percent cover ~ week number + site(random)"),
    "percent cover = -0.33(week) + 43.28",
  ))

#annotating figure
annotate_figure(
  plot,
  top = text_grob(title),
  left = text_grob(
    "Percent Cover",
    color = "black",
    rot = 90
  ),
  bottom = text_grob(
    "Weeks since first survey",
    color = "black",
    rot = 0
      )
)

2.2 Does Fucus total density change over time?
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

# Extract the prediction data frame
pred.mm <- ggpredict(mlm2.2, terms = c("week_nr"))  # this gives overall predictions for the model, adding sqr.week_nr makes a zig zap shape to the top graph

# Plot the predictions of model 
model.plot <-
  ggplot(pred.mm) +
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(
    aes(
      x = x,
      ymin = predicted - std.error,
      ymax = predicted + std.error
    ),
    fill = "lightgrey",
    alpha = 0.5
  ) +
  labs(x = "",
       y = "",
       title = "Predicted values of model") +
  geom_point(data = all,
             aes(x = week_nr, y = sqrt.no.fuc.q, colour = site)) +
  scale_colour_discrete(                                    #key modifiers
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  theme_minimal()



#add predicted values from the model to the dataset
all$predlm = predict(mlm2.2)

#then use geom_line to add fitted lines based on this new column
site <-
  ggplot(all, aes(x = week_nr, y = sqrt.no.fuc.q, colour = site)) +
  geom_point() +
  scale_colour_discrete(
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  geom_line(aes(y = predlm), size = 1) +
  labs(x = "",
       y = "",
       title = "Predicted values of model per site") +
  theme_minimal()


#arrange plots
plot<-ggarrange(model.plot, site, ncol=1, nrow=2, common.legend = TRUE, legend="right")

#making subtitle for top
title <-
  expression(atop(
    bold("sqrt(total density) ~ week + week^2 + site(random)"),
    "total density = 0.0000027(week)^4 - 0.00051(week)^3 + .053(week)^2 -2.75(week) + 73.62",
  ))

#annotating figure
annotate_figure(
  plot,
  top = text_grob(title),
  left = text_grob(
    "Total Density (number of Fucus per 0.25^2)",
    color = "black",
    rot = 90
  ),
  bottom = text_grob(
    "Weeks since first survey",
    color = "black",
    rot = 0
      )
)

I can’t figure out how to add the curve to the top graph/overall predicted values from the model. Also the intercept that I calcualted doesn’t match what’s graphed so I need to look into this more.

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

Added position dodge to the site trend lines because they were stacked (https://stackoverflow.com/questions/10866047/jitter-geom-line)
There’s no difference between sites in this model so maybe we just show the top graph

#remove rows with NA 
no.na<-all[!is.na(all$no.large.fuc.q), ]

#model
mlm2.3<-lmer(no.large.fuc.q~ week_nr + (1|site), data=no.na)

## boundary (singular) fit: see ?isSingular

# Extract the prediction data frame
pred.mm <- ggpredict(mlm2.3, terms = c("week_nr"))  # this gives overall predictions for the model

# Plot the predictions of model 
model.plot <-
  ggplot(pred.mm) +
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(
    aes(
      x = x,
      ymin = predicted - std.error,
      ymax = predicted + std.error
    ),
    fill = "lightgrey",
    alpha = 0.5
  ) +
  labs(x = "",
       y = "",
       title = "Predicted values of model") +
  geom_point(data = no.na,
             aes(x = week_nr, y = no.large.fuc.q, colour = site)) +
  scale_colour_discrete(                                    #key modifiers
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  theme_minimal()



#add predicted values from the model to the dataset
no.na$predlm = predict(mlm2.3)

#then use geom_line to add fitted lines based on this new column
site <-
  ggplot(no.na, aes(x = week_nr, y = no.large.fuc.q, colour = site)) +
  geom_point() +
  scale_colour_discrete(
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  geom_line(aes(y = predlm), size = 1, position=position_dodge(width=5)) +
  labs(x = "",
       y = "",
       title = "Predicted values of model per site") +
  theme_minimal()


#arrange plots
plot<-ggarrange(model.plot, site, ncol=1, nrow=2, common.legend = TRUE, legend="right")

## Warning: position_dodge requires non-overlapping x intervals

#making subtitle for top
title <-
  expression(atop(
    bold("large thalli density ~ week since first survey + site(random)"),
    "large thalli density = -0.10(week) + 10.61"
  ))

#annotating figure
annotate_figure(
  plot,
  top = text_grob(title),
  left = text_grob(
    "Large thalli density (thalli per 0.25m^2)",
    color = "black",
    rot = 90
  ),
  bottom = text_grob(
    "Weeks since first survey",
    color = "black",
    rot = 0
      )
)

The lines are perfectly on top of each other so I had to use a larger jitter width than usual. The table (below) shows that there’s no difference amoung sites so maybe I just show the top graph?

2.4 Does small thalli density change over time? sqrt(no.small.fuc.q ~ week_nr + week_nr^2 +((1|site)
Added position dodge to the site trend lines because they were stacked (https://stackoverflow.com/questions/10866047/jitter-geom-line)

#remove rows with NA 
no.na<-all[!is.na(all$no.small.fuc.q), ]

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

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

#model
mlm2.4<-lmer(sqrt.small.fuc.q ~ week_nr + sqr.week_nr + (1|site), data=no.na)

## Warning: Some predictor variables are on very different scales: consider
## rescaling

# Extract the prediction data frame
pred.mm <- ggpredict(mlm2.4, terms = c("week_nr"))  # this gives overall predictions for the model

# Plot the predictions of model 
model.plot <-
  ggplot(pred.mm) +
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(
    aes(
      x = x,
      ymin = predicted - std.error,
      ymax = predicted + std.error
    ),
    fill = "lightgrey",
    alpha = 0.5
  ) +
  labs(x = "",
       y = "",
       title = "Predicted values of model") +
  geom_point(data = no.na,
             aes(x = week_nr, y = no.small.fuc.q, colour = site)) +
  scale_colour_discrete(                                    #key modifiers
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  theme_minimal()



#add predicted values from the model to the dataset
no.na$predlm = predict(mlm2.4)

#then use geom_line to add fitted lines based on this new column
site <-
  ggplot(no.na, aes(x = week_nr, y = no.small.fuc.q, colour = site)) +
  geom_point() +
  scale_colour_discrete(
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  geom_line(aes(y = predlm), size = 1, position=position_dodge(width=0.2)) +   #add jitter
  labs(x = "",
       y = "",
       title = "Predicted values of model per site") +
  theme_minimal()


#arrange plots
plot<-ggarrange(model.plot, site, ncol=1, nrow=2, common.legend = TRUE, legend="right")

#making subtitle for top
title <-
  expression(atop(
    bold("sqrt(small thalli density) ~ week + week^2 + site(random)"),
    "small thalli density = 0.000014(week)^4 - 0.0025(week)^3 + 0.19(week)^2 - 6.77(week) + 105.06",
  ))

#annotating figure
annotate_figure(
  plot,
  top = text_grob(title),
  left = text_grob(
    "Small thalli density (thalli per 0.25m^2)",
    color = "black",
    rot = 90
  ),
  bottom = text_grob(
    "Weeks since first survey",
    color = "black",
    rot = 0
      )
)

Looks weird that the trend line crosses zero because we can’t have a negative density. Also the interecept in the equation is not what I’m seeing in the graph

3. Does reproduction change over time?

3.1 Does oogonia/conceptacle change over time?
oog.per.con ~ week + site(random)

#remove rows with NA in oog.per.con, was causing a problem later on with graphing. This doesn't affect the model since it was ignorning those rows anyway
no.na<-all[!is.na(all$oog.per.con), ]

#model
mlm3.1<-lmer(oog.per.con~ week_nr + (1|site), data=no.na)

# Extract the prediction data frame
pred.mm <- ggpredict(mlm3.1, terms = c("week_nr"))  # this gives overall predictions for the model

# Plot the predictions of model 
model.plot <-
  ggplot(pred.mm) +
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(
    aes(
      x = x,
      ymin = predicted - std.error,
      ymax = predicted + std.error
    ),
    fill = "lightgrey",
    alpha = 0.5
  ) +
  labs(x = "",
       y = "",
       title = "Predicted values of model") +
  geom_point(data = no.na,
             aes(x = week_nr, y = oog.per.con, colour = site)) +
  scale_colour_discrete(                                    #key modifiers
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  theme_minimal()




#add predicted values from the model to the dataset
no.na$predlm = predict(mlm3.1)

#then use geom_line to add fitted lines based on this new column
site <-
  ggplot(no.na, aes(x = week_nr, y = oog.per.con, colour = site)) +
  geom_point() +
  scale_colour_discrete(
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  geom_line(aes(y = predlm), size = 1) +
  labs(x = "",
       y = "",
       title = "Predicted values of model per site") +
  theme_minimal()


#arrange plots
plot<-ggarrange(model.plot, site, ncol=1, nrow=2, common.legend = TRUE, legend="right")

#making subtitle for top
title <-
  expression(atop(
    bold("oogonia per conceptacle ~ week + site(random)"),
    "oogonia per conceptacle = -0.18(week) + 30.12"
  ))

#annotating figure
annotate_figure(
  plot,
  top = text_grob(title),
  left = text_grob(
    "Oogonia per receptacle",
    color = "black",
    rot = 90
  ),
  bottom = text_grob(
    "Weeks since first survey",
    color = "black",
    rot = 0
      )
)

3.2 Does oogonia/thalli change over time?
oog.thalli ~ time + site(random)

#remove rows with NA in oog.per.con, was causing a problem later on with graphing. This doesn't affect the model since it was ignorning those rows anyway
no.na<-all[!is.na(all$oog.thalli), ]

#model
mlm3.2<-lmer(oog.thalli~ week_nr + (1|site), data=no.na)

# Extract the prediction data frame
pred.mm <- ggpredict(mlm3.2, terms = c("week_nr"))  # this gives overall predictions for the model

# Plot the predictions of model 
model.plot <-
  ggplot(pred.mm) +
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(
    aes(
      x = x,
      ymin = predicted - std.error,
      ymax = predicted + std.error
    ),
    fill = "lightgrey",
    alpha = 0.5
  ) +
  labs(x = "",
       y = "",
       title = "Predicted values of model") +
  geom_point(data = no.na,
             aes(x = week_nr, y = oog.thalli, colour = site)) +
  scale_colour_discrete(                                    #key modifiers
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  theme_minimal()




#add predicted values from the model to the dataset
no.na$predlm = predict(mlm3.2)

#then use geom_line to add fitted lines based on this new column
site <-
  ggplot(no.na, aes(x = week_nr, y = oog.thalli, colour = site)) +
  geom_point() +
  scale_colour_discrete(
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  geom_line(aes(y = predlm), size = 1) +
  labs(x = "",
       y = "",
       title = "Predicted values of model per site") +
  theme_minimal()


#arrange plots
plot<-ggarrange(model.plot, site, ncol=1, nrow=2, common.legend = TRUE, legend="right")

#making subtitle for top
title <-
  expression(atop(
    bold("oogonia per thalli ~ week + site(random)"),
    "oogonia per thalli = 1180(week) + 253,623"
  ))

#annotating figure
annotate_figure(
  plot,
  top = text_grob(title),
  left = text_grob(
    "Oogonia per thalli",
    color = "black",
    rot = 90
  ),
  bottom = text_grob(
    "Weeks since first survey",
    color = "black",
    rot = 0
      )
)

Brickyard Park and Paradise cay lines are very close. May need to jitter

3.3 Does conceptical/thalli change over time?
con.thalli ~ time + site(random)

#remove rows with NA 
no.na<-all[!is.na(all$con.thalli), ]

#model
mlm3.3<-lmer(con.thalli~ week_nr + (1|site), data=no.na)

# Extract the prediction data frame
pred.mm <- ggpredict(mlm3.3, terms = c("week_nr"))  # this gives overall predictions for the model

# Plot the predictions of model 
model.plot <-
  ggplot(pred.mm) +
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(
    aes(
      x = x,
      ymin = predicted - std.error,
      ymax = predicted + std.error
    ),
    fill = "lightgrey",
    alpha = 0.5
  ) +
  labs(x = "",
       y = "",
       title = "Predicted values of model") +
  geom_point(data = no.na,
             aes(x = week_nr, y = con.thalli, colour = site)) +
  scale_colour_discrete(                                    #key modifiers
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  theme_minimal()



#add predicted values from the model to the dataset
no.na$predlm = predict(mlm3.3)

#then use geom_line to add fitted lines based on this new column
site <-
  ggplot(no.na, aes(x = week_nr, y = con.thalli, colour = site)) +
  geom_point() +
  scale_colour_discrete(
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  geom_line(aes(y = predlm), size = 1) +
  labs(x = "",
       y = "",
       title = "Predicted values of model per site") +
  theme_minimal()


#arrange plots
plot<-ggarrange(model.plot, site, ncol=1, nrow=2, common.legend = TRUE, legend="right")

#making subtitle for top
title <-
  expression(atop(
    bold("conceptacle per thalli ~ week since first survey + site(random)"),
         "conceptacle per thalli = 52.72(week) + 11494.49"
      ))

#annotating figure
annotate_figure(
  plot,
  top = text_grob(title),
  left = text_grob(
    "Conceptacle per thalli",
    color = "black",
    rot = 90
  ),
  bottom = text_grob(
    "Weeks since first survey",
    color = "black",
    rot = 0
      )
)

Still need to remove outlier. Hard to read this graph

3.4 Does percent reproductive tissue change over time?
perc.rdw ~ time + site(random)

#remove rows with NA 
no.na<-all[!is.na(all$perc.rdw), ]

#model
mlm3.4<-lmer(perc.rdw~ week_nr + (1|site), data=no.na)

# Extract the prediction data frame
pred.mm <- ggpredict(mlm3.4, terms = c("week_nr"))  # this gives overall predictions for the model

# Plot the predictions of model 
model.plot <-
  ggplot(pred.mm) +
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(
    aes(
      x = x,
      ymin = predicted - std.error,
      ymax = predicted + std.error
    ),
    fill = "lightgrey",
    alpha = 0.5
  ) +
  labs(x = "",
       y = "",
       title = "Predicted values of model") +
  geom_point(data = no.na,
             aes(x = week_nr, y = perc.rdw, colour = site)) +
  scale_colour_discrete(                                    #key modifiers
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  theme_minimal()



#add predicted values from the model to the dataset
no.na$predlm = predict(mlm3.4)

#then use geom_line to add fitted lines based on this new column
site <-
  ggplot(no.na, aes(x = week_nr, y = perc.rdw, colour = site)) +
  geom_point() +
  scale_colour_discrete(
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  geom_line(aes(y = predlm), size = 1) +
  labs(x = "",
       y = "",
       title = "Predicted values of model per site") +
  theme_minimal()


#arrange plots
plot<-ggarrange(model.plot, site, ncol=1, nrow=2, common.legend = TRUE, legend="right")

#making subtitle for top
title <-
  expression(atop(
    bold("percent reproductive dry weight ~ week since first survey + site(random)"),
    "percent reproductive  dry weight = 0.023(week) + 18.84"
      ))

#annotating figure
annotate_figure(
  plot,
  top = text_grob(title),
  left = text_grob(
    "Percent reproductive dry weight",
    color = "black",
    rot = 90
  ),
  bottom = text_grob(
    "Weeks since first survey",
    color = "black",
    rot = 0
      )
)

3.5 Does percent reproductive apices change over time?
perc.ra ~ time + site(random)

#remove rows with NA 
no.na<-all[!is.na(all$perc.ra), ]

#model
mlm3.5<-lmer(perc.ra~ week_nr + (1|site), data=no.na)

# Extract the prediction data frame
pred.mm <- ggpredict(mlm3.5, terms = c("week_nr"))  # this gives overall predictions for the model

# Plot the predictions of model 
model.plot <-
  ggplot(pred.mm) +
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(
    aes(
      x = x,
      ymin = predicted - std.error,
      ymax = predicted + std.error
    ),
    fill = "lightgrey",
    alpha = 0.5
  ) +
  labs(x = "",
       y = "",
       title = "Predicted values of model") +
  geom_point(data = no.na,
             aes(x = week_nr, y = perc.ra, colour = site)) +
  scale_colour_discrete(                                    #key modifiers
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  theme_minimal()



#add predicted values from the model to the dataset
no.na$predlm = predict(mlm3.5)

#then use geom_line to add fitted lines based on this new column
site <-
  ggplot(no.na, aes(x = week_nr, y = perc.ra, colour = site)) +
  geom_point() +
  scale_colour_discrete(
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  geom_line(aes(y = predlm), size = 1) +
  labs(x = "",
       y = "",
       title = "Predicted values of model per site") +
  theme_minimal()


#arrange plots
plot<-ggarrange(model.plot, site, ncol=1, nrow=2, common.legend = TRUE, legend="right")

#making subtitle for top
title <-
  expression(atop(
    bold("percent reproductive apices ~ week since first survey + site(random)"),
         "percent reproductive apices = 0.054(week) + 26.34"
      ))

#annotating figure
annotate_figure(
  plot,
  top = text_grob(title),
  left = text_grob(
    "Percent reproductive apices",
    color = "black",
    rot = 90
  ),
  bottom = text_grob(
    "Weeks since first survey",
    color = "black",
    rot = 0
      )
)

3.6 Does conceptacles per receptacle change over time?
no.concept.recp ~ time + site(random)

#remove rows with NA 
no.na<-all[!is.na(all$no.concept.recp), ]

#model
mlm3.6<-lmer(no.concept.recp~ week_nr + (1|site), data=no.na)

# Extract the prediction data frame
pred.mm <- ggpredict(mlm3.6, terms = c("week_nr"))  # this gives overall predictions for the model

# Plot the predictions of model 
model.plot <-
  ggplot(pred.mm) +
  geom_line(aes(x = x, y = predicted)) +
  geom_ribbon(
    aes(
      x = x,
      ymin = predicted - std.error,
      ymax = predicted + std.error
    ),
    fill = "lightgrey",
    alpha = 0.5
  ) +
  labs(x = "",
       y = "",
       title = "Predicted values of model") +
  geom_point(data = no.na,
             aes(x = week_nr, y = no.concept.recp, colour = site)) +
  scale_colour_discrete(                                    #key modifiers
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  theme_minimal()



#add predicted values from the model to the dataset
no.na$predlm = predict(mlm3.6)

#then use geom_line to add fitted lines based on this new column
site <-
  ggplot(no.na, aes(x = week_nr, y = no.concept.recp, colour = site)) +
  geom_point() +
  scale_colour_discrete(
    name = "Site",
    breaks = c("pc.cc", "nd.eos", "by.rb", "hs.fp"),
    labels = c(
      "Paradise Cay",
      "Point Chauncy",
      "Brickyard Park",
      "Horseshoe Bay"
    )
  ) +
  geom_line(aes(y = predlm), size = 1) +
  labs(x = "",
       y = "",
       title = "Predicted values of model per site") +
  theme_minimal()


#arrange plots
plot<-ggarrange(model.plot, site, ncol=1, nrow=2, common.legend = TRUE, legend="right")

#making subtitle for top
title <-
  expression(atop(
    bold("conceptacles per receptacle ~ week since first survey + site(random)"),
    "conceptacles per receptacle = 0.037(week) + 292.29"
  ))

#annotating figure
annotate_figure(
  plot,
  top = text_grob(title),
  left = text_grob(
    "conceptacles per receptacle",
    color = "black",
    rot = 90
  ),
  bottom = text_grob(
    "Weeks since first survey",
    color = "black",
    rot = 0
      )
)

####Table####

Not what we’re going for but good to know https://ourcodingclub.github.io/tutorials/mixed-models/#presenting

library(stargazer)

## Warning: package 'stargazer' was built under R version 4.0.3

## 
## Please cite as:

##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.

##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer

stargazer(mlm2.1, type = "text",
          digits = 3,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digit.separator = "")

## 
## =================================================
##                          Dependent variable:     
##                     -----------------------------
##                                 cover            
## -------------------------------------------------
## week_nr                       -0.331***          
##                                (0.071)           
##                                                  
## Constant                      43.281***          
##                                (3.585)           
##                                                  
## -------------------------------------------------
## Observations                     61              
## Log Likelihood                -236.102           
## Akaike Inf. Crit.              480.204           
## Bayesian Inf. Crit.            488.647           
## =================================================
## Note:               *p<0.05; **p<0.01; ***p<0.001

Not quite what we want but I like that you can put multiple predictors andd models in the same table.
https://cran.r-project.org/web/packages/sjPlot/vignettes/tab_mixed.html

# load required packages
library(sjPlot) # table functions

## Warning: package 'sjPlot' was built under R version 4.0.5

## Learn more about sjPlot with 'browseVignettes("sjPlot")'.

library(sjmisc)

## 
## Attaching package: 'sjmisc'

## The following object is masked from 'package:purrr':
## 
##     is_empty

## The following object is masked from 'package:tidyr':
## 
##     replace_na

## The following object is masked from 'package:tibble':
## 
##     add_case

library(sjlabelled)

## 
## Attaching package: 'sjlabelled'

## The following object is masked from 'package:forcats':
## 
##     as_factor

## The following object is masked from 'package:dplyr':
## 
##     as_label

library(insight)
library(tibble)

#table
tab_model(mlm2.1)

	cover
Predictors	Estimates	CI	p
(Intercept)	43.28	36.26 – 50.31	<0.001
week_nr	-0.33	-0.47 – -0.19	<0.001
Random Effects
σ2	129.06
τ00 site	21.24
ICC	0.14
N site	4
Observations	61
Marginal R2 / Conditional R2	0.236 / 0.344

Only the confidence interval looks different. Why? https://stats.stackexchange.com/questions/146988/getting-degrees-of-freedom-from-lmer
https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html
http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#why-doesnt-lme4-display-denominator-degrees-of-freedomp-values-what-other-options-do-i-have

Why the p values?
https://cran.r-project.org/web/packages/sjPlot/vignettes/tab_mixed.html

List of what we can display in table or here https://strengejacke.github.io/sjPlot/reference/tab_model.html

?tab_model

## starting httpd help server ... done

Expanding on the previous table using: https://strengejacke.github.io/sjPlot/articles/tab_model_estimates.html

tab_model(mlm2.1, show.se = TRUE, show.std = TRUE, show.stat = TRUE)

	cover
Predictors	Estimates	std. Error	std. Beta	standardized std. Error	CI	standardized CI	Statistic	p
(Intercept)	43.28	3.58	0.00	0.20	36.26 – 50.31	-0.38 – 0.39	12.07	<0.001
week_nr	-0.33	0.07	-0.50	0.11	-0.47 – -0.19	-0.70 – -0.29	-4.64	<0.001
Random Effects
σ2	129.06
τ00 site	21.24
ICC	0.14
N site	4
Observations	61
Marginal R2 / Conditional R2	0.236 / 0.344

Multiple models on one table. Abundance (cover, density, small thalli denisty, large thalli density)

tab_model(mlm2.1,  mlm2.2, mlm2.3, mlm2.4, #models to table
  show.ci = 0.95,        # displays 95% confidence interval
  show.se = TRUE,        # displayed standard error (standard deviation also available)
  show.stat = FALSE,     # displays coefficients's test statistic
  emph.p = TRUE,         # bolds significant pvalues
  df.method = "kr",   # pick method for how df is calculted. This way is closer to the model output
  p.style = "numeric_stars", # displays numeric and asterik pvalues
  title = "Linear Mixed Models for Abundance", #main title of table
  linebreak = TRUE,       # supposed to add space between models but I'm not sure it worked
  dv.labels = c("Percent Cover", "sqrt(Total Density)", "Large thalli density", "sqrt(Small thalli density)"),  # change dependent variables (column names)
  pred.labels = c("(Intercept)", "week number", "sqrt(week number)"), #change predictor names
  CSS = css_theme("cells")    #adds lines/cells to table
)

Linear Mixed Models for Abundance

	Percent Cover				sqrt(Total Density)				Large thalli density				sqrt(Small thalli density)
Predictors	Estimates	std. Error	CI	p	Estimates	std. Error	CI	p	Estimates	std. Error	CI	p	Estimates	std. Error	CI	p
(Intercept)	43.28 ***	3.59	35.13 – 51.44	<0.001	8.58 ***	0.73	6.99 – 10.17	<0.001	10.61 ***	1.11	8.36 – 12.87	<0.001	10.25 ***	0.87	8.37 – 12.12	<0.001
week number	-0.33 ***	0.07	-0.47 – -0.19	<0.001	-0.16 ***	0.04	-0.25 – -0.08	<0.001	-0.10 ***	0.03	-0.16 – -0.05	0.001	-0.33 ***	0.05	-0.42 – -0.23	<0.001
sqrt(week number)					0.00 *	0.00	0.00 – 0.00	0.013					0.00 ***	0.00	0.00 – 0.01	<0.001
Random Effects
σ2	129.06				3.06				16.05				2.53
τ00	21.24 site				0.85 site				0.00 site				1.22 site
ICC	0.14				0.22								0.33
N	4 site				4 site				4 site				4 site
Observations	61				61				57				57
Marginal R2 / Conditional R2	0.236 / 0.344				0.326 / 0.472				0.196 / NA				0.474 / 0.646
p<0.05   ** p<0.01   *** p<0.001

Not sure these are the best labels

tab_model df.method and p.val let’s you choose the computing method:
“Method for computing degrees of freedom for p-values, standard errors and confidence intervals (CI). Only applies to mixed models. Use df.method =”wald" for a faster, but less precise computation. df.method = “kenward” (or df.method = “kr”) uses Kenward-Roger approximation for the degrees of freedom. df.method = “satterthwaite” uses Satterthwaite’s approximation and “ml1” uses a “m-l-1” heuristic see degrees_of_freedom for details). Use show.df = TRUE to show the approximated degrees of freedom for each coefficient."

auto.label

Reproduction table

#simple table
tab_model(mlm3.1, mlm3.2, mlm3.3, mlm3.4, mlm3.5, mlm3.6)

	oog.per.con			oog.thalli			con.thalli			perc.rdw			perc.ra			no.concept.recp
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	30.12	20.32 – 39.91	<0.001	253622.53	42437.20 – 464807.86	0.019	11494.49	4025.71 – 18963.27	0.003	18.84	8.79 – 28.88	<0.001	26.34	13.64 – 39.05	<0.001	292.29	215.13 – 369.46	<0.001
week_nr	-0.18	-0.39 – 0.02	0.080	-1180.50	-5595.32 – 3234.31	0.600	52.72	-114.15 – 219.59	0.536	0.02	-0.14 – 0.19	0.786	0.05	-0.16 – 0.27	0.619	0.04	-1.29 – 1.36	0.956
Random Effects
σ2	231.85			105801613226.78			151173665.91			150.26			243.93			9529.99
τ00	42.59 site			20306261789.71 site			20754818.46 site			67.94 site			107.81 site			3844.48 site
ICC	0.16			0.16			0.12			0.31			0.31			0.29
N	4 site			4 site			4 site			4 site			4 site			4 site
Observations	57			57			57			57			57			57
Marginal R2 / Conditional R2	0.044 / 0.193			0.004 / 0.164			0.006 / 0.126			0.001 / 0.312			0.003 / 0.309			0.000 / 0.287

#modified table         
tab_model(mlm3.1, mlm3.2, mlm3.3, mlm3.4, mlm3.5, mlm3.6,
  show.ci = 0.95,
  show.se = TRUE,
  show.stat = FALSE,
  emph.p = TRUE,
  df.method = "kr",
  p.style = "numeric_stars",
  title = "Linear Mixed Models for Reproduction",
  linebreak = TRUE,
  dv.labels = c(
    "oogonia/conceptacle",
    "oogonia/thalli",
    "conceptacle/thalli",
    "percent reproductive dry weight",
    "percent reproductive apices",
    "conceptacle/receptacle"
      ),
  pred.labels = c("(Intercept)", "week number", "sqrt(week number)"),
  CSS = css_theme("cells")
  )

## Length of `pred.labels` does not equal number of predictors, no labelling applied.

Linear Mixed Models for Reproduction

	oogonia/conceptacle				oogonia/thalli				conceptacle/thalli				percent reproductive dry weight				percent reproductive apices				conceptacle/receptacle
Predictors	Estimates	std. Error	CI	p	Estimates	std. Error	CI	p	Estimates	std. Error	CI	p	Estimates	std. Error	CI	p	Estimates	std. Error	CI	p	Estimates	std. Error	CI	p
(Intercept)	30.12 ***	5.00	18.68 – 41.56	<0.001	253622.53 *	107796.04	5941.38 – 501303.68	0.046	11494.49 *	3813.02	2978.54 – 20010.44	0.013	18.84 *	5.13	5.90 – 31.77	0.013	26.34 **	6.48	10.03 – 42.66	0.008	292.29 ***	39.38	194.28 – 390.31	<0.001
week_nr	-0.18	0.11	-0.40 – 0.03	0.086	-1180.50	2253.11	-5701.61 – 3340.61	0.603	52.72	85.17	-118.18 – 223.62	0.539	0.02	0.08	-0.15 – 0.19	0.787	0.05	0.11	-0.16 – 0.27	0.621	0.04	0.68	-1.32 – 1.39	0.956
Random Effects
σ2	231.85				105801613226.78				151173665.91				150.26				243.93				9529.99
τ00	42.59 site				20306261789.71 site				20754818.46 site				67.94 site				107.81 site				3844.48 site
ICC	0.16				0.16				0.12				0.31				0.31				0.29
N	4 site				4 site				4 site				4 site				4 site				4 site
Observations	57				57				57				57				57				57
Marginal R2 / Conditional R2	0.044 / 0.193				0.004 / 0.164				0.006 / 0.126				0.001 / 0.312				0.003 / 0.309				0.000 / 0.287
p<0.05   ** p<0.01   *** p<0.001