% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/walker.R
\name{walker_glm}
\alias{walker_glm}
\title{Bayesian generalized linear model with time-varying coefficients}
\usage{
walker_glm(
  formula,
  data,
  beta,
  init,
  chains,
  return_x_reg = FALSE,
  distribution,
  initial_mode = "kfas",
  u,
  mc_sim = 50,
  return_data = TRUE,
  ...
)
}
\arguments{
\item{formula}{An object of class \code{{formula}} with additional terms
\code{rw1} and/or \code{rw2} e.g. \code{y ~ x1 + rw1(~ -1 + x2)}. See details.}

\item{data}{An optional data.frame or object coercible to such, as in \code{{lm}}.}

\item{beta}{A length vector of length two which defines the
prior mean and standard deviation of the Gaussian prior for time-invariant coefficients}

\item{init}{Initial value specification, see \code{\link[rstan:stanmodel-method-sampling]{rstan::sampling()}}.
Note that compared to default in \code{rstan}, here the default is a to sample from the priors.}

\item{chains}{Number of Markov chains. Default is 4.}

\item{return_x_reg}{If \code{TRUE}, does not perform sampling, but instead returns the matrix of
predictors after processing the \code{formula}.}

\item{distribution}{Either \code{"poisson"} or \code{"binomial"}.}

\item{initial_mode}{The initial guess of the fitted values on log-scale.
Defines the Gaussian approximation used in the MCMC.
Either \code{"obs"} (corresponds to log(y+0.1) in Poisson case),
\code{"glm"} (mode is obtained from time-invariant GLM), \code{"mle"}
(default; mode is obtained from maximum likelihood estimate of the model),
or numeric vector (custom guess).}

\item{u}{For Poisson model, a vector of exposures i.e. \eqn{E(y) = u*exp(x*beta)}.
For binomial, a vector containing the number of trials. Defaults 1.}

\item{mc_sim}{Number of samples used in importance sampling. Default is 50.}

\item{return_data}{if \code{TRUE}, returns data input to \code{\link[rstan:stanmodel-method-sampling]{rstan::sampling()}}.
This is needed for \code{lfo}.}

\item{...}{Further arguments to \code{\link[rstan:stanmodel-method-sampling]{rstan::sampling()}}.}
}
\value{
A list containing the \code{stanfit} object, observations \code{y},
covariates \code{xreg_fixed}, and \code{xreg_rw}.
}
\description{
Function \code{walker_glm} is a generalization of \code{walker} for non-Gaussian
models. Compared to \code{walker}, the returned samples are based on Gaussian approximation,
which can then be used for exact-approximate analysis by weighting the sample properly. These weights
are also returned as a part of the \code{stanfit} (they are generated in the
generated quantities block of Stan model). Note that plotting functions \code{pp_check},
\code{plot_coefs}, and \code{plot_predict} resample the posterior based on weights
before plotting, leading to "exact" analysis.
}
\details{
The underlying idea of \code{walker_glm} is based on Vihola, Helske, Franks (2020).

\code{walker_glm} uses the global approximation (i.e. start of the MCMC) instead of more accurate
but slower local approximation (where model is approximated at each iteration).
However for these restricted models global approximation should be sufficient,
assuming the the initial estimate of the conditional mode of p(xbeta | y, sigma) not too
far away from the true posterior. Therefore by default \code{walker_glm} first finds the
maximum likelihood estimates of the standard deviation parameters
(using \code{\link[KFAS:KFAS]{KFAS::KFAS()}}) package, and
constructs the approximation at that point, before running the Bayesian
analysis.
}
\examples{

set.seed(1)
n <- 25
x <- rnorm(n, 1, 1)
beta <- cumsum(c(1, rnorm(n - 1, sd = 0.1)))

level <- -1
u <- sample(1:10, size = n, replace = TRUE)
y <- rpois(n, u * exp(level + beta * x))
ts.plot(y)

# note very small number of iterations for the CRAN checks!
out <- walker_glm(y ~ -1 + rw1(~ x, beta = c(0, 10), 
  sigma = c(2, 10)), distribution = "poisson", 
  iter = 200, chains = 1, refresh = 0)
print(out$stanfit, pars = "sigma_rw1") ## approximate results
if (require("diagis")) {
  weighted_mean(extract(out$stanfit, pars = "sigma_rw1")$sigma_rw1, 
    extract(out$stanfit, pars = "weights")$weights)
}
plot_coefs(out)
pp_check(out)
             
\dontrun{
data("discoveries", package = "datasets")
out <- walker_glm(discoveries ~ -1 + 
  rw2(~ 1, beta = c(0, 10), sigma = c(2, 10), nu = c(0, 2)), 
  distribution = "poisson", iter = 2000, chains = 1, refresh = 0)

plot_fit(out)

# Dummy covariate example

fit <- walker_glm(VanKilled ~ -1 + 
  rw1(~ law, beta = c(0, 1), sigma = c(2, 10)), dist = "poisson", 
   data = as.data.frame(Seatbelts), chains = 1, refresh = 0)

# compute effect * law
d <- coef(fit, transform = function(x) {
  x[, 2, 1:170] <- 0
  x
})

require("ggplot2")
d \%>\% ggplot(aes(time, mean)) +
  geom_ribbon(aes(ymin = `2.5\%`, ymax = `97.5\%`), fill = "grey90") +
  geom_line() + facet_wrap(~ beta, scales = "free") + theme_bw()
}

}
\references{
Vihola, M, Helske, J, Franks, J. (2020). Importance sampling
type estimators based on approximate marginal Markov chain Monte Carlo.
Scandinavian Journal of Statistics. 47: 1339–1376. \doi{doi:10.1111/sjos.12492}
}
\seealso{
Package \code{diagis} in CRAN, which provides functions for computing weighted
summary statistics.
}