Playing with probability-of-presence from Maxent-type models

pop_maxent.R

# converting maxent/Poisson model output to true probability of presence

# clear workspace
rm(list = ls())

# set rng seed
set.seed(1)

# number of datapoints
n <- 100000

# 'true' model
f <- function(x) pnorm(-3 * sin(x * 0.6))

# covariate
x <- runif(n, -6, 6)

# 'true' probability
p <- f(x)

# observed presence/absence
y <- rbinom(n, 1, p)

# x where at presence locations: \pi = p(x|y=1)
x_y1 <- x[y == 1]

adj <- 0.6

# Density estimate of p(x|y=1)
f1_star <- density(x_y1,
                   adjust = adj,
                   from = -4,
                   to = 4)

# Density estimate of p(x)
f_star <- density(sample(x, length(x_y1), replace = FALSE),
                  adjust = adj,
                  from = -4,
                  to = 4)

# ~~~~~~~
# Attempts to reproduce p(y=1|x)

# Quantities: 
# p(y=1) == prevalence in x == mean(p) ~= mean(y)
py1 <- mean(f(f1_star$x))

# Answer using Bayes theorem:
# p(y=1|x) = (p(x|y=1) * p(y=1)) / (p(x))

plot(f1_star$y / f_star$y * py1 ~ f1_star$x,
     col = 'blue',
     type = 'l',
     ylim = c(-0.3, 1.3))

# Answer using Dave's approximation
# Z = mean(p(x|y=1)) / p(y=1)
# p(y=1|x) = p(x|y=1) / Z

Z <- mean(f1_star$y) / py1
lines(f1_star$y / Z ~ f1_star$x,
      col = 'red')

# true answer
lines(f(f1_star$x) ~ f1_star$x,
      lty = 2)

# add in lines for 0 and 1
abline(h = c(0, 1),
       col = 'grey')