SEA.B values

[BrianHayden]

Hi Andrew,

Is it possible to quickly call the mean/mode and CI values of SEA.B for a particular model? Similar to siarhdrs() in siar.

Cheers,

B

[AndrewLJackson]

not quickly. you need to call

hdrcde::hdr()

on each column of the matrix or data frame that holds the SEA.B values

i need to add a function to do this to the new SIBER package. a couple of people have been asking

On 8 September 2015 at 19:11, Brian Hayden notifications@github.com wrote:

Hi Andrew,

Is it possible to quickly call the mean/mode and CI values of SEA.B for a particular model? Similar to siarhdrs() in siar.

Cheers,

B

— Reply to this email directly or view it on GitHub https://github.com/AndrewLJackson/SIBER/issues/9.

*You can now follow SIAR on Facebook http://www.facebook.com/pages/SIAR-Stable-Isotope-Analysis-in-R/148501811896914 or try my video podcasts http://www.tcd.ie/Zoology/research/research/theoretical/Rpodcasts.php for R

Dr Andrew Jackson Assistant Professor School of Natural Sciences Director of the Trinity Centre for Biodiversity Research Zoology Building, Trinity College Dublin, Dublin 2, Ireland Tel. + 353 1 896 2728, Fax. + 353 1 677 8094, Email. a.jackson@tcd.ie, Twitter: @yodacomplex https://twitter.com/yodacomplex http://www.tcd.ie/Zoology/research/research/theoretical/AndrewJackson.php

[BrianHayden]

Nice one, thankfully I don't have too many groups this time ;)

B

[AndrewLJackson]

hi again - I'm being more polite this time!.. sorry about that abrupt last reply :)

if the SEA.B data is in a matrix called for example x, here x is just random data.... then

x <- matrix(rnorm(200),100,2)

when calling lapply, force x to be a dataframe

lapply(data.frame(x),hdrcde::hdr)

returns a list of density region estimates, with one entry for each column

its not exactly pretty, but you can read off the credible intervals from the output easily enough if thats all you want. Or you could traverse this output list and extract the bits you want, but then its probably as easy to loop over the columns separately and extract it then.

cheers andrew

On 8 September 2015 at 19:49, Brian Hayden notifications@github.com wrote:

Nice one, thankfully I don't have too many groups this time ;)

B

— Reply to this email directly or view it on GitHub https://github.com/AndrewLJackson/SIBER/issues/9#issuecomment-138666899.

*You can now follow SIAR on Facebook http://www.facebook.com/pages/SIAR-Stable-Isotope-Analysis-in-R/148501811896914 or try my video podcasts http://www.tcd.ie/Zoology/research/research/theoretical/Rpodcasts.php for R

Dr Andrew Jackson Assistant Professor School of Natural Sciences Director of the Trinity Centre for Biodiversity Research Zoology Building, Trinity College Dublin, Dublin 2, Ireland Tel. + 353 1 896 2728, Fax. + 353 1 677 8094, Email. a.jackson@tcd.ie, Twitter: @yodacomplex https://twitter.com/yodacomplex http://www.tcd.ie/Zoology/research/research/theoretical/AndrewJackson.php

[AndrewLJackson]

need to implement something like siarhdrs()

[TorKemp]

Hi Andrew,

Is there any update on either assigning names to the matrix output from SiberEllipses() or better still calling the mode and CI values of SEA.B for a unique group/community?

I have found a way to assign names manually to the columns of the matrix, but it is clunky and open to error.

Thanks, Tor

[AndrewLJackson]

hi tor

in short no, and its not high up my list of things to do, sorry. That said, you should be able to rename it from the names stored in the siber object with something like:

BEGIN CODE

data("demo.siber.data") #

create the siber object

siber.example <- createSiberObject(demo.siber.data)

Calculate summary statistics for each group: TA, SEA and SEAc

group.ML <- groupMetricsML(siber.example)

options for running jags

parms <- list() parms$n.iter <- 2 10^4 # number of iterations to run the model for parms$n.burnin <- 1 10^3 # discard the first set of values parms$n.thin <- 10 # thin the posterior by this many parms$n.chains <- 2 # run this many chains

define the priors

priors <- list() priors$R <- 1 * diag(2) priors$k <- 2 priors$tau.mu <- 1.0E-3

fit the ellipses which uses an Inverse Wishart prior

on the covariance matrix Sigma, and a vague normal prior on the

means. Fitting is via the JAGS method.

ellipses.posterior <- siberMVN(siber.example, parms, priors)

