tests/testthat/test_predict.R

In stan-dev/rstanarm: Bayesian Applied Regression Modeling via Stan

# Part of the rstanarm package for estimating model parameters
# Copyright (C) 2015, 2016, 2017 Trustees of Columbia University
# 
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 3
# of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

suppressPackageStartupMessages(library(rstanarm))
library(betareg)
SEED <- 123
set.seed(SEED)
CHAINS <- 2
ITER <- 100
REFRESH <- 0

if (!exists("example_model")) {
  example_model <- run_example_model()
}

plink <- function(fit, nd = NULL, sef = TRUE) 
  predict(fit, newdata = nd, type = "link", se.fit = sef)
presp <- function(fit, nd = NULL, sef = TRUE) 
  predict(fit, newdata = nd, type = "response", se.fit = sef)

context("predict")

test_that("predict recommends posterior_predict for glmer models", {
  expect_error(predict(example_model), 
               "Please use the 'posterior_predict' function")
})

test_that("predict ok for binomial", {
  # example from help(predict.glm)
  ldose <- rep(0:5, 2)
  numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16)
  sex <- factor(rep(c("M", "F"), c(6, 6)))
  SF <- cbind(numdead, numalive = 20-numdead)
  
  glmfit <- glm(SF ~ sex*ldose, family = binomial)
  SW({
    stanfit <- stan_glm(SF ~ sex*ldose, family = binomial, chains = CHAINS, 
                        iter = ITER, seed = SEED, refresh = 0)
    stanfit_opt <- stan_glm(SF ~ sex*ldose, family = binomial, 
                            prior = NULL, prior_intercept = NULL, 
                            seed = SEED, refresh = 0, QR = TRUE,
                            algorithm = "optimizing")
  })
  
  
  pg <- plink(glmfit)
  ps <- plink(stanfit)
  pso <- plink(stanfit_opt)
  expect_equal(pg$fit, ps$fit, tol = 0.1)
  expect_equal(pg$fit, pso$fit, tol = 0.05)
  expect_equal(pg$se.fit, ps$se.fit, tol = 0.2)
  expect_equal(pg$se.fit, pso$se.fit, tol = 0.1)
  expect_equal(presp(glmfit)[1:2], presp(stanfit_opt), tol = 0.05)

  ld <- seq(0, 5, 0.1)
  newd <- data.frame(ldose = ld, sex = factor(rep("M", length(ld)), 
                                              levels = levels(sex)))
  pg <- plink(glmfit, newd)
  ps <- plink(stanfit, newd)
  pso <- plink(stanfit_opt, newd)
  # expect_equal(pg$fit, ps$fit, tol = 0.05)
  expect_equal(pg$fit, pso$fit, tol = 0.05)
  expect_equal(pg$se.fit, ps$se.fit, tol = 0.2)
  expect_equal(pg$se.fit, pso$se.fit, tol = 0.1)
  expect_equal(presp(glmfit, newd)[1:2], presp(stanfit_opt, newd), tol = 0.1)
})

test_that("predict ok for gaussian", {
  glmfit <- glm(mpg ~ wt, data = mtcars)
  SW({
    stanfit <- stan_glm(mpg ~ wt, data = mtcars, chains = CHAINS,
                        iter = 2 * ITER, seed = SEED, refresh = 0)
    stanfit_opt <- stan_glm(mpg ~ wt, data = mtcars,
                            prior = NULL, prior_intercept = NULL,
                            iter = 2 * ITER, seed = SEED, refresh = 0, 
                            algorithm = "optimizing")
  })
  
  pg <- plink(glmfit)
  ps <- plink(stanfit)
  pso <- plink(stanfit_opt)
  expect_equal(pg$fit, ps$fit, tol = 0.05)
  expect_equal(pg$fit, pso$fit, tol = 0.05)
  expect_equal(pg$se.fit, ps$se.fit, tol = 0.3)
  expect_equal(pg$se.fit, pso$se.fit, tol = 0.1)
  expect_equal(presp(glmfit)[1:2], presp(stanfit_opt), tol = 0.1)

  newd <- data.frame(wt = c(1,5))
  pg <- plink(glmfit, newd)
  ps <- plink(stanfit, newd)
  pso <- plink(stanfit_opt, newd)
  expect_equal(pg$fit, ps$fit, tol = 0.05)
  expect_equal(pg$fit, pso$fit, tol = 0.05)
  expect_equal(pg$se.fit, ps$se.fit, tol = 0.3)
  expect_equal(pg$se.fit, pso$se.fit, tol = 0.1)
  expect_equal(presp(glmfit, newd)[1:2], presp(stanfit_opt, newd), tol = 0.1)
})

test_that("predict ok for Poisson", {
  dat <- data.frame(counts = c(18,17,15,20,10,20,25,13,12),
                    outcome = gl(3,1,9), treatment = gl(3,3))

  glmfit <- glm(counts ~ outcome + treatment, data = dat, family = poisson())
  SW({
    stanfit <- stan_glm(counts ~ outcome + treatment, data = dat, family = poisson(), 
                        chains = CHAINS, iter = ITER, seed = SEED, refresh = 0)
    stanfit_opt <- stan_glm(counts ~ outcome + treatment, data = dat, family = poisson(), 
                            iter = ITER, seed = SEED, refresh = 0, algorithm = "optimizing")
  })
  pg <- plink(glmfit)
  ps <- plink(stanfit)
  pso <- plink(stanfit_opt)
  expect_equal(pg$fit, ps$fit, tol = 0.05)
  expect_equal(pg$fit, pso$fit, tol = 0.05)
  expect_equal(pg$se.fit, ps$se.fit, tol = 0.1)
  expect_equal(pg$se.fit, pso$se.fit, tol = 0.1)
  expect_equal(presp(glmfit)[1:2], presp(stanfit_opt), tol = 0.1)

  expect_equal(plink(stanfit, sef = FALSE), plink(glmfit, sef = FALSE), tol = 0.05)
  expect_equal(presp(stanfit, sef = FALSE), presp(glmfit, sef = FALSE), tol = 0.05)

  newd <- dat[1:2, ]
  pg <- plink(glmfit, newd)
  ps <- plink(stanfit, newd)
  pso <- plink(stanfit_opt, newd)
  expect_equal(pg$fit, ps$fit, tol = 0.05)
  expect_equal(pg$fit, pso$fit, tol = 0.05)
  expect_equal(pg$se.fit, ps$se.fit, tol = 0.1)
  expect_equal(pg$se.fit, pso$se.fit, tol = 0.1)
  expect_equal(presp(glmfit, newd)[1:2], presp(stanfit_opt, newd), tol = 0.1)
})