RegressionFactory (version 0.7.4)

fbase2.gamma.log.log: Double-Parameter Base Log-likelihood Function for Gamma GLM

Description

Vectorized, double-parameter base log-likelihood functions for Gamma GLM. The link functions map the mean and dispersion parameter to linear predictors. The base function can be supplied to the expander function regfac.expand.2par in order to obtain the full, high-dimensional log-likleihood and its derivatives.

Usage

fbase2.gamma.log.log(u, v, y, fgh = 2)

Arguments

u

First parameter of the base log-likelihood function (usually the result of applying a link function to distribution mean). This parameter is intended to be projected onto a high-dimensional space using the familiar regression transformation of u <- X%*%beta. In the typical use-case where the caller is regfac.expand.2par, a vector of values are supplied, and return objects will have the same length as u or v.

v

Second parameter of the base log-likelihood function (usually the result of applying a link function to distribution dispersion parameter). This parameter is intended to be projected onto a high-dimensional space using the familiar regression transformation of v <- Z%*%gamma. In the typical use-case where the caller is regfac.expand.2par, a vector of values are supplied, and return objects will have the same length as u or v.

y

Fixed slot of the base distribution, corresponding to the response variable in the regression model.

fgh

Integer with possible values 0,1,2. If fgh=0, the function only calculates and returns the log-likelihood vector and no derivatives. If fgh=1, it returns the log-likelihood and its first derivative in a list. If fgh=2, it returns the log-likelihood, as well as its first and second derivatives in a list.

Value

If fgh==0, the function returns -exp(-v)*(u + y*exp(-u) + v - log(y)) - log(y) - log(gamma(exp(-v))). If fgh==1, a list is returned with elements f and g, where f is the same object as in fgh==0 and g is a matrix of dimensions length(u)-by-2, with first column being the derivative of the above expression with respect to u, and the second column being the derivative of the above expression with respect to v. If fgh==2, the list will include an element named h, which is a matrix of dimensions length(u)-by-3, with the first column being the second derivative of f with respect to u, the second column being the second derivative of f with respect to v, and the third column is the cross-derivative term, i.e. the derivative of f with respect to u and v.

See Also

regfac.expand.2par

Examples

# NOT RUN {
# we use this library for univariate slice sampling
# of multivariate distributions
library(MfUSampler)
library(dglm)

# simulating data according to assumed generative model
# we assume log link functions for both mean and dispersion
# given variance function V(mu) = mu^2, we have:
# log(mu) = X%*%beta
# log(phi) = X%*%gamma
N <- 10000
K <- 5
X <- cbind(1,matrix(runif(N*(K-1), min=-0.5, max=+0.5), ncol=K-1))
beta <- runif(K, min=0.0, max=+1.0)
gamma <- runif(K, min=0.0, max=+1.0)
shape.vec <- 1 / exp(X%*%gamma)
rate.vec <- 1 / exp(X%*%gamma + X%*%beta)
y <- rgamma(N, shape = shape.vec, rate = rate.vec)
# implied dispersion:
dispersion.vec <- 1 / shape.vec

# model estimation using dglm package
reg.dglm <- dglm(y~X-1, dformula = ~X-1, family=Gamma(link="log"), dlink = "log")
beta.dglm <- reg.dglm$coefficients
gamma.dglm <- reg.dglm$dispersion.fit$coefficients

# model estimation using RegressionFactory
# (with univariate slice sampling)
# defining the log-likelihood using the expander framework
# assumng same covariates for both slots, hence we set Z=X
# slice sampler does not need derivatives, hence we set fgh=0
loglike.gamma <- function(coeff, X, y) {
  regfac.expand.2par(coeff, X=X, Z=X, y=y, fbase2=fbase2.gamma.log.log, fgh=0)
}
nsmp <- 100
coeff.smp <- array(NA, dim=c(nsmp, 2*K)) 
coeff.tmp <- rep(0.1, 2*K)
for (n in 1:nsmp) {
  coeff.tmp <- MfU.Sample(coeff.tmp, f=loglike.gamma, X=X, y=y)
  coeff.smp[n,] <- coeff.tmp
}
beta.slice <- colMeans(coeff.smp[(nsmp/2+1):nsmp, 1:K])
gamma.slice <- colMeans(coeff.smp[(nsmp/2+1):nsmp, K+1:K])

# compare results
cbind(beta.dglm, beta.slice)
cbind(gamma.dglm, gamma.slice)

# }