library(spatstat)
window <- as.owin(bei)

# Pilot model
D_pilot <- runifpoint(1e4, window)
Q_pilot <- quadscheme.logi(bei, D_pilot)
fit_pilot <- ppm(Q_pilot  ~ polynom(elev, grad, 2),
                 data = bei.extra, method = "logi")

# Pilot intensity estimate
lambda_pilot <- predict(fit_pilot, locations = window)
lambda_pilot <- 1e4 * lambda_pilot / integral(lambda_pilot)

# Second-stage model
D_local <- rpoint(1e4, lambda_pilot)
Q_local <- quadscheme.logi(bei, D_local)
rho <- Q_local$param$rho
fit_local <- ppm(Q_local ~ polynom(elev, grad, 2) +
                   offset(-log(lambda_pilot) + log(rho)),
                 data = bei.extra, method = "logi")

# Final intensity estimate
plot(predict(fit_local, locations = window) * lambda_pilot)