Wong_whelk/Knobbed_Mixing_Models.R at main · nmotz44/Wong_whelk · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
##### mixing models for knobbed whelk #####
# first run of mixing models for knobbed whelk
# loop through multiple TDFs from 0-3 by 0.5

## load data and packages
# packages
library(tidyverse)
library(simmr)
source('format_simmr_data.R')

# load data
prey = read.csv("Data/kmeans_output.csv") |>
  select(-X.1, -CommonName, -meand13C, -sdd13C, -meand15N, -sdd15N)
whelk_data = read.csv("Data/SIA_Database_2025_12_11.csv")

# rename CommonNameTissue to CommonName
prey = prey %>%
  rename("CommonName" = "CommonNameTissue") %>%
  mutate(kmeans = case_when(kmeans == 1 ~ "one",
                            kmeans == 2 ~ "two",
                            kmeans == 3 ~ "three",
                            kmeans == 4 ~ "four",
                            TRUE ~ "other"))

# subset for only knobbed whelk
whelk_data = whelk_data %>%
  filter(ProjID == "Whelk" & CommonName == "Knobbed Whelk")

# add in a k means column
whelk_data$kmeans = "Knobbed Whelk"

# put knobbed whelk in the dataframe
iso_data = rbind(prey, whelk_data)

# make kmeans categorical
iso_data$kmeans = as.character(iso_data$kmeans)

##### format mixng models ######
# make nitrogen TDF
n_tdf = c(0, 0.5, 1, 1.5, 2, 2.5, 3)

# make blank model fits
model_fits = data_frame(tdf = n_tdf, DIC = NA)

for (i in 1:length(n_tdf)) {
  # print TDF for this iteration
  print("########################################################")
  print(paste0("TDF for this run: ", n_tdf[i]))

  # format data for simmr
  simmr_object = format_simmr_data(data = iso_data,
                                   consumer.names = "Knobbed Whelk",
                                   prey.names = c("one", "two", "three", "four"),
                                   species.column = "kmeans",
                                   d13C.column = "d13C",
                                   d15N.column = "d15N",
                                   disc.factors.c = 1.35,
                                   disc.factors.c.sd = 0.3,
                                   disc.factors.n = n_tdf[i],
                                   disc.factors.n.sd = 0.5,
                                   data.type = "raw"
                                   )

  # run simmr model
  simmr.out = simmr_mcmc(simmr_object, mcmc_control = list(iter = 10000, burn = 1000, thin = 10, n.chain = 4))

  # check diagnostics
  summary(simmr.out, type = "diagnostics")

  # make bi-plot
  plot(simmr_object, title = paste0(n_tdf[i], " Bi-plot"))

  # make boxplot
  plot(simmr.out, type = "boxplot", title = paste0(n_tdf[i], " Box Plot"))

  # density plot
  plot(simmr.out, type = "density", title = paste0(n_tdf[i], " Box Plot"))

  # correlation plot
  plot(simmr.out, type = "matrix", title = paste0(n_tdf[i], " Box Plot"))


  # add model fit into dataframe
  model_fits$DIC[model_fits$tdf == n_tdf[i]] = simmr.out$output$`Knobbed Whelk`$BUGSoutput$DIC

  # add in correlation plot

  # print
  print("#####################################################################")
}


###### rerun models without group 4 #####
for (i in 1:length(n_tdf)) {
  # print TDF for this iteration
  print("########################################################")
  print(paste0("TDF for this run: ", n_tdf[i]))

  # format data for simmr
  simmr_object = format_simmr_data(data = iso_data,
                                   consumer.names = "Knobbed Whelk",
                                   prey.names = c("one", "two", "three"),
                                   species.column = "kmeans",
                                   d13C.column = "d13C",
                                   d15N.column = "d15N",
                                   disc.factors.c = 1.35,
                                   disc.factors.c.sd = 0.3,
                                   disc.factors.n = n_tdf[i],
                                   disc.factors.n.sd = 0.5,
                                   data.type = "raw"
  )

  # run simmr model
  simmr.out = simmr_mcmc(simmr_object, mcmc_control = list(iter = 10000, burn = 1000, thin = 10, n.chain = 4))

  # check diagnostics
  summary(simmr.out, type = "diagnostics")

  # make bi-plot
  plot(simmr_object, title = paste0(n_tdf[i], " Bi-plot"))

  # make boxplot
  plot(simmr.out, type = "boxplot", title = paste0(n_tdf[i], " Box Plot"))

  # density plot
  plot(simmr.out, type = "density", title = paste0(n_tdf[i], " Box Plot"))

  # correlation plot
  plot(simmr.out, type = "matrix", title = paste0(n_tdf[i], " Box Plot"))


  # add model fit into dataframe
  model_fits$DIC[model_fits$tdf == n_tdf[i]] = simmr.out$output$`Knobbed Whelk`$BUGSoutput$DIC

  # add in correlation plot

  # print
  print("#####################################################################")
}

# models were slightly better with group 4 included