nathanmct / Beaufort-Sea-Lagoons

Supplementary Materials for Harris et al. (2018)
1 stars 0 forks source link

Issues to obtain 95% credible interval #1

Open GuillaumeBridier opened 4 years ago

GuillaumeBridier commented 4 years ago

Dear Nathan,

Thank you so much to share your R codes based on the analyses performed in Harris et al. 2018. I really liked your paper and your approach on stable isotope analyses. I would like also to make some bayesian mixing models on my own macrofauna data (from NE Greenland) and make a figure similar as your fig.5. I first tried to run our r script with your data. Everything was running well until the part "Visually asses convergence of chains and make note of Gelman diagnostics". Just below this line of code, I faced several issues:

I cannot plot convergence of chains and calculate Gelman diagnostics. R return me this error messages : > plot(susp_out, type = 'convergence') Error in match.arg(type, several.ok = TRUE) : 'arg' should be one of “isospace”, “histogram”, “density”, “matrix”, “boxplot” > summary(susp_out, type = "diagnostics") _Summary for 1 Error in apply(outall, 2, "quantile", probs = c(0.025, 0.25, 0.5, 0.75, : dim(X) doit avoir un longueur positive [dim(X) should have a positive length]

I cannot plot OM source correlations: > plot(susp_out, type='matrix', group=1:max(grp.susp)) _Error in colnames<-(*tmp*, value = x$input$sourcenames) : attempt to set 'colnames' on an object with less than two dimensions

And more importantly, I cannot obtain 95% credible intervals. The summary.simmr_output function was first not find by R (it's a bit weird as I supposed that it should be included inside the simmr package...). But I loaded the code of this function from this page: https://rdrr.io/cran/simmr/src/R/summary.simmr_output.R I then tried to run again the code but R return me this error message: > susp.summ <- summary.simmr_output(susp_out, type=c("quantiles"), group=c(1:max(susp$Code))) _Error in summary.simmr_output(suspout, type = c("quantiles"), group = c(1:max(susp$Code))) : l'argument 3 correspond à plusieurs arguments formels [arguments 3 corresponds to several formal arguments]

Do you have any idea on what's wrong with my analysis? Thank you very much! :)

All the best, Guillaume

nathanmct commented 4 years ago

Hi Guillaume,

Thanks so much for the correspondence! I'm glad the paper and code are helpful to you!

Unfortunately I'm not exactly sure what is causing the problem, but I have a series of ideas that you could try: -make sure you have the most recent version of the package installed. Since the 'summary.simmr_output' function wasn't found, it makes me think you have an old or partial version installed. install.packages('simmr') should do the trick. -make sure you have JAGS installed. It runs in the background and does the Bayesian computing. I downloaded it from: https://sourceforge.net/projects/mcmc-jags/files/ . It doesn't need to be open, just installed. -do you have the 'simmrmcmc.r' file in your working directory? This is the function that does the computing in the step "susp_out = simmr_mcmc(simmr_groups_susp)". It sounds like the errors are saying the objects you are trying to plot() and summary() are the wrong length, which could be because they weren't created correctly by the "simmr_mcmc" step. -you can also follow the vignette for the package: https://cran.r-project.org/web/packages/simmr/vignettes/simmr.html The package creator Andrew Parnell has more info on his github site too. -Can you download the data for the Harris et al. paper and run the script on that? That will determine if it's the code that's not working, or if your own input data is not structured correctly. Make sure yours is the exact same as mine. If it's different, say in the number of columns, then you'll have to through the code and make sure the correct columns are called.

Let me know if any of this helps. If it doesn't we can go back to the drawing board.

Best regards, Nathan

On Tue, May 26, 2020 at 11:50 AM GuillaumeBridier notifications@github.com wrote:

Dear Nathan,

Thank you so much to share your R codes based on the analyses performed in Harris et al. 2018. I really liked your paper and your approach on stable isotope analyses. I would like also to make some bayesian mixing models on my own macrofauna data (from NE Greenland) and make a figure similar as your fig.5. I first tried to run our r script with your data. Everything was running well until the part "Visually asses convergence of chains and make note of Gelman diagnostics". Just below this line of code, I faced several issues:

I cannot plot convergence of chains and calculate Gelman diagnostics. R return me this error messages : > plot(susp_out, type = 'convergence')

Error in match.arg(type, several.ok = TRUE) : 'arg' should be one of “isospace”, “histogram”, “density”, “matrix”, “boxplot” > summary(susp_out, type = "diagnostics")

Summary for 1 Error in apply(out_all, 2, "quantile", probs = c(0.025, 0.25, 0.5, 0.75, : dim(X) doit avoir un longueur positive [dim(X) should have a positive length]

I cannot plot OM source correlations: > plot(susp_out, type='matrix', group=1:max(grp.susp))

Error in colnames<-(tmp, value = x$input$source_names) : attempt to set 'colnames' on an object with less than two dimensions

And more importantly, I cannot obtain 95% credible intervals. The summary.simmr_output function was first not find by R (it's a bit weird as I supposed that it should be included inside the simmr package...). But I loaded the code of this function from this page: https://rdrr.io/cran/simmr/src/R/summary.simmr_output.R I then tried to run again the code but R return me this error message: > susp.summ <- summary.simmr_output(susp_out, type=c("quantiles"), group=c(1:max(susp$Code)))

Error in summary.simmr_output(susp_out, type = c("quantiles"), group = c(1:max(susp$Code))) : l'argument 3 correspond à plusieurs arguments formels [arguments 3 corresponds to several formal arguments]

Do you have any idea on what's wrong with my analysis? Thank you very much! :)

All the best, Guillaume

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/nathanmct/Beaufort-Sea-Lagoons/issues/1, or unsubscribe https://github.com/notifications/unsubscribe-auth/AD2CTMJZN6JQB2UUZO3AL3LRTQMPHANCNFSM4NLATCFQ .

GuillaumeBridier commented 4 years ago

Hello Nathan,

Thank you so much for your quick reply! :)

I tried yesteday to run the script on the data from Harris et al. 2018 and defined your Beaufort-Sea-Lagoons file as my working directory (with all R and csv files). However, I didn't re-install all the packages (almost all packages were already on my computer (e.g. siar, SIBER, rjags...). JAGS was also already installed on my computer.

I tried today to run again the script with the new R version (R4.0.0) and I re-install all the package and the JAGS software from your URL: https://sourceforge.net/projects/mcmc-jags/files/ It's really strange as now I face new issues:

In the "SIBER Table" section, I cannot compare the probability that g1 isopace is greater than g2 isopace > comp_SEAb(SEA.Jags.Trophic$SEA.B, iso$Feeding.Type) _Error in combn(seqalong(s), 2) : n < m As a consequence, all the other lines of code from this section are not running...

In the "Need similar table with statistics for all species..." section, I face also some other issues met2 <- data.frame(SEA.Jags.iso$SEA) names(met2) <- levels(iso.SIBER.genus$Genus) met2 <- melt(met2) No id variables; using all as measure variables > names(met2) <- c("Genus", "value") > hh <- function(x){ h <- hdr(x$value,h = bw.nrd0(x$value)) data.frame(mean = mean(x$value), sd=sd(x$value), mode = h$mode,lo95 = h$hdr[2,1], hi95 = h$hdr[2,2]) } > tSEAb <- ddply(met2,.(Genus), hh) Error in hdr(x$value, h = bw.nrd0(x$value)) : Insufficient data

And finally, in the "simmr (stable isotope mixing model in R)" section, I cannot create objects for each trophic guilds... (although it was appararently running yesterday... ^^) simmr_groups_susp <- simmr_load(mixtures=mix.susp, source_names=s_names, source_means=s_means, source_sds=s_sds, correction_means=c_means, correction_sds=c_sds, group=grp.susp) _Error in simmr_load(mixtures = mix.depo, source_names = s_names, source_means = s_means, : sourcenames must be a vector

It's really strange... Have you ever faced this issue before? Maybe it's still linked to my new package versions?

Best regards, Guillaume

nathanmct commented 4 years ago

Hi Guillaume,

This is quite peculiar. I will open an R environment and re-run all the code as-is from GitHub to try an replicate the problem. I'll get back to you soon.

Best, Nathan

On Wed, May 27, 2020 at 1:42 AM GuillaumeBridier notifications@github.com wrote:

Hello Nathan,

Thank you so much for your quick reply! :)

I tried yesteday to run the script on the data from Harris et al. 2018 and defined your Beaufort-Sea-Lagoons file as my working directory (with all R and csv files). However, I didn't re-install all the packages (almost all packages were already on my computer (e.g. siar, SIBER, rjags...). JAGS was also already installed on my computer.

I tried today to run again the script with the new R version (R4.0.0) and I re-install all the package and the JAGS software from your URL: https://sourceforge.net/projects/mcmc-jags/files/ It's really strange as now I face new issues:

In the "SIBER Table" section, I cannot compare the probability that g1 isopace is greater than g2 isopace > comp_SEAb(SEA.Jags.Trophic$SEA.B, iso$Feeding.Type) Error in combn(seq_along(s), 2) : n < m As a consequence, all the other lines of code from this section are not running...

In the "Need similar table with statistics for all species..." section, I face also some other issues

met2 <- data.frame(SEA.Jags.iso$SEA) names(met2) <- levels(iso.SIBER.genus$Genus) met2 <- melt(met2) No id variables; using all as measure variables > names(met2) <- c("Genus", "value") > hh <- function(x){ h <- hdr(x$value,h = bw.nrd0(x$value)) data.frame(mean = mean(x$value), sd=sd(x$value), mode = h$mode,lo95 = h$hdr[2,1], hi95 = h$hdr[2,2]) } > tSEAb <- ddply(met2,.(Genus), hh) Error in hdr(x$value, h = bw.nrd0(x$value)) : Insufficient data

And finally, in the "simmr (stable isotope mixing model in R)" section, I cannot create objects for each trophic guilds... (although it was appararently running yesterday... ^^)

simmr_groups_susp <- simmr_load(mixtures=mix.susp, source_names=s_names, source_means=s_means, source_sds=s_sds, correction_means=c_means, correction_sds=c_sds, group=grp.susp)

Error in simmr_load(mixtures = mix.depo, source_names = s_names, source_means = s_means, : source_names must be a vector

It's really strange... Have you ever faced this issue before? Maybe it's still linked to my new package versions?

Best regards, Guillaume

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/nathanmct/Beaufort-Sea-Lagoons/issues/1#issuecomment-634549819, or unsubscribe https://github.com/notifications/unsubscribe-auth/AD2CTMO7K73Y6WUAQMDKNWTRTTOA7ANCNFSM4NLATCFQ .

nathanmct commented 4 years ago

Hi Guillaume,

I've replicated the error you describe and unfortunately I am at a loss for why this is occurring. I've emailed Andrew Parnell who created simmr to query the loss of functionality. Thanks for your patience.

I noticed that the simmr_out objects look correctly made, so you could write your own summary functions to get the quantiles. I'll let you know when Andrew gets back to me.

Best, Nathan

On Wed, May 27, 2020 at 10:35 AM Nathan McTigue mctigue@utexas.edu wrote:

Hi Guillaume,

This is quite peculiar. I will open an R environment and re-run all the code as-is from GitHub to try an replicate the problem. I'll get back to you soon.

Best, Nathan

On Wed, May 27, 2020 at 1:42 AM GuillaumeBridier notifications@github.com wrote:

Hello Nathan,

Thank you so much for your quick reply! :)

I tried yesteday to run the script on the data from Harris et al. 2018 and defined your Beaufort-Sea-Lagoons file as my working directory (with all R and csv files). However, I didn't re-install all the packages (almost all packages were already on my computer (e.g. siar, SIBER, rjags...). JAGS was also already installed on my computer.

I tried today to run again the script with the new R version (R4.0.0) and I re-install all the package and the JAGS software from your URL: https://sourceforge.net/projects/mcmc-jags/files/ It's really strange as now I face new issues:

In the "SIBER Table" section, I cannot compare the probability that g1 isopace is greater than g2 isopace > comp_SEAb(SEA.Jags.Trophic$SEA.B, iso$Feeding.Type) Error in combn(seq_along(s), 2) : n < m As a consequence, all the other lines of code from this section are not running...

In the "Need similar table with statistics for all species..." section, I face also some other issues

met2 <- data.frame(SEA.Jags.iso$SEA) names(met2) <- levels(iso.SIBER.genus$Genus) met2 <- melt(met2) No id variables; using all as measure variables > names(met2) <- c("Genus", "value") > hh <- function(x){ h <- hdr(x$value,h = bw.nrd0(x$value)) data.frame(mean = mean(x$value), sd=sd(x$value), mode = h$mode,lo95 = h$hdr[2,1], hi95 = h$hdr[2,2]) } > tSEAb <- ddply(met2,.(Genus), hh) Error in hdr(x$value, h = bw.nrd0(x$value)) : Insufficient data

And finally, in the "simmr (stable isotope mixing model in R)" section, I cannot create objects for each trophic guilds... (although it was appararently running yesterday... ^^)

simmr_groups_susp <- simmr_load(mixtures=mix.susp, source_names=s_names, source_means=s_means, source_sds=s_sds, correction_means=c_means, correction_sds=c_sds, group=grp.susp)

Error in simmr_load(mixtures = mix.depo, source_names = s_names, source_means = s_means, : source_names must be a vector

It's really strange... Have you ever faced this issue before? Maybe it's still linked to my new package versions?

Best regards, Guillaume

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/nathanmct/Beaufort-Sea-Lagoons/issues/1#issuecomment-634549819, or unsubscribe https://github.com/notifications/unsubscribe-auth/AD2CTMO7K73Y6WUAQMDKNWTRTTOA7ANCNFSM4NLATCFQ .

GuillaumeBridier commented 4 years ago

Hi Nathan,

Thank you so much for trying to fix this issue! At least I know that I'm not the only one with this problem ;)

I stay tuned to the chat and look forward to Andrew's answers.

Thank you very much for your help! Best, Guillaume

nathanmct commented 4 years ago

Hi Guillaume,

I heard back from Andrew Parnell about simmr updates. The error is that my code calls the simmrmcmc.R function before running susp_out = simmr_mcmc(simmr_groups_susp). If you don't use that outside function but instead use the internal simmr_mcmc() function that's built into the simmr package, the object is then made correctly for the subsequent code to work. As far as the summary.simmr_output() function, that has been removed and you only need to use summary() now.

Were you planning on using the ggplot code too? The new version of simmr makes the list objects in slightly different order now, so the for-loops that grab the quantiles to create the boxplots are off by a row. I will update the code from Harris et al. on GitHub soon, but wanted to let you know that there are a few things that need to be tweaked before it will all work again.

I'll let you know when I've updated it.

Cheers, Nathan

On Wed, May 27, 2020 at 10:47 PM GuillaumeBridier notifications@github.com wrote:

Hi Nathan,

Thank you so much for trying to fix this issue! At least I know that I'm not the only one with this problem ;)

I stay tuned to the chat and look forward to Andrew's answers.

Thank you very much for your help! Best, Guillaume

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/nathanmct/Beaufort-Sea-Lagoons/issues/1#issuecomment-635141906, or unsubscribe https://github.com/notifications/unsubscribe-auth/AD2CTMK7ZM4VHX5HQYJWZ6LRTYCI3ANCNFSM4NLATCFQ .

nathanmct commented 4 years ago

Hi Guillaume,

Actually, I just edited the code in Beaufort_lagoons_simmr.R for compatibility with the new version of simmr. Go ahead and try it out and let me know if it works. Let me know if it works for you now, and I will close this issue.

Cheers, Nathan

GuillaumeBridier commented 4 years ago

Hi Nathan,

Thank you so much for your reply and for yours code edits! I just tried your script and everything is working fine now! :)

I would like to use your lines of codes to run the same mixing model on my own data for publication and present the results as you did in figure 5. Would you agree with this? Off course, I will cite the URL of your github R codes in the Mat & Met part and mention your help in the acknowledgment part! :)

Thank you again for your help and availability! Cheers, Guillaume

nathanmct commented 4 years ago

Hi Guillaume,

I'm glad we worked out the issue! Knowing you are using the code I generated so you don't have to re-create the wheel makes the whole process worthwhile. This is science! Open data, open code, open source!

Yes, I like Figure 5 from the Harris et al paper. Please mirror the analysis and visualization to whatever degree you would like. And please cite the Harris et al paper. Also check out another paper where I used simmr - McTigue & Dunton 2017 (attached; if the attachment doesn't work through this GitHub chain then I can send it to your private email address; let me know). There is no Figure 5 in that pub, but maybe some of the other food web analysis could be helpful to you.

I will soon be conducting a pan-Arctic stable isotope synthesis project. If you are willing, I would like to include your work once it is published. It will focus mostly on primary producers and grazers, but might include some well represented higher trophic organisms as well. Please keep me posted on your progress (mctigue@utexas.edu).

All the best, Nathan

On Sun, Jun 14, 2020 at 12:53 AM GuillaumeBridier notifications@github.com wrote:

Hi Nathan,

Thank you so much for your reply and for yours code edits! I just tried your script and everything is working fine now! :)

I would like to use your lines of codes to run the same mixing model on my own data for publication and present the results as you did in figure 5. Would you agree with this? Off course, I will cite the URL of your github R codes in the Mat & Met part and mention your help in the acknowledgment part! :)

Thank you again for your help and availability! Cheers, Guillaume

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/nathanmct/Beaufort-Sea-Lagoons/issues/1#issuecomment-643737962, or unsubscribe https://github.com/notifications/unsubscribe-auth/AD2CTMKTRBMCM5NHWTVJUC3RWSFXJANCNFSM4NLATCFQ .

GuillaumeBridier commented 4 years ago

Hi Nathan,

Thank you for your agreement for using R codes. Yes of course, I will cite also your paper Harris et al. 2018. I have dowload the McTigue & Dunton (2017) paper on the Chukchi sea! I will have a look on it ;)

Yes of course! I will be glad if my data may be used for another project! I will let you know about the outcome of my paper! If you ever need to contact me, don't hesitate to send me a message on my email address (guillaume.bridier@univ-brest.fr).

Again many thanks for your help! All the best,

Guillaume