tests/testthat/test-deviance-residuals.R

In tinyVAST: Multivariate Spatio-Temporal Models using Structural Equations

test_that("deviance residuals for Gamma match glm", {
  skip_on_cran()
  #skip_on_ci()

  set.seed(101)
  x = rnorm(100)
  mu = exp(1 + 0.5*x)
  y = rgamma( n=100, shape=1/0.5^2, scale=mu*0.5^2 )

  # simulated model
  mytiny = tinyVAST( y ~ 1 + x,
            data = data.frame(y=y, x=x),
            family = Gamma(link = "log") )
  resid1 = residuals(mytiny, type="deviance")

  # Null model
  mytiny0 = tinyVAST( y ~ 1,
            data = data.frame(y=y),
            family = Gamma(link = "log") )
  resid0 = residuals(mytiny0, type="deviance")

  #
  myglm = glm( y ~ 1 + x, family=Gamma(link="log"))
  resid2 = residuals( myglm, type="deviance" )
  expect_equal( as.numeric(resid1), as.numeric(resid2),
                tolerance=1e-3 )

  # Compare percent-deviance-explained
  PDEglm = with(summary(myglm), 1 - deviance/null.deviance)
  PDEtiny = (sum(resid0^2)-sum(resid1^2)) / sum(resid0^2)
  expect_equal( PDEglm, PDEtiny,
                tolerance=1e-3 )
})

test_that("nbinom1 and nbinom2 family matches glmmTMB", {
  skip_on_cran()
  #skip_on_ci()
  library(glmmTMB)

  set.seed(101)
  X = rnorm(100)
  eps = 0.5 * rnorm(100)
  p = 2 + 1 * X + eps
  Y = rpois( length(X), lambda = exp(p) )
  data = data.frame(X=X, Y=Y)

  #
  glmmTMB1 = glmmTMB( Y ~ 1 + X, family = glmmTMB::nbinom1(), data = data )
  tinyVAST1 = tinyVAST( Y ~ 1 + X, family = nbinom1(), data = data )
  expect_equal( as.numeric(glmmTMB1$fit$par),
                as.numeric(tinyVAST1$opt$par),
                tolerance=1e-3 )
  expect_equal( as.numeric(resid(glmmTMB1,type="deviance")),
                as.numeric(resid(tinyVAST1,type="deviance")),
                tolerance=1e-3 )

  #
  glmmTMB2 = glmmTMB( Y ~ 1 + X, family = glmmTMB::nbinom2(), data = data )
  tinyVAST2 = tinyVAST( Y ~ 1 + X, family = nbinom2(), data = data )
  expect_equal( as.numeric(glmmTMB2$fit$par),
                as.numeric(tinyVAST2$opt$par),
                tolerance=1e-3 )
  expect_equal( as.numeric(resid(glmmTMB2,type="deviance")),
                as.numeric(resid(tinyVAST2,type="deviance")),
                tolerance=1e-3 )
  if( FALSE ){
    # nbinom2 matches
    plot( x = tinyVAST2$rep$devresid_i,
          y = resid(glmmTMB2, type="deviance") )
    # nbinom1 does match
    plot( x = tinyVAST1$rep$devresid_i,
          y = resid(glmmTMB1, type="deviance") )

    # Manually calculate
    y = Y
    dnbinom = function( x, mu, var, log=FALSE ){
      size = mu^2 / (var - mu)  # number of successful trials
      prob = size / (size + mu)
      stats::dnbinom( x=x, size=size, prob=prob, log=log )
    }

    # Matches
    mu = predict(tinyVAST1)
    phi = exp(tinyVAST1$opt$par[3])
    theta = mu / phi           # log_theta = log_mu - log_phi
    logp1_i = dnbinom(x=y, mu=(y+1e-10), var = (y+1e-10)*(1 + (y+1e-10)/theta), log=TRUE)
    logp2_i = dnbinom(x=y, mu=mu, var = mu*(1 + mu/theta), log=TRUE)
    dev_i = 2 * ( logp1_i - logp2_i )
    devresid_i = sign( y - mu ) * sqrt(dev_i)
    plot( x = resid(glmmTMB1, type="deviance"),
          y =  devresid_i )

    # Matches
    mu = predict(tinyVAST2)
    theta = exp(tinyVAST2$opt$par[3])
    logp1_i = dnbinom(x=y, mu=(y+1e-10), var = (y+1e-10)*(1 + (y+1e-10)/theta), log=TRUE)
    logp2_i = dnbinom(x=y, mu=mu, var = mu*(1 + mu/theta), log=TRUE)
    dev_i = 2 * ( logp1_i - logp2_i )
    devresid_i = sign( y - mu ) * sqrt(dev_i)
    plot( x = resid(glmmTMB2, type="deviance"),
          y =  devresid_i )
    cbind( y, mu, logp1_i, logp2_i, dev_i, devresid_i )

  }

})

