| fbase2.gamma.log.log | R Documentation |
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.
fbase2.gamma.log.log(u, v, y, fgh = 2)
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 |
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 |
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 |
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.
Alireza S. Mahani, Mansour T.A. Sharabiani
regfac.expand.2par
## 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)
## End(Not run)