Check posterior probabilities of RCPs for each model fit

Source: R/regional_mix_s3-class.R

This function checks for the minimum posterior probabilities of each RCP across site fitted in the model. Basically sometime you can get effectively empty or single site RCPs, this function will help you check for these and remove them from model selection.

check_RCP_posteriors(object, min_membership, ...)

Arguments

object

An RCP fit or a multifit RCP object.

min_membership

The minimum number of sites (by highest posterior membership) an RCP needs to have to not be flagged as effectively empty.

...

other arguments

Value

For a regional_mix object, a table of the number of sites hard-assigned (by maximum posterior probability) to each RCP, with a warning if any RCP has fewer than min_membership sites. For a regional_mix.multifit object, a list of such tables, one per random start.

Examples

if (FALSE) {
set.seed( 151)
n <- 100
S <- 10
nRCP <- 3
my.dist <- "negative.binomial"
X <- as.data.frame( cbind( x1=runif( n, min=-10, max=10),
                          x2=runif( n, min=-10, max=10)))
Offy <- log( runif( n, min=30, max=60))
pols <- list()
pols[[1]] <- poly( X$x1, degree=3)
pols[[2]] <- poly( X$x2, degree=3)
X <- as.matrix( cbind( 1, X, pols[[1]], pols[[2]]))
colnames( X) <- c("const", 'x1', 'x2', paste( "x1",1:3,sep='.'),
paste( "x2",1:3,sep='.'))
p.x <- ncol( X[,-(2:3)])
p.w <- 3
W <- matrix(sample( c(0,1), size=(n*p.w), replace=TRUE), nrow=n, ncol=p.w)
colnames( W) <- paste( "w",1:3,sep=".")
alpha <- rnorm( S)
tau.var <- 0.5
b <- sqrt( tau.var/2)
tau <- matrix( rexp( n=(nRCP-1)*S,rate=1/b) - rexp( n=(nRCP-1)*S, rate=1/b),
 nrow=nRCP-1, ncol=S)
beta <- 0.2 * matrix( c(-1.2, -2.6, 0.2, -23.4, -16.7, -18.7, -59.2,
 -76.0,-14.2, -28.3, -36.8, -17.8, -92.9,-2.7), nrow=nRCP-1, ncol=p.x)
gamma <- matrix( rnorm( S*p.w), ncol=p.w, nrow=S)
logDisp <- log( rexp( S, 1))
set.seed(121)
simDat <- regional_mix.simulate( nRCP=nRCP, S=S, p.x=p.x, p.w=p.w, n=n,
alpha=alpha, tau=tau, beta=beta, gamma=gamma, X=X[,-(2:3)], W=W,
family=my.dist, logDisp=logDisp, offset=Offy)
form.RCP <- paste( paste( paste('cbind(', paste( paste( 'spp', 1:S, sep=''),
 collapse=','), sep=''),')',sep=''),'~x1.1+x1.2+x1.3+x2.1+x2.2+x2.3',sep='')
form.spp <- ~w.1+w.2+w.3
fm_regional_mix <- regional_mix(rcp_formula=form.RCP,species_formula=form.spp,
                                data=simDat, family='negative.binomial', nRCP=5)
check_RCP_posteriors(fm_regional_mix, min_membership=5)
}