test_that("deviance residuals for lognormal match glm", {
  skip_on_cran()
  #skip_on_ci()

  set.seed(101)
  x = rnorm(100)
  logmu = 1 + 0.5*x
  y = rlnorm( n=100, meanlog=logmu, sdlog=0.5 )

  # simulated model
  mytiny = tinyVAST( y ~ 1 + x,
            data = data.frame(y=y, x=x),
            family = lognormal(link = "log") )
  resid1 = residuals(mytiny, type="deviance")

  # Null model
  mytiny0 = tinyVAST( y ~ 1,
            data = data.frame(y=y),
            family = lognormal(link = "log") )
  resid0 = residuals(mytiny0, type="deviance")

  #
  myglm = glm( log(y) ~ 1 + x, family=gaussian(link="identity"))
  resid2 = residuals( myglm, type="deviance" )
  expect_equal( as.numeric(resid1), as.numeric(resid2),
                tolerance=1e-3 )

  # Compare percent-deviance-explained
  PDEglm = with(summary(myglm), 1 - deviance/null.deviance)
  PDEtiny = (sum(resid0^2)-sum(resid1^2)) / sum(resid0^2)
  expect_equal( PDEglm, PDEtiny,
                tolerance=1e-3 )
})

test_that("deviance residuals for tweedie match mgcv", {
  skip_on_cran()
  #skip_on_ci()
  skip_if_not_installed("tweedie")
  skip_if_not_installed("mgcv")

  set.seed(101)
  y = tweedie::rtweedie( n=100, mu=2, phi=1, power=1.5 )

  #
  mytiny = tinyVAST( y ~ 1,
            data = data.frame(y=y),
            family = tweedie(link = "log") )
  resid1 = residuals(mytiny, type="deviance")

  if( FALSE ){
    mu = predict(mytiny)
    phi = exp(mytiny$opt$par[2])
    power = 1 + plogis(mytiny$opt$par[3])
    mytiny$rep$devresid_i^2
    2 * log( dtweedie(y=y, mu=y+1e-10, phi=phi, power=power) / dtweedie(y=y, mu=mu, phi=phi, power=power) )
  }

  #
  library(mgcv)
  mygam = gam( y ~ 1, family=tw(link="log"))
  resid2 = residuals( mygam, type="deviance" )
  expect_equal( as.numeric(resid1), as.numeric(resid2),
                tolerance=1e-3 )
})

test_that("deviance residuals for poisson works", {
  skip_on_cran()
  #skip_on_ci()

  # simulate data
  set.seed(101)
  x = rnorm(100)
  mu = exp(1 + 0.5*x)
  y = rpois( n=100, lambda=mu )

  # delta-gamma in tinyVAST
  mytiny = tinyVAST( y ~ 1 + x,
                     data = data.frame(x=x, y=y),
                     family = poisson("log") )
  resid1 = residuals( mytiny, type="deviance" )
  if( FALSE ){
    resid1^2
    mu = mytiny$rep$mu_i;
    2*y*log(y/mu) - (y-mu)
    sign(y - mu) * pow(2*(y*log((1e-10 + y)/mu) - (y-mu)), 0.5)
    resid1
    pow = function(a,b) a^b
    Type = function(z)z
    sign(y - mu) * pow(2*(y*log((1e-10 + y)/mu) - (y-mu)), 0.5)

    # Compare log-ratio vs.
    plot( x = 2 * log(dpois(y, lambda=y) / dpois(y, lambda=mu)),
          y = mytiny$rep$devresid_i^2 )
  }

  #
  myglm = glm( y ~ 1 + x,
                     data = data.frame(x=x, y=y),
                     family = poisson("log") )
  resid2 = residuals( myglm, type="deviance" )
  expect_equal( as.numeric(resid1), as.numeric(resid2),
                tolerance=1e-3 )
})

test_that("deviance residuals for Bernoulli works", {
  skip_on_cran()
  #skip_on_ci()

  # simulate data
  set.seed(101)
  x = rnorm(100)
  mu = plogis(1 + 0.5*x)
  y = rbinom( n=100, prob=mu, size=1 )

  # delta-gamma in tinyVAST
  mytiny = tinyVAST( y ~ 1 + x,
                     data = data.frame(x=x, y=y),
                     family = binomial("logit") )
  resid1 = residuals( mytiny, type="deviance" )

  #
  myglm = glm( y ~ 1 + x,
                     data = data.frame(x=x, y=y),
                     family = binomial("logit") )
  resid2 = residuals( myglm, type="deviance" )
  expect_equal( as.numeric(resid1), as.numeric(resid2),
                tolerance=1e-3 )
})

test_that("delta-gamma works", {
  skip_on_cran()
  #skip_on_ci()

  # simulate data
  set.seed(101)
  x = rnorm(100)
  mu = exp(1 + 0.5*x)
  y = rgamma( n=100, shape=1/0.5^2, scale=mu*0.5^2 ) *
      rbinom( n=100, size=1, prob=0.5 )

  # delta-gamma in tinyVAST
  mytiny = tinyVAST( y ~ 1 + x,
                     data = data.frame(x=x, y=y),
                     family = delta_gamma(),
                     delta_options = list(
                       formula = ~ 1 + x
                     ) )

  # separate GLMs
  myglm1 = glm( present ~ 1 + x,
                data = data.frame(x=x, present=y>0),
                family = binomial(link="logit") )
  myglm2 = glm( y ~ 1 + x,
                data = subset(data.frame(x=x, y=y),y>0),
                family = Gamma(link="log") )
  
  # relative deviance
  expect_equal( mytiny$rep$deviance, 
                as.numeric(myglm1$deviance + myglm2$deviance),
                tolerance=1e-2 )

  # Compare them
  glmpar = c( coef(myglm1), coef(myglm2) )
  tinypar = mytiny$opt$par[1:4]
  expect_equal( as.numeric(glmpar), as.numeric(tinypar),
                tolerance=1e-3 )
})

test_that("Poisson-link delta-gamma works", {
  skip_on_cran()
  #skip_on_ci()
  skip_if_not_installed("sdmTMB")

  # simulate data
  set.seed(101)
  x = rnorm(100)
  mu = exp(1 + 0.5*x)
  y = rgamma( n=100, shape=1/0.5^2, scale=mu*0.5^2 ) *
      rbinom( n=100, size=1, prob=0.5 )

  # delta-gamma in tinyVAST
  mytiny = tinyVAST( y ~ 1 + x,
                     data = data.frame(x=x, y=y),
                     family = delta_gamma(type="poisson-link"),
                     delta_options = list(
                       formula = ~ 1 + x
                     ),
                     control = tinyVASTcontrol(verbose=FALSE) )

  # delta-gamma in sdmTMB
  mysdmTMB = sdmTMB::sdmTMB( list( y ~ 1 + x, y ~ 1 + x ),
                     data = data.frame(x=x, y=y),
                     family = delta_gamma(type="poisson-link"),
                     spatial = "off" )

  # separate GLMs
  myglm1 = glm( present ~ 1 + x,
                data = data.frame(x=x, present=y>0),
                family = binomial(link="cloglog") )

  # Compare them
  glmpar = c( coef(myglm1) )
  tinypar = mytiny$opt$par[1:4]
  sdmTMBpar = mysdmTMB$model$par[1:4]
  expect_equal( as.numeric(tinypar), as.numeric(sdmTMBpar),
                tolerance=1e-3 )
})