glib <- function (x,y,n=rep(1,nrow(x)),error="poisson",link="log", scale=1,models=matrix(rep(1,ncol(as.matrix(x))),nrow=1),phi=c(1,1.65,5), psi=1,nu=0, pmw=rep(1,nrow(models)),glimest=T,glimvar=F, output.priorvar=F,post.bymodel=T,output.postvar=F, priormean,priorvar) { # Version 2.0: Revised November 20, 1996 by Chris Volinsky # (Version 1.1: May 18, 1993, by Adrian E. Raftery) # S-PLUS function to evaluate Bayes factors and account for model # uncertainty in generalized linear models. This also calculates # posterior distributions from a set of reference proper priors. # Revision (1.2) fixes bug concerning binomial data where n>1 and allows for prior # mean and variance specification by the user. # Developers: Adrian Raftery (raftery@stat.washington.edu) and Chris T. Volinsky # (volinsky@stat.washington.edu). # This software is not formally maintained, but we will be happy to # hear from people who have problems with it, although we cannot guarantee to # to solve them! # Permission is hereby given to StatLib to redistribute this software. # The software can be freely used for non-commercial purposes, and can # be freely distributed for non-commercial purposes only. # The copyright is retained by the developer. # Copyright 1993 Adrian E. Raftery # Copyright 1996 by Chris T. Volinsky and Adrian E. Raftery # Thanks to Mark Becker and Tim Futing Liao for pointing out bugs and to # Jon Wakefield for helpful suggestions. # INPUTS: # ------ # x,y,n,error,link,scale are as in "glim". Note that the x matrix # must be explicitly provided (this is the function "glim", not "glm"). # NOTE: Currently, link must be one of "log", "logit" or "loglog", # and error must be one "binomial" or "poisson". This will soon # be generalized to other links and error functions. # models is a nmodel x p matrix in which each row represents a model. # The corresponding entry in the row is 1 if that variable is included # in the model; 0 if not. Default: a vector of ones, corresponding to # the full model only. If the null model is to be included in the output, # a row of zeros must be included in models. # phi is the vector of phi values. Default: c(1,1.65,5) # psi is a scalar prior parameter. Default: 1. # nu is a scalar prior parameter. Default: 0 # pmw = prior model weights. These must be positive, but do not # have to sum to one. The prior model probabilities are given by # pmw/sum(pmw). Default:a nrow(models)-vector of ones. # The following inputs control the output: # glimest: If T, GLIM estimates and standard errors are output for each model # (default T). # glimvar: If T, GLIM variance matrices are output for each model (default F). # output.priorvar: If T, the prior variance is output for each model and value of phi. # (default F). # post.bymodel: If T, the posterior mean and sd are output for each model # and value of phi (default T). # output.postvar: If T, the posterior variance matrix is output for each model and # value of phi (default F) # The last two inputs allow the user to specify an arbitrary prior mean and variance matrix: # (Note that this is optional). # priormean: Optional user specified prior mean on the variables (including the intercept), # a (p+1) x 1 vector, where p=number of independent variables. # priorvar: Optional user specified prior variance matrix, a (p+1) x (p+1) matrix. # Default has the prior variance estimated as in Biometrika paper. # NOTE: The outputs controlled by glimvar, output.priorvar and output.postvar can take up # a lot of space, which is why these control parameters are F by default. # VALUE: # ----- # A list of six lists, as follows (some of these are output only if specified): # inputs: List echoing the inputs (x,y,n,error,link,models,phi,psi,nu) # bf: List containing the model comparison results, namely # twologB10= a nmodelxnphi matrix whose [i,j] element is 2logB10 for # model i against the null model with phi=phi[j]. My best Laplace # approximation (one-step Newton) is used. # postprob= a nmodelxnphi containing the posterior probabilities of the # models for each value of phi. # deviance= a nmodel-vector containing the deviances for the models. # chi2=a nmodel-vector containing the (DV0-DV1)/scale for the models # npar=a nmodel-vector containing the number of parameters estimated for # each model. # scale=the estimated or assigned scale used (see "glim") # posterior: List containing the Bayesian model mixing results, namely # prob0= a ncol(x) x nphi matrix whose [k.j] element is the posterior # probability that the parameter corresponding to the k-th column of x # is zero, for the j-th value of phi. # mean= a ncol(x) x nphi matrix whose [k,j] element is the posterior # mean of the parameter corresponding to the k-th column of x, for the # j-th value of phi. # sd= as mean, but for the posterior standard deviation. # NOTE: Both mean and sd are CONDITIONAL on the parameter being non-zero. # Unlike the components of glim.est, posterior.bymodel and prior, # posterior$mean and posterior$sd do NOT include the intercept. # glim.est: List containing the GLIM estimates for the different models: # coef: An nmodel-list, each of whose elements is the coef value from # "glim" for one of the models. # se: as coef, but contains standard errors. # var: as coef, but contains variance matrices of the estimates. # posterior.bymodel: List containing model-specific posterior means and sds: # mean= a nmodel-list, each of whose elements is a npar[i]xnphi matrix, # containing the posterior means of the npar[i] parameters of model i, # for each value of phi. # sd= as mean, but for posterior standard deviations. # var= a nmodel-list, each of whose elements is a npar[i] x npar[i] x nphi # array, containing the posterior variance matrix of the parameters of # model i for each value of phi. # prior: List containing the prior distributions. # mean= prior mean for the biggest model (this doesn't depend on phi) # var= similar to corresponding member of posterior.bymodel. # # References: # ---------- # Raftery, A.E. (1988). Approximate Bayes factors for generalized linear # models. Technical Report no. 121, Department of Statistics, University # of Washington. # Raftery, Adrian E. (1995). Bayesian model selection in social research # (with Discussion). Sociological Methodology 1995 (Peter V. Marsden, ed.), # pp. 111-196, Cambridge, Mass.: Blackwells. # An earlier version of this was issued as Working Paper 94-12, Center for Studies in # Demography and Ecology, University of Washington, and this is available as a # Postscript file at http://www.stat.washington.edu/tech.reports/bic.ps. # It is also available via ftp using the following commands: # ftp ftp.stat.washington.edu (or 128.95.17.34) # login as anonymous # enter your email address as your password # ftp>cd /pub/tech.reports # ftp>get bic.ps # ftp>bye # Raftery, A.E. (1996). Approximate Bayes factors and accounting for # model uncertainty in generalized linear models. Biometrika (83: 251-266). # An earlier version of this was issued as Technical Report no. 255, Department of # Statistics, University of Washington, and is available as a Postscript file at # http://www.stat.washington.edu/tech.reports/tr255.ps. # It is also available via ftp in the same way as the previous paper, replacing # "bic.ps" by "tr255.ps". # # BUGS: # ---- # # (1) When error="binomial" and scale="pearson", the program stops saying # "missing argument n in var(yhat)". This appears to be a bug in # the S-PLUS function "glim" itself and not in this function "glib". # This does not appear to happen when scale="deviance", or for other # error distributions. #--------- # Note: In what follows, I try to adopt the following conventions for # dummy indices: i for models (i=1,...,nmodel); j for values of phi # (j=1,...,nphi); k for parameters (k=1,...,npar[i]). prior.spec _ list(mean=F,var=F) if (!missing(priormean)) prior.spec$mean_priormean if (!missing(priorvar)) { prior.spec$var_priorvar phi_1 } #if (missing(priormean)) priormean_rep(0,ncol(x)+1) #if (missing(priorvar)) priorvar_diag(rep(1,ncol(x)+1)) # Do calculations for the null model 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=F, 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) # Initialize quantities for the non-null model calculations 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) # Calculate prior variances 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 { #Current model is not the null model 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 } } # Non-null model calculations for (i in (1:nmodel)) { if (sum(models[i,])==0) { #Current model is null model 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 { #Current model is not the null model 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] } # Construct posterior.bymodel for model i 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])) } } } # Construct BFs and model mixing posterior (lists bf and posterior) for (j in (1:nphi)) { # Calculate posterior model probabilities tmp <- app1[,j] - max(app1[,j]) tmp <- pmw*exp(0.5*tmp) postprob[,j] <- tmp/sum(tmp) # Calculate the posterior probability that each parameter is zero for (k in (1:ncol(x))) prob0[k,j] <- sum(postprob[models[,k]==0,j]) # Calculate the posterior mean and sd (by model mixing) for each parameter 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 #This is the position of parameter k # in the parameter vector of model i 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) } } } # Construct lists for output 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==T) { if (glimvar==F) glim.est <- list(coef=glim.coef,se=glim.se) if (glimvar==T) glim.est <- list(coef=glim.coef,se=glim.se,var=glim.var) } if (glimest==F) glim.est <- NULL if (post.bymodel==T) { if (output.postvar==F) posterior.bymodel <- list(mean=postbym.mean,sd=postbym.sd) if (output.postvar==T) posterior.bymodel <- list(mean=postbym.mean,sd=postbym.sd, var=postbym.var) } if (post.bymodel==F) posterior.bymodel <- NULL if (output.priorvar==F) prior <- list(mean=pmean1) if (output.priorvar==T) prior <- list(mean=pmean1,var=prior.var) # Output list (inputs=inputs,bf=bf,posterior=posterior,glim.est=glim.est, posterior.bymodel=posterior.bymodel,prior=prior) } #======================================================================== glim.pmean.A _ function (x,y,n,error,link,scale,nu,prior.spec) { # Oct 21, 1992 # Calculate prior mean and the matrix A used in the prior, # as given on (C5624.19), # using the reference model defined by x. For now, this is defined # only for the Poisson and Binomial variance functions, and for the # log, logit and cloglog link functions (note that the cloglog link is # referred to as "loglog" in the S function "glim"). # Inputs: as required by the function "glim". plus # nu = prior mean of beta0 in the standardized normal case # Value: # pmean = prior mean vector reflecting nu=0 # A = matrix A for the prior, as on (C5624.19) 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 #Linear predictor # Calculate mean function and d(eta)/d(mu) if (link=="identity") { mu <- eta deta <- 1 } if (link=="log") { mu _ exp(eta) deta _ 1/mu } #deta = d(eta)/d(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)) } #Calculate variance function if (error=="gaussian") v <- 1 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 # Calculate the A matrix and the prior mean w_1/(deta^2 * v) #weights in the GLIM fit 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(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) { # Oct 21, 1992 # Calculates the prior variance given the A matrix calculated by # glim.pmean.A and the prior parameters phi,psi # Inputs: # model is a vector indicating which of the independent variables # in the full x matrix are included in the current model # phi is the prior sd of beta_j in the standardized normal model # In this function, this is a scalar # psi is the prior sd of beta_0 in the standardized normal model # Aot is the output from glim.pmean.A # Value: # phi and psi are echoed from the input # pvar is the prior variance matrix 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) { # Oct 21, 1992 # Calculate E1 for the 4 approximations. The approximate BF is # defined as chi2 + (E1-E0), where E0 refers to the null model with # only an intercept (see the function "E0"). # Inputs: # model is a vector indicating which of the independent variables # in the full x matrix are included in the current model # This is a vector containing the indices of the variables in the # present model. Thus, for example, if the present model includes # variables 1,5,6 in the main matrix x (which does not include the # column of ones), then model=c(1,5,6). # phi is the prior sd of beta_j in the standardized normal model # In this function, this is a vector # psi is the prior sd of beta_0 in the standardized normal model # glimot1 is the output from the relevant glim run # Aot is the output from glim.pmean.A # N is the sample size, needed for the BIC approximation # Value: # the inputs are echoed # app1 is a vector containing the value of E1 for each value of phi # Initialize quantities needed for all the calculations 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)) # Calculate the E1 values 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) { # Oct 21, 1992 # Calculate E0, where the approximate 2 log B10 for M1 against M0 is # chi2 + (E1-E0). # Inputs: # psi is the prior sd of beta_0 in the standardized normal model # glimot0 is the output from the relevant glim run # Aot is the output from glim.pmean.A # Value: # the input psi is echoed # app1 is a vector containing the value of E0 for the approximation 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) } #============================================================ #EXAMPLES: #Example 1: Output from glib function for Finney data #---------------------------------------------------- #> finney.glib <- glib (x,y,n,error="binomial",link="logit",models=finney.models, # glimvar=T,output.priorvar=T,output.postvar=T) #> finney.glib #$inputs: #$inputs$x: # [,1] [,2] # [1,] 3.70 0.825 # [2,] 3.50 1.090 # [3,] 1.25 2.500 # [4,] 0.75 1.500 # [5,] 0.80 3.200 # [6,] 0.70 3.500 # [7,] 0.60 0.750 # [8,] 1.10 1.700 # [9,] 0.90 0.750 #[10,] 0.90 0.450 #[11,] 0.80 0.570 #[12,] 0.55 2.750 #[13,] 0.60 3.000 #[14,] 1.40 2.330 #[15,] 0.75 3.750 #[16,] 2.30 1.640 #[17,] 3.20 1.600 #[18,] 0.85 1.415 #[19,] 1.70 1.060 #[20,] 1.80 1.800 #[21,] 0.40 2.000 #[22,] 0.95 1.360 #[23,] 1.35 1.350 #[24,] 1.50 1.360 #[25,] 1.60 1.780 #[26,] 0.60 1.500 #[27,] 1.80 1.500 #[28,] 0.95 1.900 #[29,] 1.90 0.950 #[30,] 1.60 0.400 #[31,] 2.70 0.750 #[32,] 2.35 0.030 #[33,] 1.10 1.830 #[34,] 1.10 2.200 #[35,] 1.20 2.000 #[36,] 0.80 3.330 #[37,] 0.95 1.900 #[38,] 0.75 1.900 #[39,] 1.30 1.625 #$inputs$y: # [1] 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 0 1 0 0 0 0 1 0 1 0 1 0 1 0 0 1 1 1 0 0 #[39] 1 #$inputs$n: # [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 #[39] 1 #$inputs$error: #[1] "binomial" #$inputs$link: #[1] "logit" #$inputs$models: # [,1] [,2] #[1,] 1 0 #[2,] 0 1 #[3,] 1 1 #$inputs$phi: #[1] 1.00 1.65 5.00 #$inputs$psi: #[1] 1 #$inputs$nu: #[1] 0 #$bf: #$bf$twologB10: # [,1] [,2] [,3] #[1,] 2.48399941 1.5156518 -0.6843579 #[2,] -0.01912596 -0.9868003 -3.1864267 #[3,] 17.26706590 15.7700444 11.6178158 #$bf$postprob: # [,1] [,2] [,3] #[1,] 0.0006159617 0.0008021391 0.0021253384 #[2,] 0.0001762004 0.0002295351 0.0006082901 #[3,] 0.9992078378 0.9989683258 0.9972663715 #$bf$deviance: #[1] 46.98938 49.65482 29.77230 #$bf$df: #[1] 37 37 36 #$bf$chi2: #[1] 7.050458 4.385013 24.267532 #$bf$npar: #[1] 2 2 3 #$bf$scale: #[1] 1 #$posterior: #$posterior$prob0: # [,1] [,2] [,3] #[1,] 0.0001762004 0.0002295351 0.0006082901 #[2,] 0.0006159617 0.0008021391 0.0021253384 #$posterior$mean: # [,1] [,2] [,3] #[1,] 3.557023 3.754823 3.862704 #[2,] 2.439686 2.567756 2.638947 #$posterior$sd: # [,1] [,2] [,3] #[1,] 1.3739701 1.4071282 1.4274537 #[2,] 0.8791384 0.9001849 0.9122723 #$glim.est: #$glim.est$coef: #$glim.est$coef[[1]]: # Intercept X1 # -1.662583 1.335594 #$glim.est$coef[[2]]: # Intercept X1 # -1.353489 0.8439545 #$glim.est$coef[[3]]: # Intercept X1 X2 # -9.529558 3.882138 2.649111 #$glim.est$se: #$glim.est$se[[1]]: #[1] 0.8092169 0.6124761 #$glim.est$se[[2]]: #[1] 0.7921867 0.4389710 #$glim.est$se[[3]]: #[1] 3.2275840 1.4261653 0.9128501 #$glim.est$var: #$glim.est$var[[1]]: # [,1] [,2] #[1,] 0.6548320 -0.4469873 #[2,] -0.4469873 0.3751269 #$glim.est$var[[2]]: # [,1] [,2] #[1,] 0.6275597 -0.3143267 #[2,] -0.3143267 0.1926955 #$glim.est$var[[3]]: # [,1] [,2] [,3] #[1,] 10.417298 -4.3093613 -2.7199119 #[2,] -4.309361 2.0339475 0.9950569 #[3,] -2.719912 0.9950569 0.8332953 #$posterior.bymodel: #$posterior.bymodel$mean: #$posterior.bymodel$mean[[1]]: # [,1] [,2] [,3] #[1,] -1.647979 -1.657155 -1.661944 #[2,] 1.323403 1.331098 1.335115 #$posterior.bymodel$mean[[2]]: # [,1] [,2] [,3] #[1,] -1.3385830 -1.3485054 -1.3536955 #[2,] 0.8340087 0.8400859 0.8432647 #$posterior.bymodel$mean[[3]]: # [,1] [,2] [,3] #[1,] -8.754038 -9.229292 -9.495917 #[2,] 3.558400 3.756769 3.868091 #[3,] 2.439970 2.568153 2.640042 #$posterior.bymodel$sd: #$posterior.bymodel$sd[[1]]: # [,1] [,2] [,3] #[1,] 0.8056591 0.8075636 0.8085558 #[2,] 0.6096254 0.6113954 0.6123171 #$posterior.bymodel$sd[[2]]: # [,1] [,2] [,3] #[1,] 0.7881395 0.7904827 0.7917056 #[2,] 0.4363792 0.4379662 0.4387941 #$posterior.bymodel$sd[[3]]: # [,1] [,2] [,3] #[1,] 3.0960167 3.1768838 3.2213618 #[2,] 1.3731901 1.4059085 1.4239112 #[3,] 0.8789382 0.8998828 0.9114071 #$posterior.bymodel$var: #$posterior.bymodel$var[[1]]: #, , 1 # [,1] [,2] #[1,] 0.6490865 -0.4431225 #[2,] -0.4431225 0.3716431 #, , 2 # [,1] [,2] #[1,] 0.6521590 -0.4456993 #[2,] -0.4456993 0.3738043 #, , 3 # [,1] [,2] #[1,] 0.6537626 -0.4470442 #[2,] -0.4470442 0.3749322 #$posterior.bymodel$var[[2]]: #, , 1 # [,1] [,2] #[1,] 0.6211639 -0.3109171 #[2,] -0.3109171 0.1904268 #, , 2 # [,1] [,2] #[1,] 0.6248630 -0.3131826 #[2,] -0.3131826 0.1918144 #, , 3 # [,1] [,2] #[1,] 0.6267978 -0.3143677 #[2,] -0.3143677 0.1925402 #$posterior.bymodel$var[[3]]: #, , 1 # [,1] [,2] [,3] #[1,] 9.585319 -3.9619743 -2.4979929 #[2,] -3.961974 1.8856510 0.9037059 #[3,] -2.497993 0.9037059 0.7725324 #, , 2 # [,1] [,2] [,3] #[1,] 10.092591 -4.1748682 -2.6339864 #[2,] -4.174868 1.9765787 0.9596574 #[3,] -2.633986 0.9596574 0.8097890 #, , 3 # [,1] [,2] [,3] #[1,] 10.377172 -4.2942946 -2.7102850 #[2,] -4.294295 2.0275230 0.9910891 #[3,] -2.710285 0.9910891 0.8306629 #$prior: #$prior$mean: #[1] 0.01940434 0.00000000 0.00000000 #$prior$var: #$prior$var[[1]]: #, , 1 # [,1] [,2] #[1,] 74.02733 -51.04583 #[2,] -51.04583 40.66962 #, , 2 # [,1] [,2] #[1,] 184.3868 -138.9723 #[2,] -138.9723 110.7231 #, , 3 # [,1] [,2] #[1,] 1611.692 -1276.146 #[2,] -1276.146 1016.741 #$prior$var[[2]]: #, , 1 # [,1] [,2] #[1,] 61.85684 -29.40020 #[2,] -29.40020 16.65493 #, , 2 # [,1] [,2] #[1,] 151.25263 -80.04205 #[2,] -80.04205 45.34304 #, , 3 # [,1] [,2] #[1,] 1307.430 -735.0050 #[2,] -735.005 416.3731 #$prior$var[[3]]: #, , 1 # [,1] [,2] [,3] #[1,] 125.92620 -51.04583 -29.40020 #[2,] -51.04583 40.66962 0.00000 #[3,] -29.40020 0.00000 16.65493 #, , 2 # [,1] [,2] [,3] #[1,] 325.68148 -138.9723 -80.04205 #[2,] -138.97228 110.7231 0.00000 #[3,] -80.04205 0.0000 45.34304 #, , 3 # [,1] [,2] [,3] #[1,] 2909.164 -1276.146 -735.0050 #[2,] -1276.146 1016.741 0.0000 #[3,] -735.005 0.000 416.3731 #=================================================== #Example 2: Test run of glib with Yates (teeth) data. #---------------------------------------------------- #> glib.yates <- glib (x,y,n,models=models,glimvar=T,output.priorvar=T,output.postvar=T) # # Error, link, scale, phi, psi, nu and pmw are all defaults #> glib.yates #$inputs: #$inputs$x: # [,1] [,2] [,3] #[1,] 0 0 0 #[2,] 0 1 0 #[3,] 1 0 0 #[4,] 1 1 1 #$inputs$y: #[1] 4 16 1 21 #$inputs$n: #[1] 1 1 1 1 #$inputs$error: #[1] "poisson" #$inputs$link: #[1] "log" #$inputs$models: # [,1] [,2] [,3] #[1,] 1 1 0 #[2,] 1 1 1 #$inputs$phi: #[1] 1.00 1.65 5.00 #$inputs$psi: #[1] 1 #$inputs$nu: #[1] 0 #$bf: #$bf$twologB10: # [,1] [,2] [,3] #[1,] 20.59868 19.28636 15.22839 #[2,] 21.30901 19.72433 14.26971 #$bf$postprob: # [,1] [,2] [,3] #[1,] 0.4121309 0.4454715 0.6175917 #[2,] 0.5878691 0.5545285 0.3824083 #$bf$deviance: #[1] 2.509921e+00 -8.896602e-16 #$bf$df: #[1] 1 0 #$bf$chi2: #[1] 27.65769 30.16761 #$bf$npar: #[1] 3 4 #$bf$scale: #[1] 1 #$posterior: #$posterior$prob0: # [,1] [,2] [,3] #[1,] 0.0000000 0.0000000 0.0000000 #[2,] 0.0000000 0.0000000 0.0000000 #[3,] 0.4121309 0.4454715 0.6175917 #$posterior$mean: # [,1] [,2] [,3] #[1,] -0.2592935 -0.433626 -0.4384975 #[2,] 1.6539873 1.688133 1.7693294 #[3,] 0.7488731 1.109594 1.5687607 #$posterior$sd: # [,1] [,2] [,3] #[1,] 0.6434585 0.8400730 0.9840144 #[2,] 0.5262522 0.5601034 0.5839693 #[3,] 0.7306174 0.9284488 1.1308835 #$glim.est: #$glim.est$coef: #$glim.est$coef[[1]]: # Intercept X1 X2 # 0.8675006 0.09531018 2.00148 #$glim.est$coef[[2]]: # Intercept X1 X2 X3 # 1.386294 -1.386294 1.386294 1.658228 #$glim.est$se: #$glim.est$se[[1]]: #[1] 0.4755254 0.3089512 0.4764052 #$glim.est$se[[2]]: #[1] 0.4999992 1.1180116 0.5590163 1.1662199 #$glim.est$var: #$glim.est$var[[1]]: # [,1] [,2] [,3] #[1,] 0.22612445 -4.999816e-02 -1.999349e-01 #[2,] -0.04999816 9.545084e-02 4.670940e-17 #[3,] -0.19993488 4.670940e-17 2.269619e-01 #$glim.est$var[[2]]: # [,1] [,2] [,3] [,4] #[1,] 0.2499992 -0.2499992 -0.2499992 0.2499992 #[2,] -0.2499992 1.2499499 0.2499992 -1.2499499 #[3,] -0.2499992 0.2499992 0.3124992 -0.3124992 #[4,] 0.2499992 -1.2499499 -0.3124992 1.3600689 #$posterior.bymodel: #$posterior.bymodel$mean: #$posterior.bymodel$mean[[1]]: # [,1] [,2] [,3] #[1,] 0.97121027 0.90804149 0.87349618 #[2,] 0.08994549 0.09326691 0.09508333 #[3,] 1.88885031 1.95858306 1.99671801 #$posterior.bymodel$mean[[2]]: # [,1] [,2] [,3] #[1,] 1.2881117 1.3087990 1.372029 #[2,] -0.5041305 -0.8568968 -1.300234 #[3,] 1.4893345 1.4708715 1.402096 #[4,] 0.7488731 1.1095940 1.568761 #$posterior.bymodel$sd: #$posterior.bymodel$sd[[1]]: # [,1] [,2] [,3] #[1,] 0.4619521 0.4695197 0.4736071 #[2,] 0.3001304 0.3056216 0.3085833 #[3,] 0.4628066 0.4712722 0.4758381 #$posterior.bymodel$sd[[2]]: # [,1] [,2] [,3] #[1,] 0.4550705 0.4749823 0.4950572 #[2,] 0.7040615 0.8918755 1.0843833 #[3,] 0.5051261 0.5304290 0.5547354 #[4,] 0.7306174 0.9284488 1.1308835 #$posterior.bymodel$var: #$posterior.bymodel$var[[1]]: #, , 1 # [,1] [,2] [,3] #[1,] 0.21339971 -4.718393e-02 -1.886843e-01 #[2,] -0.04718393 9.007824e-02 1.613246e-12 #[3,] -0.18868427 1.613443e-12 2.141900e-01 #, , 2 # [,1] [,2] [,3] #[1,] 0.22044875 -4.892629e-02 -1.956501e-01 #[2,] -0.04892629 9.340456e-02 1.734678e-12 #[3,] -0.19565013 1.734794e-12 2.220975e-01 #, , 3 # [,1] [,2] [,3] #[1,] 0.22430368 -4.987916e-02 -1.994596e-01 #[2,] -0.04987916 9.522366e-02 1.803016e-12 #[3,] -0.19945958 1.803061e-12 2.264219e-01 #$posterior.bymodel$var[[2]]: #, , 1 # [,1] [,2] [,3] [,4] #[1,] 0.20708913 -0.09516072 -0.20175323 0.08592216 #[2,] -0.09516072 0.49570253 0.07581017 -0.46255709 #[3,] -0.20175323 0.07581017 0.25515236 -0.12546749 #[4,] 0.08592216 -0.46255709 -0.12546749 0.53380184 #, , 2 # [,1] [,2] [,3] [,4] #[1,] 0.2256082 -0.1563093 -0.2234279 0.1511332 #[2,] -0.1563093 0.7954419 0.1442973 -0.7749410 #[3,] -0.2234279 0.1442973 0.2813549 -0.2000331 #[4,] 0.1511332 -0.7749410 -0.2000331 0.8620171 #, , 3 # [,1] [,2] [,3] [,4] #[1,] 0.2450817 -0.2346970 -0.2458941 0.2338902 #[2,] -0.2346970 1.1758871 0.2327084 -1.1724998 #[3,] -0.2458941 0.2327084 0.3077314 -0.2941947 #[4,] 0.2338902 -1.1724998 -0.2941947 1.2788976 #$prior: #$prior$mean: #[1] 2.710514 0.000000 0.000000 0.000000 #$prior$var: #$prior$var[[1]]: #, , 1 # [,1] [,2] [,3] #[1,] 3.7922108 -0.8382782 -3.353113 #[2,] -0.8382782 1.6003492 0.000000 #[3,] -3.3531127 0.0000000 3.806236 #, , 2 # [,1] [,2] [,3] #[1,] 9.636706 -2.282212 -9.128849 #[2,] -2.282212 4.356951 0.000000 #[3,] -9.128849 0.000000 10.362478 #, , 3 # [,1] [,2] [,3] #[1,] 85.2[2,] -20.95695 40.00873 0.00000 #[3,] -83.82782 0.00000 95.15590 #$prior$var[[2]]: #, , 1 # [,1] [,2] [,3] [,4] #[1,] 4.1913909 -0.8382782 -3.353113 -0.7983602 #[2,] -0.8382782 1.6003492 0.000000 0.0000000 #[3,] -3.3531127 0.0000000 3.806236 0.0000000 #[4,] -0.7983602 0.0000000 0.000000 1.5967203 #, , 2 # [,1] [,2] [,3] [,4] #[1,] 10.723474 -2.282212 -9.128849 -2.173536 #[2,] -2.282212 4.356951 0.000000 0.000000 #[3,] -9.128849 0.000000 10.362478 0.000000 #[4,] -2.173536 0.000000 0.000000 4.347071 #, , 3 # [,1] [,2] [,3] [,4] #[1,] 95.20445 -20.95695 -83.82782 -19.95900 #[2,] -20.95695 40.00873 0.00000 0.00000 #[3,] -83.82782 0.00000 95.15590 0.00000 #[4,] -19.95900 0.00000 0.00000 39.91801