Consider adding the LSW patch size distribution

[knstr]

Hi,
a suggestion to add an option to fit the Lifshitz-Slyozov-Wagner distribution to patch size data.

This must be done numerically, as described in:
Tiryakioğlu, et al., On evaluating fit of the Lifshitz–Slyozov–Wagner (LSW) distribution to particle size data. Mater. Sci. Eng. A 527, 1636–1639 (2010)

Patch size distributions are expected to converge to the LSW distribution as cover approaches 0. That is, as conditions become harsher. This is expected in systems with a conservative resource (or consumer), whose mobility is controlled by an ecosystem engineer. See:
Siteur, et al., 2023. Phase-separation physics underlies new theory for the resilience of patchy ecosystems. [PNAS 120 (2): 1-12].

Happy to help out, if you decide to include this option.

Groet!
Koen

[alexgenin]

Hi Koen, thanks for your interest ! This would be a great addition to spatialwarnings indeed ! In fact your new paper is already near the top of my reading list ;)

What would really help is to have a reference implementation for this type of things, so that the spatialwarnings results can be checked against something for correctness. If you have any pointer to an existing piece of code that implements the distribution fitting it would be great ! Another useful contribution can be datasets (samples), with the distribution parameter values expected from a fit.

Thanks,
Alex

[knstr]

Great,
I'll have to look up my code (which is in Matlab), but it boils down to this:
The LSW distribution is given by:
$$g(r;u)=\frac{4r^2}{9\mu^3}\left(\frac{3\mu}{3\mu+r}\right)^{7/3}\left(\frac{3\mu}{3\mu-2r}\right)^{11/3}\exp\left(\frac{2r}{2r-3\mu}\right)$$
for $0\lt r\lt 3\mu /2$ and $g(r;\mu)=0$ for $r\geq 3\mu /2$.
Here $r$ is the patch radius and $\mu$ is the mean patch radius. The patch radius is given by $\sqrt{a/2\pi}$, where $a$ is the patch area.
So the distribution has only one parameter $\mu$. A rough first estimate of $\mu$ is the average patch radius of the sample. However, this estimate does not yield the best fit. The (log) likelihood is maximized by finding the value for $\mu$ at which $dL/d\mu=0$. This means that one has to numerically solve:

$$\frac{dL}{d\mu}=-3n+ \frac{7}{3}\sum_{i=1}^n\frac{r_i}{3\mu+r_i} -\frac{11}{3}\sum_{i=1}^n\frac{2r_i}{3\mu-2r_i} +3\mu \sum_{i=1}^n\frac{2r_i}{(2r_i-3\mu)^2}=0$$

See the paper of Tiryakioğlu. I guess in R this could be done with "optim" or another function?

I'll try to find my code and sample data.

In the end the relevant indicators are

1. a decline in cover,
2. a decline in skewness (towards the theoretical limit of -0.92)
3. a shift from a log normal to the LSW distribution.

See the figure below from the Supplementary Material.

[knstr]

Hi!
I made some code to fit the LSW distribution in R.
I hope this helps.

Groet!
Koen

fit_LSW.zip

[alexgenin]

Hi Koen, this is great ! I just created a first draft of what things could look like. Right now we just return the difference in AIC, cover, and skewness. This is trivial to do in spatialwarnings so I could do it now, but giving the user the ability to plot distributions, etc requires a bit more of code-writing. I don't have a dataset to test this on, so it'll have to wait a little bit, but I wanted you to have a look at what it could look like.

devtools::install_github("spatial-ews/spatialwarnings") 
library(spatialwarnings)
indic <- spatialwarnings:::lsw_sews(forestgap) # makes no real sense for this dataset, but the code runs
plot(indic)

see also https://github.com/spatial-ews/spatialwarnings/blob/master/R/lsw_sews.R

[JohanvandeKoppel]

Great initiative!

[knstr]

Hi Alex,
great work. If neccessary I can do a number of modelruns to get a dataset that is more appropriate than the forestgap pattern?

Groet,
Koen

[alexgenin]

