R: Double-Parameter Base Log-likelihood Function for Gaussian...

fbase2.gaussian.identity.log

R Documentation

Double-Parameter Base Log-likelihood Function for Gaussian GLM

Description

Vectorized, double-parameter base log-likelihood functions for Gaussian GLM. The link functions map the mean and variance 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.gaussian.identity.log(u, v, y, fgh = 2)

Arguments

`u`	First parameter of the base log-likelihood function (usually the 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.1par`, a vector of values are supplied, and return objects will have the same length as `u`.
`v`	Second parameter of the base log-likelihood function (usually the mean). 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.1par`, a vector of values are supplied, and return objects will have the same length as `u`.
`y`	Fixed slot of the base distribution, corresponding to the response variable in the regression model. For `binomial` family, it must be an integer vector with values between `0` and `n`.
`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 -0.5*(v + exp(-v)*(u-y)*(u-y)) (ignoring additive terms that are independent of u,v). It will therefore be of the same length as u. If fgh==1, a list is returned with elements f and g, where f is the same object as in fgh==1 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.

Author(s)

Alireza S. Mahani, Mansour T.A. Sharabiani

Examples

## Not run: 
library(sns)
library(MfUSampler)
library(dglm)

# defining log-likelihood function
# vd==FALSE leads to constant-dispersion model (ordinary linear regression)
# while vd==TRUE produces varying-dispersion model
loglike.linreg <- function(coeff, X, y, fgh, block.diag = F, vd = F) {
  if (vd) regfac.expand.2par(coeff = coeff, X = X, Z = X, y = y
    , fbase2 = fbase2.gaussian.identity.log, fgh = fgh, block.diag = block.diag)
  else regfac.expand.2par(coeff = coeff, X = X, y = y
    , fbase2 = fbase2.gaussian.identity.log, fgh = fgh, block.diag = block.diag)
}

# simulating data according to generative model
N <- 1000 # number of observations
K <- 5 # number of covariates
X <- matrix(runif(N*K, min=-0.5, max=+0.5), ncol=K)
beta <- runif(K, min=-0.5, max=+0.5)
gamma <- runif(K, min=-0.5, max=+0.5)
mean.vec <- X%*%beta
sd.vec <- exp(X%*%gamma)
y <- rnorm(N, mean.vec, sd.vec)

# constant-dispersion model
# estimation using glm
est.glm <- lm(y~X-1)
beta.glm <- est.glm$coefficients
sigma.glm <- summary(est.glm)$sigma
# estimation using RegressionFactory
# (we set rnd=F in sns to allow for better comparison with glm)
nsmp <- 20
coeff.smp <- array(NA, dim=c(nsmp, K+1)) 
coeff.tmp <- rep(0, K+1)
for (n in 1:nsmp) {
  coeff.tmp <- sns(coeff.tmp, fghEval=loglike.linreg
    , X=X, y=y, fgh=2, block.diag = F, vd = F, rnd = F)
  coeff.smp[n,] <- coeff.tmp
}
beta.regfac.cd <- colMeans(coeff.smp[(nsmp/2+1):nsmp, 1:K])
sigma.regfac.cd <- sqrt(exp(mean(coeff.smp[(nsmp/2+1):nsmp, K+1])))
# comparing glm and RegressionFactory results
# beta's must match exactly between glm and RegressionFactory
cbind(beta, beta.glm, beta.regfac.cd)
# sigma's won't match exactly
cbind(mean(sd.vec), sigma.glm, sigma.regfac.cd)

# varying-dispersion model
# estimation using dglm
est.dglm <- dglm(y~X-1, dformula = ~X-1, family = "gaussian", dlink = "log")
beta.dglm <- est.dglm$coefficients
gamma.dglm <- est.dglm$dispersion.fit$coefficients
# estimation using RegressionFactory
coeff.smp <- array(NA, dim=c(nsmp, 2*K)) 
coeff.tmp <- rep(0, 2*K)
for (n in 1:nsmp) {
  coeff.tmp <- sns(coeff.tmp, fghEval=loglike.linreg
    , X=X, y=y, fgh=2, block.diag = F, vd = T, rnd = F)
  coeff.smp[n,] <- coeff.tmp
}
beta.regfac.vd <- colMeans(coeff.smp[(nsmp/2+1):nsmp, 1:K])
gamma.regfac.vd <- colMeans(coeff.smp[(nsmp/2+1):nsmp, K+1:K])
# comparing dglm and RegressionFactory results
# neither beta's nor gamma's will match exactly
cbind(beta, beta.dglm, beta.regfac.vd)
cbind(gamma, gamma.dglm, gamma.regfac.vd)


## End(Not run)