OveDisPoiReg: Regression part

In fblues/Test: Prediction region based on over-dispersed Poisson Regression

Description

Function handles the regression part

Usage

OveDisPoiReg(Y, X, X0 = NULL, xival = NULL, alpha = 0.05, tol = 10^(-5))

Arguments

Y	data
X	data
X0	prediction
xival	over-dispersion parameter, xi
alpha	alpha
tol	approximation tolerance

Details

Okay

Value

list

Author(s)

E. A. Pena

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (Y, X, X0 = NULL, xival = NULL, alpha = 0.05, tol = 10^(-5)) 
{
    p = dim(X)[2]
    n = dim(X)[1]
    if (is.null(X0)) {
        X0 = X
    }
    m = dim(X0)[1]
    Xmean = apply(X, 2, "mean")
    Xstd = apply(X, 2, "sd")
    XS = X
    X0S = X0
    for (j in 2:p) {
        XS[, j] = (X[, j] - Xmean[j])/Xstd[j]
        X0S[, j] = (X0[, j] - Xmean[j])/Xstd[j]
    }
    outglm = glm(Y ~ XS[, -1], family = poisson(link = log))
    outglmsumm = summary(outglm)
    theta = outglm$coefficients
    if (is.null(xival)) {
        xi = NewtRaphXi(1, theta, Y, XS, tol)$xi
    }
    else {
        xi = Inf
    }
    XiCov = NewtRaph(theta, xi, Y, XS, tol)$XiCov
    XiCov11 = XiCov[1:p, 1:p]
    rho0 = exp(X0S %*% theta)
    zcrit = qnorm(1 - alpha/2)
    PE = rep(0, m)
    PE2 = rep(0, m)
    for (i in 1:m) {
        PE[i] = zcrit * sqrt(rho0[i] + rho0[i] * (1 + rho0[i])/xi + 
            (1/n) * (rho0[i]^2) * (t(X0S[i, ]) %*% XiCov11) %*% 
                X0S[i, ])
    }
    PILow = rho0 - PE
    PIHig = rho0 + PE
    for (i in 1:m) {
        PILow[i] = max(0, PILow[i])
    }
    return(list(outglm = outglm, outglmsumm = outglmsumm, theta = theta, 
        xi = xi, XiCov = XiCov, Xmean = Xmean, Xstd = Xstd, FitPred = data.frame(rho0, 
            PILow, PIHig, X0)))
  }

fblues/Test documentation built on June 27, 2020, 12:18 a.m.