Hi Koen - it would be brilliant. We could include it in the package (if it's not too big) as a demo dataset. Don't forget to save the model parameters as well!

Cheers,
Alex

[knstr]

Hi Alex,
please find attached a list of 1000x1000 logical matrices, representing an area of 100x100 dimensionless space.

u_cover_DDAmodel.zip

I made the matrices with the DDA model in my paper, with parameters $\alpha=4$, $\beta=1$ and $\delta=0.01$. What you see is the cover of state variable $u$, which was obtained by applying a threshold at $u=7.11$. This threshold seems arbitrary, but it's the unstable equilibrium of the model for this parameter setting ($\tilde u_-$ in the paper). I changed $\langle \tau \rangle$ from 3.25 to 10. The names of the matrices in the list give the value of $\langle \tau \rangle$.

Unfortunately I find some unexpected behavior in the "aic_diff" plot, as you can see below:

It may be caused by a bug I found in the code I sent you. It appears that entering mu=min(r) into the function dLdmu(r,mu) returns NA. Therefore, the lower limit of the interval in the uniroot function should not be equal to min(r).

One solution is adding a small amount to the lower limit of the interval:

fitLSW <- function(r){
  mu_fit <- uniroot(function(mu){dLdmu(r,mu)},
                    interval=c(min(r)+0.001,
                               max(r)-0.001))$root
  AIC_fit <- 2-2*sum(log(LSW(r,mu_fit)))
  return(list(mu=mu_fit,AIC=AIC_fit))
}

I haven't tested if this solves the issue with the "aic_diff" plot.

[alexgenin]

Thanks, it works now !

I think the numerical issues stem from the fact that the log-likelihood is infinite if any observed patch has a radius above 3mu/2. Restricting the optim/uniroot search to values above a minimum 2*max(r)/3 + { a little number found by bisection } guarantees we will always work on the range of values for which the LL is finite. Note that in the current version I just used a direct negative LL minimization instead of finding the zero of the derivative, but both produce the same results (I'll add that to the unit tests later).

The skewness produces values below -0.92, but I'm not sure what to think about it, if that's a mistake or a sampling effect.

So it seems everything works well ! The last remaining thing to do is to add a little bit of documentation. lsw_sews needs a description saying in which context these indicators would be useful, a little bit of the background + references behind them -- anything the user needs to know/pay attention to before actually using it. You're probably the best person to do it, if you want to. An example can be found in the Details section of the function that computes generic indicators.

This will get us a releasable working version of lsw_ews to compute skewness/proba of LSW distribution and cover (i.e. produce the graph above). It would be nice next to extend the function to draw the fitted distributions, but that can come later.

The simulation dataset is distributed as dda/dda.pars, e.g.

install_github("spatial-ews/spatialwarnings")
indics <- lsw_sews(dda)
plot(indics, along = dda.pars[ ,"tau"]) + 
  labs(x = expression(tau))

We'll probably have to reduce its size for CRAN, but that's something for later.

[knstr]

Great work!
I'll write some documentation when I find the time.
Great that the model output it's now a demoset. Maybe we can reduce the number of matrices instead of the matrix size, for CRAN?

About the skewness below 0.92. To be honest, this is the first time to see a skewness below this value for this model. Three possible explanations could be:

1. sample size
2. the LSW distribution was derived for the Cahn Hilliard model, which has a potential function (P(u)=u^3-u) with the unstable equilibrium exactly in between the two stable equilibria.
3. the LSW distribution was derived for three spatial dimensions.

I wouldn't know what the cause could be, but I guess it's not a problem for ecological purposes.

Groet!
Koen

[alexgenin]

Hi Koen,

I have spent some time improving the code related to all this. Users will have access to lsw_sews to obtain the plots above, as well as some utility functions to work with and fit LSW distributions (LSW_fit, d/pLSW). I wrote some documentation based off my understanding of your paper, but it would probably benefit from a quick read by you, could you have a read to the documentation in ./R/lsw_sews.R? Thanks !

If you have also relevant papers/books regarding the LSW distribution, that would be great.

This will make it into the next release of spatialwarnings (mid-january I hope). So your input would be very much appreciated.

Cheers,

Alex