0

This might be quite basic but I'm struggling.

I'm trying to take a published ecological model and implement it in the Turing.jl package using Julia. Say I have a number of species, $S$, and I want to find out the probability of any of them eating another based off of a series of characteristics (traits) of each (i.e. body size, habitat preference etc). The probability of a predator $i$ eating a prey item $j$ is given as:

$$P(i,j,\theta) = \alpha \prod^D_{d=1} \exp\left(- \left|\frac{n_{d,j}-c_{d,i}}{r_{d,i}/2}\right|^e\right).$$

where $\theta = \{ n_{1,1}...n_{D,S},c_{1,1}...c_{D,S},r_{1,1}...r_{D,S},e\}$, the parameter $n_{d,j}$ is the trait value of $j$, the prey, in dimension $d$, $c_{d,i}$ is the optimal feeding position of species $i$ in dimension $d$; $r_{d,i}$ is the feeding range of speices $i$ in dimension $d$, $e$ just varies the shape of the probability function, and $\alpha$ is the probability of $i$ eating $j$ if the trait value of the prey $n_{d,j}$ is equal to the optimal value of the predator $c_{d,i}$.

To reduce the number of free parameters, we can set the range $r_{d,i}$ to a large value (i.e. 10) and $c_{d,i}$ to the centre of the range (i.e. 5) and then simply fit $n_{j,i}$.

The data, $\textbf{X}$ are $S \times S$ matrices containing an observation $X_{ij}$ for each instance $i,j$. ($X_{ij}$ = 1 means $i$ eats $j$, $X_{ij}$ = 0 means $i$ does not eat $j$. The matrices are not especially sparse, (roughly 2/3 1s and 1/3 0s).

Obviously there can only be one outcome (species $j$ is eaten by species $i$ or not eaten AKA 1 or 0), so my goal is to formulate the problem as a Bayesian logistic regression, which takes the form $$ Y \sim Bernoulli(p)$$ but I am unclear how to formulate the probability equation above into some form that I can pass as $p$.

Extra (Possibly Useful?) Info:

In the published version of the model the authors state for the dataset $\textbf{X}$ they estimate the the maximum likelihood parameter set for each link $X_ij$ according to the model, and use simulated annealing.

They state their log-likelihood is:

$$\ell(\textbf{X}|\theta) = \sum_i \sum_j \text{ln}\left\{\begin{array}{lr} P(i,j|\theta), & \text{if } X_{ij}= 1\\ 1-P(i,j|\theta), & \text{if } X_{ij}= 0 \end{array}\right\}$$

I am hoping to frame this in a Bayesian logistic regression as mentioned above.

colebrookson
  • 113
  • 4
  • 2
    Which of the terms $\eta_j,$ $c_i,$ $r_i,$ and $\alpha$ are known and which are ones you need to estimate? This makes such a *huge* difference in the solution that we need clarification before we can proceed. – whuber Jan 12 '22 at 23:03
  • 1
    This is an interesting question, but it looks like it needs more work in its formulation because its parameters are not identifiable. Specifically, if you think you have estimates of $(r_i,c_i,\eta_j),$ then for any nonzero number $\lambda$ the values $(\lambda r_i,\lambda c_i,\lambda \eta_j)$ work equally well because the exponents of $\alpha$ are unchanged. I suggest that you consider changing your question to one in which you describe your substantive practical problem and ask for suggestions about modeling it effectively. – whuber Jan 13 '22 at 15:57
  • Moreover, it sounds like your data will be predominantly zeros, which will make it impossible to pin down any of the associated parameters, creating questions concerning what your data look like and how you have selected them. That information belongs in your question, too, for otherwise you risk getting good mathematical, algorithmic, or programming answers that, when applied, give you terribly wrong information. – whuber Jan 13 '22 at 15:59
  • 1
    @whuber thank you so much for these comments, very helpful. I'm going to reformat the question to include the information you mention here, and also, going back to the paper that I'm trying to reproduce the results from, I found a trick to reduce the number of free parameters. I'll make edits to the question and re-post. Thanks again!! – colebrookson Jan 13 '22 at 18:45
  • 1
    The problem is becoming more complicated! It sounds like one key aspect is whether you can assume--even approximately--that all the responses are independent. For instance, to what extent might an observation of $i$ eating $j$ change the chances of other observations of what $i$ eats? After all, if $i$ needs a balanced diet and that requires variation in $j,$ over short periods one would suppose there is some negative association between the events. If you *can* assume independence, you're basically done, because you just plug your formula into the expression for the log likelihood. – whuber Jan 13 '22 at 22:21
  • Without getting into the nitty gritty, I definitely *can* assume independence (I verified this a while ago with supervisors etc), so it sounds like that's okay? Maybe I'm still being slow on the uptake, but so does that mean I can essentially write my sampling of $Y \sim Bernoulli(p)$ where p is just the equation for $P(i,j,\theta)$? – colebrookson Jan 13 '22 at 22:44
  • 1
    Yes, that's right. "$\theta$" stands for all the parameters you wish to estimate. – whuber Jan 13 '22 at 23:07
  • great, thank you!!! – colebrookson Jan 13 '22 at 23:20

0 Answers0