The posterior estimates of the ellipses for each group can be used to

calculate the SEA.B for each group.

SEA.B <- siberEllipses(ellipses.posterior)

ccc <- paste(siber.example$all.communities, unlist(siber.example$group.names), sep = "|”)

colnames(SEA.B) <- cc

head(SEA.B)

END CODE

best wishes

andrew

-- –––––––––––– Dr Andrew Jackson, PhD, FTCD Associate Professor School of Natural Sciences, Department of Zoology Trinity College Dublin, the University of Dublin Dublin 2, Ireland.

+353 1 896 2728 | a.jackson@tcd.iemailto:a.jackson@tcd.ie Twitter: @yodacomplexhttps://twitter.com/yodacomplex http://www.tcd.ie/Zoology/research/research/theoretical/AndrewJackson.php

Trinity College Dublin, the University of Dublin is ranked 1st in Ireland and in the top 100 world universities by the QS World University Rankings.

On 26 Oct 2017, at 10:12, TorKemp notifications@github.com<mailto:notifications@github.com> wrote:

Hi Andrew,

Is there any update on either assigning names to the matrix output from SiberEllipses() or better still calling the mode and CI values of SEA.B for a unique group/community?

I have found a way to assign names manually to the columns of the matrix, but it is clunky and open to error.

Thanks, Tor

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/AndrewLJackson/SIBER/issues/9#issuecomment-339602893, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ADK1sEbVLQrgLFQ1L5auSkzJMjtUdgaTks5swEzxgaJpZM4F5trw.

[TorKemp]

Brilliant, thank you Andrew!

[TorKemp]

Hi Andrew, The code above appears to blindly repeat the community levels (2 in the demo data) along the length of the unlisted group names (6 in demo data).
"1|1" "2|2" "1|3" "2|1" "1|2" "2|3"

I think actually the code needs to acknowledge the listing structure, as below, when assigning the communities:

Siber.example$group.names

[1] 1 2 3 Levels: 1 2 3

[[2]] [1] 1 2 3 Levels: 1 2 3

Because of constraint of minimum of 5 individuals in a group to construct an ellipse, I have removed some groups and therefore have communities with different assemblages. I think this is heightening the problem I am having with naming.

But would the following code give the correct order of names as you can see:

cc<-names(ellipses.posterior) colnames(SEA.B) <- cc

Thanks for any thoughts on this! Tor

[costavale]

Hi all, I was trying to extract mode and CI from the SEA.B matrix using the code that Andrew provided and using the suggestion of Tor to keep the name of the Community.Group. Everything seems to work except that when I call the matrix with the hdr values, I have multiple columns (2 minimum, named [,1], [,2]) for the same Community.Group$hdr. I don't understand what is the meaning of the different columns, any suggestions? did I miss some math beyond that (probably)?

Thanks Valentina

[AndrewLJackson]

hi Valentina

im not quite sure what your code is, but i suspect your are calling hdrcde::hdr() which will return a list object with several matrices or vectors

for example

plot(density(faithful$eruptions))

hdrcde::hdr(faithful$eruptions)

will return a list of length = 3, with a matrix of hdr values. There are two values for each of the intervals because it is multimodal (see the figure). It returns the global mode by default (one value) but you can ask for all modes., and it returns the vector of alphas for the corresponding intervals.

does that explain why you are getting more than one value? An interval is defined by two values, and if you have a multimodal distribution then you may get four (or more) values.

best wishes

andrew

–––––––––––– Dr Andrew Jackson, PhD, FTCD Associate Professor School of Natural Sciences, Department of Zoology Trinity College Dublin, the University of Dublin Dublin 2, Ireland.

+353 1 896 2728 | a.jackson@tcd.iemailto:a.jackson@tcd.ie Twitter: @yodacomplexhttps://twitter.com/yodacomplex http://www.tcd.ie/Zoology/research/http://www.tcd.ie/Zoology/research/groups/jackson/groups/jackson/

Trinity College Dublin, the University of Dublin is ranked 1st in Ireland and in the top 100 world universities by the QS World University Rankings.

On 17 Apr 2018, at 14:54, Valentina notifications@github.com<mailto:notifications@github.com> wrote:

Hi all, I was trying to extract mode and CI from the SEA.B matrix using the code that Andrew provided and using the suggestion of Tor to keep the name of the Community.Group. Everything seems to work except that when I call the matrix with the hdr values, I have multiple columns (2 minimum, named [,1], [,2]) for the same Community.Group$hdr. I don't understand what is the meaning of the different columns, any suggestions? did I miss some math beyond that (probably)?

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/AndrewLJackson/SIBER/issues/9#issuecomment-382001127, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ADK1sE-pWE_2AoGH0lS9RsJ_L5omfP2Gks5tpfQmgaJpZM4F5trw.

[costavale]

Really thank you Andrew. Exactly what I was missing! So since it's a multimodal distribution (I cannot see the figure but I think that I got it) I will have several values related to each distribution. So now I have another question (sorry), when we plot the SEAb or the Layman metrics (after calculating the distribution), in the graph the minimum and the maximum of the boxplot are the min-max through the multimodal distribution?

Thank you Valentina

[AndrewLJackson]

Well that’s an interesting one... I suspect I have not coded it to handle multimodal distributions. You will need to check each boxplot against histograms to be certain. And I or you might have a little coding to do to implement it.

––––––––––––

Dr Andrew Jackson, PhD, FTCD Associate Professor School of Natural Sciences, Department of Zoology Trinity College Dublin, the University of Dublin Dublin 2, Ireland.

+353 1 896 2728 | a.jackson@tcd.iemailto:a.jackson@tcd.ie

Twitter: @yodacomplexhttps://twitter.com/yodacomplex http://www.tcd.ie/Zoology/research/research/theoretical/AndrewJackson.php

Trinity College Dublin, the University of Dublin is ranked 1st in Ireland and in the top 100 world universities by the QS World University Rankings.

On 18 Apr 2018, at 15:37, Valentina notifications@github.com<mailto:notifications@github.com> wrote:

Really thank you Andrew. Exactly what I was missing! So since it's a multimodal distribution (I cannot see the figure but I think that I got it) I will have several values related to each distribution. So now I have another question (sorry), when we plot the SEAb or the Layman metrics (after calculating the distribution), in the graph the minimum and the maximum of the boxplot are the min-max through the multimodal distribution?

Thank you Valentina

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/AndrewLJackson/SIBER/issues/9#issuecomment-382409680, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ADK1sG0L3Acbyh2p4IDFoaKGU9WeU6C5ks5tp0-igaJpZM4F5trw.

[costavale]

I will take a look at all the data and compare histograms and boxplot to understand which one is plotted (I'm not sure I know how to implement your code but I will try). But there was already a warning in your code in the function siberdensityplot.R:

' @section Warning:: This function will not currently recognise and plot

' multimodal distributions, unlike \code{\link[hdrcde]{hdr.boxplot}}. You

' should take care, and plot basic histograms of each variable (column in the

' object you are passing) to \code{siardensityplot} and check that they are

' indeed unimodal as expected.

I should have checked before. Thank you Valentina

[AndrewLJackson]

ah yes, i remember that now! It always annoyed me the way i had coded that siberdensitplot() function. With the benefit of hindsight, i think i should be able to implement a multimodal version easily enough, but i cant promise when i will get a chance to do this. i will add a development note as an issue and will get around to it at some point

thanks for your help and reminding me of this

---

Dr Andrew Jackson, PhD, FTCD Associate Professor School of Natural Sciences, Department of Zoology Trinity College Dublin, the University of Dublin Dublin 2, Ireland.

+353 1 896 2728 | a.jackson@tcd.iemailto:a.jackson@tcd.ie Twitter: @yodacomplexhttps://twitter.com/yodacomplex http://www.tcd.ie/Zoology/research/http://www.tcd.ie/Zoology/research/groups/jackson/groups/jackson/

Trinity College Dublin, the University of Dublin is ranked 1st in Ireland and in the top 100 world universities by the QS World University Rankings.

On 19 Apr 2018, at 08:53, Valentina notifications@github.com<mailto:notifications@github.com> wrote:

I will take a look at all the data and compare histograms and boxplot to understand which one is plotted (I'm not sure I know how to implement your code but I will try). But there was already a warning in your code in the function siberdensityplot.R:

' @Sectionhttps://github.com/Section Warning:: This function will not currently recognise and plot

' multimodal distributions, unlike \code{\link[hdrcde]{hdr.boxplot}}. You

' should take care, and plot basic histograms of each variable (column in the

' object you are passing) to \code{siardensityplot} and check that they are

' indeed unimodal as expected.

I should have checked before. Valentina

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/AndrewLJackson/SIBER/issues/9#issuecomment-382644432, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ADK1sD491VuW5kxs9TpcB5KKZIRJFnUTks5tqEJ8gaJpZM4F5trw.

[AndrewLJackson]

closing this issue and opening another re multimodal density plots