Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Self-attribution from model #24

Open
domenico-simone opened this issue Sep 22, 2022 · 3 comments
Open

Self-attribution from model #24

domenico-simone opened this issue Sep 22, 2022 · 3 comments

Comments

@domenico-simone
Copy link

Hi,

is there a way to get P of assignments of isolates to their own clusters? I.e., self-attribution results to assess the model reliability?

Thank you,

Domenico

@domenico-simone domenico-simone changed the title self-attribution from model Self-attribution from model Sep 22, 2022
@jmarshallnz
Copy link
Owner

jmarshallnz commented Sep 26, 2022

If what you're after is P(Source | ST) from the 'genotype distribution' part of the model (asymmetric island or dirichlet) then you can do this by using the likelihood returned by st_fit() along with a prior. st_fit() for the AI model returns a leave one out estimate of log(P(ST | Source)). So you can convert that to P(Source | ST) via Bayes Thm.,

summary() gives you the posterior mean of each source given ST assuming apriori each source equally likely. You could adapt the code there to give you the whole posterior. One thing to watch is in the exponentiation of logP - it pays to do it by first factoring out the largest value so you don't lose precision.

Code:

library(islandR)
library(tidyverse)
st = st_fit(formula = Source ~ ST,
            non_primary = "Human",
            data = manawatu,
            method="island",
            sequences = ~ ASP + GLN + GLT + GLY + PGM + TKT + UNC)

prior <- tibble(source = setdiff(unique(manawatu$Source), "Human"),
                prior = 1/4)

exp_sum_log <- function(x) { y = x - max(x); exp(y)/sum(exp(y)) }

st |> as.data.frame() |>
  left_join(prior) |>
  mutate(log_pp = log(prior) + log_p) |>
  group_by(type, iteration) |>
  mutate(posterior = exp_sum_log(log_pp)) |>
  select(type, iteration, posterior)

@jmarshallnz
Copy link
Owner

Note that this is not including the 'attribution' part of the model, which infers the mixing (here the fixed prior) of the sources.

Other SA models may be utilising that 'inferred mixing' when doing self attribution. This might be particularly important for imbalanced sources, though you may be able to incorporate that via the prior above.

@jmarshallnz
Copy link
Owner

See 22e3215

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants