For models fit using MCMC only, the log_lik method returns the \(S\) by \(N\) pointwise log-likelihood matrix, where \(S\) is the size of the posterior sample and \(N\) is the number of data points, or in the case of the stanmvreg method (when called on stan_jm model objects) an \(S\) by \(Npat\) matrix where \(Npat\) is the number of individuals.

# S3 method for stanreg
log_lik(object, newdata = NULL, offset = NULL, ...)

# S3 method for stanmvreg
log_lik(object, m = 1, newdata = NULL, ...)

# S3 method for stanjm
log_lik(object, newdataLong = NULL, newdataEvent = NULL,
  ...)

Arguments

object

A fitted model object returned by one of the rstanarm modeling functions. See stanreg-objects.

newdata

An optional data frame of new data (e.g. holdout data) to use when evaluating the log-likelihood. See the description of newdata for posterior_predict.

offset

A vector of offsets. Only required if newdata is specified and an offset was specified when fitting the model.

...

Currently ignored.

m

Integer specifying the number or name of the submodel

newdataLong, newdataEvent

Optional data frames containing new data (e.g. holdout data) to use when evaluating the log-likelihood for a model estimated using stan_jm. If the fitted model was a multivariate joint model (i.e. more than one longitudinal outcome), then newdataLong is allowed to be a list of data frames. If supplying new data, then newdataEvent should also include variables corresponding to the event time and event indicator as these are required for evaluating the log likelihood for the event submodel. For more details, see the description of newdataLong and newdataEvent for posterior_survfit.

Value

For the stanreg and stanmvreg methods an \(S\) by \(N\) matrix, where \(S\) is the size of the posterior sample and \(N\) is the number of data points. For the stanjm method an \(S\) by \(Npat\) matrix where \(Npat\) is the number of individuals.

Examples

roaches$roach100 <- roaches$roach1 / 100 fit <- stan_glm( y ~ roach100 + treatment + senior, offset = log(exposure2), data = roaches, family = poisson(link = "log"), prior = normal(0, 2.5), prior_intercept = normal(0, 10), iter = 500 # to speed up example )
#> #> SAMPLING FOR MODEL 'count' NOW (CHAIN 1). #> #> Gradient evaluation took 0.000116 seconds #> 1000 transitions using 10 leapfrog steps per transition would take 1.16 seconds. #> Adjust your expectations accordingly! #> #> #> Iteration: 1 / 500 [ 0%] (Warmup) #> Iteration: 50 / 500 [ 10%] (Warmup) #> Iteration: 100 / 500 [ 20%] (Warmup) #> Iteration: 150 / 500 [ 30%] (Warmup) #> Iteration: 200 / 500 [ 40%] (Warmup) #> Iteration: 250 / 500 [ 50%] (Warmup) #> Iteration: 251 / 500 [ 50%] (Sampling) #> Iteration: 300 / 500 [ 60%] (Sampling) #> Iteration: 350 / 500 [ 70%] (Sampling) #> Iteration: 400 / 500 [ 80%] (Sampling) #> Iteration: 450 / 500 [ 90%] (Sampling) #> Iteration: 500 / 500 [100%] (Sampling) #> #> Elapsed Time: 0.211255 seconds (Warm-up) #> 0.106999 seconds (Sampling) #> 0.318254 seconds (Total) #> #> #> SAMPLING FOR MODEL 'count' NOW (CHAIN 2). #> #> Gradient evaluation took 4.5e-05 seconds #> 1000 transitions using 10 leapfrog steps per transition would take 0.45 seconds. #> Adjust your expectations accordingly! #> #> #> Iteration: 1 / 500 [ 0%] (Warmup) #> Iteration: 50 / 500 [ 10%] (Warmup) #> Iteration: 100 / 500 [ 20%] (Warmup) #> Iteration: 150 / 500 [ 30%] (Warmup) #> Iteration: 200 / 500 [ 40%] (Warmup) #> Iteration: 250 / 500 [ 50%] (Warmup) #> Iteration: 251 / 500 [ 50%] (Sampling) #> Iteration: 300 / 500 [ 60%] (Sampling) #> Iteration: 350 / 500 [ 70%] (Sampling) #> Iteration: 400 / 500 [ 80%] (Sampling) #> Iteration: 450 / 500 [ 90%] (Sampling) #> Iteration: 500 / 500 [100%] (Sampling) #> #> Elapsed Time: 0.140815 seconds (Warm-up) #> 0.135274 seconds (Sampling) #> 0.276089 seconds (Total) #> #> #> SAMPLING FOR MODEL 'count' NOW (CHAIN 3). #> #> Gradient evaluation took 4.4e-05 seconds #> 1000 transitions using 10 leapfrog steps per transition would take 0.44 seconds. #> Adjust your expectations accordingly! #> #> #> Iteration: 1 / 500 [ 0%] (Warmup) #> Iteration: 50 / 500 [ 10%] (Warmup) #> Iteration: 100 / 500 [ 20%] (Warmup) #> Iteration: 150 / 500 [ 30%] (Warmup) #> Iteration: 200 / 500 [ 40%] (Warmup) #> Iteration: 250 / 500 [ 50%] (Warmup) #> Iteration: 251 / 500 [ 50%] (Sampling) #> Iteration: 300 / 500 [ 60%] (Sampling) #> Iteration: 350 / 500 [ 70%] (Sampling) #> Iteration: 400 / 500 [ 80%] (Sampling) #> Iteration: 450 / 500 [ 90%] (Sampling) #> Iteration: 500 / 500 [100%] (Sampling) #> #> Elapsed Time: 0.148195 seconds (Warm-up) #> 0.117762 seconds (Sampling) #> 0.265957 seconds (Total) #> #> #> SAMPLING FOR MODEL 'count' NOW (CHAIN 4). #> #> Gradient evaluation took 5e-05 seconds #> 1000 transitions using 10 leapfrog steps per transition would take 0.5 seconds. #> Adjust your expectations accordingly! #> #> #> Iteration: 1 / 500 [ 0%] (Warmup) #> Iteration: 50 / 500 [ 10%] (Warmup) #> Iteration: 100 / 500 [ 20%] (Warmup) #> Iteration: 150 / 500 [ 30%] (Warmup) #> Iteration: 200 / 500 [ 40%] (Warmup) #> Iteration: 250 / 500 [ 50%] (Warmup) #> Iteration: 251 / 500 [ 50%] (Sampling) #> Iteration: 300 / 500 [ 60%] (Sampling) #> Iteration: 350 / 500 [ 70%] (Sampling) #> Iteration: 400 / 500 [ 80%] (Sampling) #> Iteration: 450 / 500 [ 90%] (Sampling) #> Iteration: 500 / 500 [100%] (Sampling) #> #> Elapsed Time: 0.18563 seconds (Warm-up) #> 0.11219 seconds (Sampling) #> 0.29782 seconds (Total) #>
ll <- log_lik(fit) dim(ll)
#> [1] 1000 262
all.equal(ncol(ll), nobs(fit))
#> [1] TRUE
# using newdata argument nd <- roaches[1:2, ] nd$treatment[1:2] <- c(0, 1) ll2 <- log_lik(fit, newdata = nd, offset = c(0, 0)) head(ll2)
#> 1 2 #> [1,] -7.208039 -3.379939 #> [2,] -7.538332 -3.645910 #> [3,] -6.475576 -3.415374 #> [4,] -6.957932 -3.350174 #> [5,] -5.320412 -3.461726 #> [6,] -6.218462 -3.344113
dim(ll2)
#> [1] 1000 2
all.equal(ncol(ll2), nrow(nd))
#> [1] TRUE