BMA: Bayesian Model Averaging

"as.bic.glm" <-
function (g, ...) 
{
    UseMethod("as.bic.glm")
}


summary.glib <-
function (object, n.models = 5, digits = max(3, getOption("digits") - 
    3), conditional = FALSE, index.phi=1, ...) 
{
   x<- object
   gobj<- as.bic.glm.glib(x, index.phi=index.phi) 
   summary(gobj, n.models=n.models, digits=digits, conditional=conditional, display.dropped=FALSE, ...)

}

"as.bic.glm.glib" <-
function (g, index.phi = 1, ...) 
{
    probne0 <- 100*(1-g$posterior$prob0[, index.phi])
    postmean <- c(NA, g$posterior$mean[, index.phi])
    postsd <- c(NA, g$posterior$sd[, index.phi])
    size <- g$bf$npar
    bic <- g$bf$twologB10[, index.phi]
    postprob <- g$bf$postprob[, index.phi]
    deviance <- g$bf$deviance
    disp <- 1
    which <- g$models == 1
    family <- g$inputs$error
    if (!is.null(g$glm.out)) 
        reduced <- g$glm.out$reduced
    else reduced <- FALSE
    if (!is.null(g$glm.out)) 
        dropped <- g$glm.out$dropped
    else dropped <- NULL
    call <- g$call
    n.models <- length(postprob)
    n.vars <- length(probne0)
    factor.type <- FALSE
    x <- g$inputs$x
    y <- g$input$y

    if (is.null(colnames(x)))
	colnames(x)<- paste("X",1:ncol(x),sep="")

    nms <- c("Intercept", colnames(x))
    assign <- list()
    for (i in 1:length(nms)) assign[[i]] <- i
    names(assign) <- nms
    output.names <- list()
    for (i in 2:length(nms)) output.names[[i - 1]] <- NA
    names(output.names) <- nms[-1]
    namesx <- colnames(x)
    if (!is.null(g$glm.out)) 
        design <- g$glm.out$design
    else design <- NULL
    label <- rep(NA, times = n.models)
    for (i in 1:n.models) label[i] <- paste(namesx[which[i, ]], 
        sep = "", collapse = ",")
    mle <- matrix(rep(0, times = n.models * (n.vars + 1)), ncol = n.vars + 
        1)
    for (i in 1:n.models) mle[i, c(FALSE, which[i, ])] <- rbind(g$posterior.bymodel$mean[[i]][-1, 
        index.phi])
    se <- matrix(rep(0, times = n.models * (n.vars + 1)), ncol = n.vars + 
        1)
    for (i in 1:n.models) se[i, c(FALSE, which[i, ])] <- rbind(g$posterior.bymodel$sd[[i]][-1, 
        index.phi])

    mle[,1]<- rep(NA, times=nrow(mle))
    se[,1]<- rep(NA, times=nrow(se))

    nests <- n.vars + 1
    CSDbi <- rep(NA, n.vars)
    CEbi <- CSDbi
    for (i in (1:ncol(x))) {
        sel <- which[, i]
        if (sum(sel) > 0) {
            cpp <- rbind(postprob[sel]/sum(postprob[sel]))
            CEbi[i + 1] <- as.numeric(cpp %*% mle[sel, i + 1])
            CSDbi[i + 1] <- sqrt(cpp %*% (se[sel, i + 1]^2) + 
                cpp %*% ((mle[sel, i + 1] - CEbi[i + 1])^2))
        }
    }
    prior.model.weights <- NULL
    prior.param <- NULL
    result <- list(postprob = postprob, label = label, deviance = deviance, 
        size = size, bic = bic, prior.param = prior.param, prior.model.weights = prior.model.weights, 
        family = family, disp = disp, which = which, probne0 = probne0, 
        postmean = postmean, postsd = postsd, condpostmean = CEbi, 
        condpostsd = CSDbi, mle = mle, se = se, namesx = namesx, 
        reduced = reduced, dropped = dropped, call = call, n.models = n.models, 
        n.vars = n.vars, nests = nests, output.names = output.names, 
        assign = assign, factor.type = factor.type, design = design, 
        x = x, y = y)
    class(result) <- "bic.glm"
    return(result)
}

