= 1:20
x = 2
a = 1
b = a .+ b * x
y_det = rand.(Normal.(y_det, 2))
y
= lm(hcat(ones(20),x),y)
y_ols
scatter(x,y,labels=false)
plot!(x,y_det,labels="true")
abline!(coef(y_ols)[2], coef(y_ols)[1],labels="best fit") Plots.
Notes on Ch. 5 of “Ecological Models and Data in R”
This chapter is about stochastic simulation and power analysis. We start with examples of stochastic simulation.
using Distributions
using GLM
using StatsFuns
using Optim
using Plots
Linear model example
This is hardly worthy of the term stochastic simulation but here it is.
The model is \(Y \sim N(a+bx,\sigma^2)\).
Hyperbolic negative binomial
The second model is less trivial. The implementation of the negative binomial distribution in Distributions.jl uses a “successes, probability” parameterization of that distribution, so we define a function nbinom_pars
to go to the more “ecological” parameterization in terms of mean and dispersion (arising from a Gamma-Poisson mixture).
The stochastic model is \(Y \sim NegBin(\mu = ab/(b+x), k)\). The dispersion parameter is \(k\) (as \(k \rightarrow \infty\) it becomes a Poisson distribution), \(a\) is an “intercept term” such that when \(x=0\) we get \(\mu = a\).
= 20
a = 1
b = 5
k = rand(Uniform(0,5), 50)
x = sort(x)
x = @. a * b/(b + x)
y_det
# go from (mu,disperson) to (n,p) parameterization
function nbinom_pars(μ,k)
= @. k/(k+μ)
p = @. μ * p/(1-p)
n return n,p
end
# go to Distribution's (n,p) parameterization
= nbinom_pars(y_det, k)
n, p
= rand.(NegativeBinomial.(n,p))
y
# fit the model
# negative loglike of model
function mfun(θ)
= exp.(θ)
a, b, k = @. a/(b+x)
μ = nbinom_pars(μ, k)
n, p return -sum(logpdf.(NegativeBinomial.(n, p), y))
end
= optimize(mfun, log.([20.0, 1.0, 1.0]))
m1
= exp.(Optim.minimizer(m1))
a_opt, b_opt, k_opt = @. a_opt * b_opt/(b_opt + x)
y_opt
scatter(x,y,labels=false)
plot!(x,y_det,labels="true")
plot!(x,y_opt,labels="best fit")
We can do the same thing but with “groups” (fixed effects).
= fill(1,50)
g sample(1:50,25,replace=false)] .= 2
g[
= [20, 10]
a = [1, 2]
b = 5
k
= @. a[g]/(b[g] + x)
y_det
# go to Distribution's (n,p) parameterization
= @. k/(k+y_det)
p = @. y_det * p/(1-p)
n
= rand.(NegativeBinomial.(n,p))
y
# "true" relationship without noise
= @. a[1]/(b[1] + x)
y_grp1 = @. a[2]/(b[2] + x)
y_grp2
scatter(x[g.==1],y[g.==1],labels="group 1",seriescolor=1)
scatter!(x[g.==2],y[g.==2],labels="group 2",seriescolor=2)
plot!(x,y_grp1,label="",seriescolor=1)
plot!(x,y_grp2,label="",seriescolor=2)
Reef fish settlement
This model has “two levels” of stochastic effects. The first is random variation in settlement (density of larvae arriving on a given anemone), and subsequent random variation in density-dependent recruitment (number of settlers surviving on the anemone).
The first level is simulated from a zero-inflated negative binomial model (implemented as rzinbinom
). Then given some number \(S\) of settlers on each anemone, the number of recruits is \(Y \sim Bin(n = S, p = \frac{a}{1+(a/b)S})\). The parameter \(p\) is a decreasing function of density \(S\), hence the density-dependent moniker.
= 603
N = 0.696
a = 9.79
b = 25.32
mu = 0.123
zprob = 0.932
k
function recrprob(S)
/(1+(a/b)*S)
aend
function rzinbinom(N, mu, k, zprob)
= nbinom_pars(mu, k)
n, p = zeros(Int64, N)
samp = rand(Binomial(N,zprob))
nzero = sample(1:N, N-nzero, replace = false)
nonzero .= rand(NegativeBinomial(n,p), N-nzero)
samp[nonzero] return samp
end
= rzinbinom(N,mu,k,zprob)
settlers = @. rand(Binomial(settlers, recrprob(settlers)))
recr
= histogram(settlers,legend=false,xlabel="Settlers",ylabel="Frequency",breaks=40)
p1
= scatter(settlers,recr,xlabel="Settlers",ylabel="Recruits",legend=false)
p2 = plot!(x->a*x/(1+(a/b)*x))
p2
plot(p1,p2)
Pigweed Poisson cluster process
The most complex stochastic simulation considered was a process where locations of parent plants were given as a Poisson process, each with the same number offspr_per_parent
of offspring, which were distributed according to an isotropic distribution around each parent, where displacement followed an exponential distribution.
Each plant has a “competition index” \(C\), counting the number of other plants within a 2m neighborhood. This is used to calculate the end of year mass \(M\), which decays with increased density. The relationship was \(M \sim Gamma(shape = m/(1+C), scale = \alpha)\).
= 30
L = 50
nparents = 10
offspr_per_parent = nparents*offspr_per_parent
noffspr = 2
dispdist = rand(Uniform(0,L), nparents)
parent_x = rand(Uniform(0,L), nparents)
parent_y = rand(Uniform(0,2π), noffspr)
angle = rand(Exponential(dispdist), noffspr)
dist
= vcat([repeat([x],offspr_per_parent) for x in parent_x]...)
offspr_x += @. cos(angle)*dist
offspr_x
= vcat([repeat([y],offspr_per_parent) for y in parent_y]...)
offspr_y += @. sin(angle)*dist
offspr_y
= hcat(offspr_x,offspr_y)
pos = [sqrt(sum((pos[i,:] - pos[j,:]).^2)) for i in axes(pos,1), j in axes(pos,1)]
ndist = [sum(ndist[i,:] .< 2)-1 for i in axes(pos,1)]
nbrcrowd = nbrcrowd*3
ci = 2.3
M = 0.49
alpha = @. M/(1+ci)
mass_det = @. rand(Gamma(alpha, mass_det))
mass = 271.6
b = 0.569
k = b*mass
seed_det
= nbinom_pars(seed_det,k)
n, p = @. rand(NegativeBinomial(n,p))
seed
# θ are the parameters on unconstrained scale
function negbin_ll(θ,y)
= exp(θ[1])
n = logistic(θ[2])
p -loglikelihood(NegativeBinomial(n,p),y)
end
= optimize(x -> negbin_ll(x,nbrcrowd), [log(5),logit(0.5)])
f = exp(Optim.minimizer(f)[1])
fn = logistic(Optim.minimizer(f)[2])
fp
# p1
= scatter(offspr_x,offspr_y,legend=false,alpha=0.5)
p1
# p2
= histogram(nbrcrowd,bins=max(nbrcrowd...),xlabel="Number of neighbors",ylabel="Proportion",normalize=true,legend=false,seriescolor=1)
p2 = scatter!(0:max(nbrcrowd...), map(x -> pdf(NegativeBinomial(fn,fp), x),0:max(nbrcrowd...)),seriescolor=2)
p2 = plot!(0:max(nbrcrowd...), map(x -> pdf(NegativeBinomial(fn,fp), x),0:max(nbrcrowd...)),seriescolor=2)
p2
# p3
= scatter(ci,mass,yscale = :log10,xlabel="Competition index",ylabel="Biomass (g)",legend=false,alpha=0.5)
p3 = plot!(x -> M/(1+x)*alpha)
p3
# p4
function negbin_qt(μ,k,q)
= nbinom_pars(μ,k)
n, p return quantile(NegativeBinomial(n,p), q)
end
= scatter(mass,1 .+ seed,xscale=:log10,yscale=:log10,xlabel="Mass",ylabel="1+Seed set",legend=false,alpha=0.5)
p4 = plot!(x -> b*x+1,xscale=:log10,yscale=:log10,seriescolor=2)
p4 = plot!(x -> negbin_qt(b*x,k,0.975)+1,linestyle=:dash,xscale=:log10,yscale=:log10,seriescolor=2)
p4 = plot!(x -> negbin_qt(b*x,k,0.025)+1,linestyle=:dash,xscale=:log10,yscale=:log10,seriescolor=2)
p4
plot(p1,p2,p3,p4,layout=4)