simulation_3p_IRT.R

library(brms)
library(tidyverse)

# we're going to have 5 animals
nspp <- 5
# in 100 sites
nsites <- 100

#abilities are between -1 and 1, from a normal distribution:
thetas <- rnorm(nspp, 0, 1.2)

# site difficulty decreases with x
b_x <- - 0.7
# the average site is an OK place to live, not great but OK
b_0 <- -1

# correlation between slopes and intercepts for sites
site_cor <- matrix(c(1, -0.3, -0.3, 1), nrow = 2)
# site variance in slopes and inters
sigmas <- c(0.8, 0.2)

site_cov <- diag(sigmas) %*% site_cor %*% diag(sigmas)

site_varying_effects <- MASS::mvrnorm(nsites, mu = c(0, 0), Sigma = site_cov)

# site x-values
xs <- runif(nsites, min = -3, max = 3)
hist(xs)

# site "difficulty"
beta_site <- b_0 + xs * b_x + site_varying_effects[,1]

# site discrimination
a0 <- 0.1

alpha_site <- exp(a0 + site_varying_effects[,2])

projected_response <- expand.grid(site = 1:nsites, species = 1:nspp) %>% 
  tbl_df %>% 
  mutate(
    site_x = xs[site],
    site_beta = beta_site[site],
    site_alpha = alpha_site[site],
    spp_theta       = thetas[species],
    response = site_alpha * (spp_theta - site_beta),
    probability = rethinking::inv_logit(response),
    pres_abs = rbinom(length(probability), 1, prob = probability))


projected_response %>% 
  ggplot(aes(x = site_x, y = pres_abs, colour = species, group = species)) + geom_point() +
  stat_smooth(method = "glm", method.args = list(family = "binomial"))

inv_logit <- function(x) 1 / (1 + exp(-x))

irt_3p <- bf(pres_abs ~ inv_logit(alpha * (theta - beta)), 
   alpha ~ 0 + (1|site),
   beta ~ 0 + intercept + site_x + (1|site),
   theta ~ 0 + (1 | species),
   nl = TRUE)

get_prior(irt_3p, data = projected_response)

my_priors <- c(set_prior("cauchy(0, 3)", class = "sd", nlpar = "alpha"),
               set_prior("cauchy(0, 3)", class = "sd", nlpar = "beta"),
               set_prior("cauchy(0, 3)", class = "sd", nlpar = "theta"),
               set_prior("normal(0, 2)", class = "b", coef = "intercept", nlpar = "beta"),
               set_prior("normal(0, 2)", class = "b", coef = "site_x", nlpar = "beta") # lb = 0)
)

fit_ir3 <- brm(irt_3p,
               data = projected_response, family = bernoulli("identity"), 
               prior = my_priors)

summary(fit_ir3)