"glib" <-
function (x, ...) 
UseMethod("glib")


"glib.data.frame" <-
function (x, y, n = rep(1, nrow(x)), error = "poisson", link = "log", 
    scale = 1, models = NULL, phi = c(1, 1.65, 5), psi = 1, nu = 0, 
    pmw = rep(1, nrow(models)), glimest = TRUE, glimvar = FALSE, output.priorvar = FALSE, 
    post.bymodel = TRUE, output.postvar = FALSE, priormean = NULL, priorvar = NULL, 
    nbest = 150, call = NULL, ...) 
{
    glim.pmean.A <- function(x, y, n, error, link, scale, nu, 
        prior.spec) {
        glimot <- glim(x, y, n, error = error, link = link, scale = scale)
        coef <- matrix(glimot$coef, ncol = 1)
        eta <- cbind(rep(1, nrow(x)), x) %*% coef
        if (link == "identity") {
            mu <- eta
            deta <- rep(1,length(mu))
        }
        if (link == "log") {
            mu <- exp(eta)
            deta <- 1/mu
        }
        if (link == "logit") {
            mu <- exp(eta)/(1 + exp(eta)) * n
            deta <- 1/mu + 1/(n - mu)
        }
        if (link == "probit") {
            mu <- pnorm(eta)
            deta <- 1/(dnorm(qnorm(mu)))
        }
        if (link == "sqrt") {
            mu <- eta^2
            deta <- 1/(2 * sqrt(mu))
        }
        if (link == "inverse") {
            mu <- 1/eta
            deta <- -1/(mu^2)
        }
        if (link == "loglog") {
            mu <- 1 - exp(-exp(eta))
            deta <- -1/((1 - mu) * log(1 - mu))
        }
        if (error == "gaussian") 
            v <- rep(1,length(mu))
        if (error == "poisson") 
            v <- mu
        if (error == "binomial") 
            v <- mu * (1 - mu/n)
        if (error == "gamma") 
            v <- mu^2
        if (error == "inverse gaussian") 
            v <- mu^3
        w <- 1/(deta^2 * v)
        w <- w/sum(w)
        z <- eta + (y - mu) * deta
        xbar <- t(x) %*% w
        xbar2 <- t(x^2) %*% w
        s2x <- xbar2 - xbar^2
        zbar <- as.numeric(weighted.mean(z, w))
        zbar2 <- t(z^2) %*% w
        if (error == "binomial") 
            s2z <- weighted.mean((eta - zbar)^2 + 2 * (eta - 
                zbar) * (deta * n) * (y - mu)/n + (deta * n)^2 * 
                (y - 2 * mu * y/n + mu^2/n)/n, w)
        else s2z <- zbar2 - zbar^2
        varz <- as.numeric(s2z)
        if (is.numeric(prior.spec$mean)) 
            pmean <- prior.spec$mean
        else pmean <- c((zbar + nu * sqrt(varz)), rep(0, ncol(x)))
        if (is.numeric(prior.spec$var)) 
            A <- prior.spec$var
        else {
            A <- diag(as.vector(1/sqrt(s2x)), nrow = length(s2x))
            A <- rbind(t(-xbar/sqrt(s2x)), A)
            Acol1 <- c(1, rep(0, ncol(x)))
            A <- sqrt(varz) * cbind(Acol1, A)
        }
        list(x = x, y = y, n = n, error = error, link = link, 
            pmean = pmean, A = A, scale = glimot$scale, varz = varz)
    }
    glim.pvar <- function(model, phi = 1, psi = 1, Aot, prior.spec) {
        model1 <- c(1, model + 1)
        if (is.numeric(prior.spec$var)) {
            sigma <- Aot$A
            sig11 <- sigma[model1, model1]
            sig12 <- sigma[model1, -model1]
            sig21 <- sigma[-model1, model1]
            sig22 <- sigma[-model1, -model1]
            if (length(model1) != ncol(Aot$A)) 
                pvar <- sig11 - sig12 %*% solve(sig22) %*% sig21
            else pvar <- sigma
        }
        else {
            A <- Aot$A[model1, model1]
            p <- length(model)
            B <- diag(c(psi^2, rep(phi^2, p)), nrow = (p + 1))
            pvar <- A %*% B %*% t(A)
        }
        list(model = model, phi = phi, psi = psi, pvar = pvar)
    }
    E1 <- function(model, phi, psi, glimot1, Aot, prior.spec) {
        V <- glimot1$var
        A <- solve(V)
        logdetV <- sum(log(eigen(V)$values))
        pmean <- matrix(Aot$pmean[c(1, model + 1)], ncol = 1)
        thetahat <- matrix(glimot1$coef, ncol = 1)
        nphi <- length(phi)
        app1 <- rep(0, nphi)
        npar <- length(model) + 1
        pvar <- array(rep(0, npar * npar * nphi), dim = c(npar, 
            npar, nphi))
        for (j in (1:nphi)) {
            phij <- phi[j]
            pvar[, , j] <- glim.pvar(model, phij, psi, Aot, prior.spec)$pvar
            B <- solve(pvar[, , j])
            C <- solve(A + B)
            I <- diag(rep(1, length(model) + 1), nrow = length(model) + 
                1)
            logdetW <- sum(log(eigen(pvar[, , j])$values))
            F1 <- B %*% C %*% A %*% C %*% B + (t(I - C %*% B)) %*% 
                B %*% (I - C %*% B)
            app1[j] <- -(t(thetahat - pmean)) %*% F1 %*% (thetahat - 
                pmean) - logdetW - sum(log(eigen(A + B)$values))
        }
        list(model = model, phi = phi, psi = psi, app1 = app1)
    }
    E0 <- function(psi = 1, glimot0, Aot) {
        V <- glimot0$var
        A <- 1/V
        logV <- log(V)
        pmean <- Aot$pmean[1]
        thetahat <- glimot0$coef
        pvar <- Aot$varz * psi^2
        B <- 1/pvar
        C <- 1/(A + B)
        logW <- log(pvar)
        F0 <- B^2 * C^2 * A + (1 - C * B)^2 * B
        app1 <- F0 * (thetahat - pmean)^2 - logW - log(A + B)
        list(psi = psi, app1 = app1)
    }
    glim <- function(x, y, n, error = "gaussian", link = "identity", 
        wt = 1, resid = "none", init, intercept = TRUE, scale, offset = 0, 
        sequence, eps = 1e-04, iter.max = 10) {
        error.int <- charmatch(error, c("gaussian", "poisson", 
            "binomial", "gamma", "inverse gaussian"))
        if (is.na(error.int)) 
            stop("Invalid error type")
        else if (error.int == 0) 
            stop("Ambiguous error type")
        resid.int <- charmatch(resid, c("none", "pearson", "deviance"))
        if (is.na(resid.int)) 
            stop("Invalid residual type")
        if (!missing(scale) && !is.numeric(scale)) {
            temp <- charmatch(scale, c("pearson", "deviance"))
            if (is.na(temp)) 
                stop("Invalid scale option")
        }
        if (!(is.numeric(y) && is.numeric(x))) 
            stop("Invalid y or x values")
        x <- as.matrix(x)
        nn <- nrow(x)
        nvar <- ncol(x)
        if (length(y) != nn) 
            stop("Dimensions of x and y don't match")
        if (!missing(wt)) {
            if (length(wt) != nn) 
                stop("Dimensions of x and wt don't match")
            if (any(wt < 0, na.rm = TRUE)) 
                stop("Weights must be >=0")
        }
        if (!missing(offset) && (length(offset) != nn)) 
            stop("Dimensions of x and offset don't match")
        if (error.int != 3) 
            n <- rep(1, nn)
        else {
            if (missing(n)) 
                stop("Binomial fits require the n vector")
            else if (length(n) != nn) 
                stop("Length of n vector is incorrect")
            n <- as.vector(n)
        }
        nomiss <- !(is.na(y) | is.na(wt) | is.na(n) | is.na(x %*% 
            rep(1, nvar)))
        if (sum(nomiss) < nn) {
            warning(paste(nn - sum(nomiss), "deleted due to missing values"))
            x <- x[nomiss, , drop = FALSE]
            y <- y[nomiss]
            n <- n[nomiss]
            if (!missing(wt)) 
                wt <- wt[nomiss]
            if (!missing(offset)) 
                offset <- offset[nomiss]
            nn <- sum(nomiss)
        }
        if (missing(sequence)) 
            sequence <- c(0, nvar)
        else {
            if (max(sequence) > nvar) 
                stop("Invalid sequence argument")
            if (min(sequence) < 0) 
                stop("Invalid sequence argument")
            if (any(diff(sequence) <= 0)) 
                stop("Invalid sequence argument")
        }
        xn <- dimnames(x)[[2]]
        if (is.null(xn)) 
            xn <- paste("X", 1:nvar, sep = "")
        if (intercept) {
            x <- cbind(rep(1, nn), x)
            xn <- c("Intercept", xn)
            nvar <- nvar + 1
            sequence <- sequence + 1
        }
        dimnames(x) <- list(dimnames(x)[[1]], xn)
        if (error.int != 1 && any(y < 0)) 
            stop("Illegal y values")
        if (error.int >= 4 && any(y == 0)) 
            stop("Illegal y values")
        if (error.int == 3) {
            y <- y/n
            if (any(y > 1)) 
                stop("Illegal y values")
        }
        var <- switch(error.int, function(x, n) rep(1, length(x)), 
            function(x, n) x, function(x, n) (x - x^2)/n, function(x, 
                n) x^2, function(x, n) x^3)
        dev <- switch(error.int, function(y, mu, n) (y - mu)^2, 
            function(y, mu, n) 2 * (y * log(ifelse(y == 0, 1, 
                y/mu)) - (y - mu)), function(y, mu, n) 2 * n * 
                (y * log(ifelse(y == 0, 1, y/mu)) + (1 - y) * 
                  log(ifelse(y == 1, 1, (1 - y)/(1 - mu)))), 
            function(y, mu, n) -2 * (log(y/mu) - (y - mu)/mu), 
            function(y, mu, n) (y - mu)^2/(y * mu^2))
        link.int <- charmatch(link, c("identity", "log", "logit", 
            "sqrt", "inverse", "probit", "loglog"))
        if (is.na(link.int)) 
            stop("Invalid link type")
        else if (link.int == 0) 
            stop("Ambiguous link type")
        f <- switch(link.int, function(x) x, function(x) exp(x), 
            function(x) {
                z <- exp(ifelse(x > 80, 80, ifelse(x < -80, -80, 
                  x)))
                z/(z + 1)
            }, function(x) x^2, function(x) 1/x, pnorm, function(x) {
                z <- exp(ifelse(x > 80, 80, ifelse(x < -80, -80, 
                  x)))
                1 - exp(-z)
            })
        deriv <- switch(link.int, function(x) rep(1, length(x)), 
            function(x) exp(x), function(x) exp(x)/((1 + exp(x))^2), 
            function(x) 2 * x, function(x) -1/(x^2), function(x) 0.3989422 * 
                exp(-0.5 * x^2), function(x) exp(x) * exp(-exp(x)))
        temp.y <- switch(error.int, y, y + 0.01, (n * y + 0.5)/(n + 
            1), y, y)
        eta <- switch(link.int, temp.y, log(temp.y), log(temp.y/(1 - 
            temp.y)), sqrt(temp.y), 1/temp.y, qnorm(temp.y), 
            log(-log(1 - temp.y)))
        if (sum(is.na(eta)) > 0) 
            stop("Illegal value(s) of y for this link function")
        if (missing(init)) {
            tempd <- deriv(eta)
            z <- eta + (y - temp.y)/tempd - offset
            w <- sqrt((wt * tempd^2)/var(temp.y, n))
            fit <- qr(x * w)
            if (fit$rank < nvar) 
                stop("X matrix is not full rank")
            init.qr <- fit$qr[1:nvar, ]
            init.qty <- qr.qty(fit, z * w)[1:nvar]
        }
        else if (length(init) != nvar) 
            stop("Argument 'init' is the wrong length")
        if (!missing(wt)) 
            nn <- sum(wt > 0)
        deviance <- df <- NULL
        if (sequence[1] == 1 && intercept && missing(offset)) {
            if (missing(wt)) 
                yhat <- sum(y * n)/sum(n)
            else yhat <- sum((y * n * wt)[wt > 0])/sum((n * wt)[wt > 
                0])
            if ((yhat > 1) && (link.int == 3 || link.int == 6 || 
                link.int == 7)) 
                yhat <- 1
            if ((yhat < 0) && link.int > 1 && link.int != 5) 
                yhat <- 0
            deviance <- sum(dev(y, yhat, n) * wt)
            df <- nn - 1
            sequence <- sequence[-1]
        }
        else if (sequence[1] == 0) {
            deviance <- sum(dev(y, f(offset), n) * wt)
            df <- nn
            sequence <- sequence[-1]
        }
        for (modnum in sequence) {
            model <- 1:modnum
            if (missing(init)) 
                coef <- backsolve(init.qr, init.qty, k = modnum)[model]
            else coef <- init[model]
            tempx <- as.matrix(x[, model, drop = FALSE])
            nvar <- length(model)
            eta <- c(tempx %*% coef)
            yhat <- f(eta + offset)
            olddev <- sum(dev(y, yhat, n) * wt)
            fini <- FALSE
            for (i in seq(length = iter.max)) {
                tempd <- deriv(eta + offset)
                if (any(is.na(yhat) | is.na(tempd))) {
                  warning("Coef vector is diverging, iteration truncated")
                  break
                }
                leave.in <- ((1 + tempd) != 1)
                z <- eta + (y - yhat)/tempd
                w <- sqrt((wt * tempd^2)/var(yhat, n))
                fit <- qr(tempx[leave.in, ] * w[leave.in])
                if (fit$rank < nvar) {
                  warning(paste("Weighted X matrix no longer of full rank at iteration", 
                    i))
                  break
                }
                coef <- qr.coef(fit, (z * w)[leave.in])
                eta <- c(tempx %*% coef)
                yhat <- f(eta + offset)
                newdev <- sum(dev(y, yhat, n) * wt)
                if (abs((olddev - newdev)/(newdev + 1)) < eps) {
                  fini <- TRUE
                  break
                }
                else olddev <- newdev
            }
            if (fini == FALSE) 
                warning(paste("Model with", length(model), "variables did not converge"))
            deviance <- c(deviance, newdev)
            df <- c(df, nn - nvar)
        }
        coef <- as.vector(coef)
        names(coef) <- dimnames(tempx)[[2]]
        if (missing(scale)) 
            scale <- switch(error.int, newdev/(nn - nvar), 1, 
                1, newdev/(nn - nvar), newdev/(nn - nvar))
        else if (!is.numeric(scale)) 
            scale <- switch(charmatch(scale, c("pearson", "deviance")), 
                sum((y - yhat)^2/(var(yhat) * (nn - nvar))), 
                newdev/(nn - nvar))
        varmat <- backsolve(fit$qr[1:nvar,,drop=F], diag(nvar))
        varmat <- varmat %*% t(varmat)
        temp <- list(coef = coef, var = varmat * scale, deviance = deviance, 
            df = df, scale = scale, intercept = intercept)
        if (resid.int == 1) 
            temp
        else if (resid.int == 3) {
            resid <- rep(NA, length(nomiss))
            resid.good <- sqrt(dev(y, yhat, n))
            resid[nomiss] <- ifelse(y > yhat, resid.good, -resid.good)
            c(temp, list(resid = resid))
        }
        else {
            resid <- rep(NA, length(nomiss))
            resid[nomiss] <- (y - yhat)/ifelse(y == yhat, 1, 
                sqrt(var(yhat, n)))
            c(temp, list(resid = resid))
        }
    }


    ##### start of function proper #####

    if (is.null(call)) 
        cl <- match.call()
    else cl <- call
    glm.out <- NULL
    prior.spec <- list(mean = FALSE, var = FALSE)
    if (!is.null(priormean)) 
        prior.spec$mean <- priormean
    if (!is.null(priorvar)) {
        prior.spec$var <- priorvar
        phi <- 1
    }
    if (is.null(models)) {
        glm.out <- bic.glm(x = x, y = y, glm.family = error, 
            nbest = nbest, dispersion = NULL, occam.window = FALSE, strict = FALSE, 
            factor.type = FALSE, ...)

        tempdata <- data.frame(glm.out$x, glm.out$y)
        mm <- model.matrix(terms(y ~ ., data = tempdata), data = tempdata)
        df <- data.frame(mm)
        df <- df[, -ncol(df)]
        df <- df[, -1]
        xin <- as.matrix(df)
        y <- glm.out$y
        models <- glm.out$which + 0
        whichin <- matrix(rep(NA, times = nrow(models) * ncol(xin)), 
            ncol = ncol(xin))
        for (i in 1:ncol(models)) whichin[, glm.out$assign[[i + 
            1]] - 1] <- models[, i]
        models <- whichin
        family <- glm.out$family
        if (is.character(family)) 
            family <- get(family, mode = "function", envir = parent.frame())
        if (is.function(family)) 
            family <- family()
        error <- family$family
        link <- family$link
	x<- xin

    }
    if (is.data.frame(x)) {
        nx <- names(x)
        x <- as.matrix(x)
        dimnames(x) <- list(NULL, nx)
    }
    glimot0 <- glim(rep(1, length(y)), y, n, error = error, link = link, 
        intercept = FALSE, scale = scale)
    Aot <- glim.pmean.A(x, y, n, error, link, scale, nu, prior.spec)
    scale <- Aot$scale
    pmean1 <- Aot$pmean
    E0ot <- E0(psi, glimot0, Aot)
    nmodel <- nrow(models)
    nphi <- length(phi)
    chi2 <- rep(0, nmodel)
    npar <- rep(0, nmodel)
    app1 <- matrix(rep(0, nmodel * nphi), ncol = nphi, nrow = nmodel)
    glim.coef <- as.list(rep(0, nmodel))
    glim.var <- as.list(rep(0, nmodel))
    glim.se <- as.list(rep(0, nmodel))
    prior.var <- as.list(rep(0, nmodel))
    postbym.mean <- as.list(rep(0, nmodel))
    postbym.sd <- as.list(rep(0, nmodel))
    postbym.var <- as.list(rep(0, nmodel))
    deviance <- rep(0, nmodel)
    df <- rep(0, nmodel)
    prob0 <- matrix(rep(0, (ncol(x) * nphi)), ncol = nphi)
    post.mean <- matrix(rep(0, (ncol(x) * nphi)), ncol = nphi)
    post.sd <- matrix(rep(0, (ncol(x) * nphi)), ncol = nphi)
    postprob <- matrix(rep(0, (nmodel * nphi)), ncol = nphi)
    prior.var <- as.list(rep(0, nmodel))
    for (i in (1:nmodel)) {
        if (sum(models[i, ]) == 0) {
            npar[i] <- 1
            prior.var[[i]] <- array(rep(0, nphi), dim = c(1, 
                1, nphi))
            for (j in (1:nphi)) prior.var[[i]][, , j] <- psi^2 * 
                Aot$varz
        }
        else {
            model <- (1:ncol(x))[models[i, ] == 1]
            npar[i] <- length(model) + 1
            prior.var[[i]] <- array(rep(0, npar[i] * npar[i] * 
                nphi), dim = c(npar[i], npar[i], nphi))
            for (j in (1:nphi)) prior.var[[i]][, , j] <- glim.pvar(model, 
                phi[j], psi, Aot, prior.spec)$pvar
        }
    }
    for (i in (1:nmodel)) 
    {
        if (sum(models[i, ]) == 0) 
	{
            chi2[i] <- 0
            npar[i] <- 1
            app1[i, ] <- 0
            betahat <- glimot0$coef
            glim.coef[[i]] <- betahat
            V <- glimot0$var
            glim.var[[i]] <- V
            glim.se[[i]] <- sqrt(diag(glim.var[[i]]))
            deviance[i] <- glimot0$deviance[2]
            df[i] <- glimot0$df[2]
        }
        else 
	{
            model <- (1:ncol(x))[models[i, ] == 1]
            glimot1 <- glim(x[, model], y, n = n, error = error, 
                link = link, scale = scale)
            chi2[i] <- (glimot0$deviance[2] - glimot1$deviance[2])/scale
            npar[i] <- sum(models[i, ]) + 1
            E1ot <- E1(model, phi, psi, glimot1, Aot, prior.spec)
            app1[i, ] <- chi2[i] + (E1ot$app1 - E0ot$app1)
            betahat <- glimot1$coef
            glim.coef[[i]] <- betahat
            V <- glimot1$var
            glim.var[[i]] <- V
            glim.se[[i]] <- sqrt(diag(glim.var[[i]]))
            deviance[i] <- glimot1$deviance[2]
            df[i] <- glimot1$df[2]
        }


        if (sum(models[i, ]) == 0) 
	{
            prior.mean <- pmean1[1]
            postbym.mean[[i]] <- matrix(rep(0, nphi), ncol = nphi)
            postbym.sd[[i]] <- matrix(rep(0, nphi), ncol = nphi)
            postbym.var[[i]] <- array(rep(0, nphi), dim = c(1, 
                1, nphi))
            for (j in (1:nphi)) 
	    {
                W <- prior.var[[i]][, , j]
                postbym.mean[[i]][, j] <- betahat - V %*% solve(V + 
                  W) %*% (betahat - prior.mean)
                postbym.var[[i]][, , j] <- solve(solve(V) + solve(W))
                postbym.sd[[i]][, j] <- sqrt(postbym.var[[i]][, 
                  , j])
            }
        }
        else 
	{
            prior.mean <- c(pmean1[1], rep(0, (npar[i] - 1)))
            postbym.mean[[i]] <- matrix(rep(0, (npar[i] * nphi)), 
                ncol = nphi)
            postbym.sd[[i]] <- matrix(rep(0, (npar[i] * nphi)), 
                ncol = nphi)
            postbym.var[[i]] <- array(rep(0, (npar[i] * npar[i] * 
                nphi)), dim = c(npar[i], npar[i], nphi))
            for (j in (1:nphi)) 
	    {
                W <- prior.var[[i]][, , j]
                postbym.mean[[i]][, j] <- betahat - V %*% solve(V + 
                  W) %*% (betahat - prior.mean)
                postbym.var[[i]][, , j] <- solve(solve(V) + solve(W))
                postbym.sd[[i]][, j] <- sqrt(diag(postbym.var[[i]][, 
                  , j]))
            }
        }
    }

    for (j in (1:nphi)) 
    {
        tmp <- app1[, j] - max(app1[, j])
        tmp <- pmw * exp(0.5 * tmp)
        postprob[, j] <- tmp/sum(tmp)
        for (k in (1:ncol(x))) prob0[k, j] <- sum(postprob[models[, k] == 0, j])
 
        for (k in (1:ncol(x))) 
	{
            tmp1 <- 0
            tmp2 <- 0
            for (i in (1:nmodel)) 
	    {
                if (models[i, k] == 1) 
		{
                    pos <- sum(models[i, (1:k)]) + 1
                    tmp1 <- tmp1 + postbym.mean[[i]][pos, j] * postprob[i, j]
                    tmp2 <- tmp2 + (postbym.sd[[i]][pos, j]^2 + 
                            postbym.mean[[i]][pos, j]^2) * postprob[i, j]
                }
                post.mean[k, j] <- tmp1/(1 - prob0[k, j])
                post.sd[k, j] <- sqrt(tmp2/(1 - prob0[k, j]) - post.mean[k, j]^2)
            }
        }
    }
    inputs <- list(x = x, y = y, n = n, error = error, link = link, 
        models = models, phi = phi, psi = psi, nu = nu)
    if (is.numeric(prior.spec$var)) 
        inputs <- inputs[-(7:8)]
    if (is.numeric(prior.spec$mean)) 
        inputs <- inputs[-match("nu", names(inputs))]
    bf <- list(twologB10 = app1, postprob = postprob, deviance = deviance, 
        df = df, chi2 = chi2, npar = npar, scale = scale)
    posterior <- list(prob0 = prob0, mean = post.mean, sd = post.sd)
    if (glimest == TRUE) {
        if (glimvar == FALSE) 
            glim.est <- list(coef = glim.coef, se = glim.se)
        if (glimvar == TRUE) 
            glim.est <- list(coef = glim.coef, se = glim.se, 
                var = glim.var)
    }
    if (glimest == FALSE) 
        glim.est <- NULL
    if (post.bymodel == TRUE) {
        if (output.postvar == FALSE) 
            posterior.bymodel <- list(mean = postbym.mean, sd = postbym.sd)
        if (output.postvar == TRUE) 
            posterior.bymodel <- list(mean = postbym.mean, sd = postbym.sd, 
                var = postbym.var)
    }
    if (post.bymodel == FALSE) 
        posterior.bymodel <- NULL
    if (output.priorvar == FALSE) 
        prior <- list(mean = pmean1)
    if (output.priorvar == TRUE) 
        prior <- list(mean = pmean1, var = prior.var)
    result <- list(inputs = inputs, bf = bf, posterior = posterior, 
        glim.est = glim.est, posterior.bymodel = posterior.bymodel, 
        prior = prior, x = x, models = models, glm.out = glm.out, 
        call = cl)
    class(result) <- "glib"
    result
}




glib.matrix <- glib.data.frame

glib.bic.glm<- function (x, scale = 1, phi = 1, psi = 1, nu = 0, 
    glimest = TRUE, glimvar = FALSE, output.priorvar = FALSE, post.bymodel = TRUE, 
    output.postvar = FALSE, priormean = NULL, priorvar = NULL, call = NULL, ...) 
{
    if (is.null(call)) 
        cl <- match.call()
    else cl <- call
    tempdata <- data.frame(x$x, x$y)
    mm <- model.matrix(terms(y ~ ., data = tempdata), data = tempdata)
    df <- data.frame(mm)
    df <- df[, -ncol(df)]
    df <- df[, -1]
    xin <- as.matrix(df)
    y <- x$y
    models <- x$which + 0
    whichin <- matrix(rep(NA, times = nrow(models) * ncol(xin)), 
        ncol = ncol(xin))
    for (i in 1:ncol(models)) whichin[, x$assign[[i + 1]] - 1] <- models[, 
        i]
    models <- whichin
    family <- x$family
    if (is.character(family)) 
        family <- get(family, mode = "function", envir = parent.frame())
    if (is.function(family)) 
        family <- family()
    error <- family$family
    link <- family$link
    glib(x = xin, y = y, error = error, link = link, models = models, 
        scale = scale, phi = phi, psi = psi, nu = nu, glimest = glimest, 
        glimvar = glimvar, output.priorvar = output.priorvar, 
        post.bymodel = post.bymodel, output.postvar = output.postvar, 
        priormean = priormean, priorvar = priorvar, call = cl, ...)
}