direct call of function rmultireg for testing Linux installation

test_rmultireg.r

require(plyr)
#require(siar)

t1 <- data.frame(SizeClass=c("Meiofauna","Smaller Macrofauna","Bigger Macrofauna"),
           num=rpois(3,10),
           md13C=rnorm(3,-10,10),
           sd13C=rlnorm(3,1,1),
           md15N=rnorm(3,5,5) ,
           sd15N=rlnorm(3,1,1))
t1

hh <- function(x) {
  with(x,data.frame(SizeClass=SizeClass,
             #             d13C=rnorm(num,md13C,sd13C),d15N=rlnorm(num,log(md15N),log(sd15N)))
             d13C=rnorm(num,md13C,sd13C),d15N=rnorm(num,md15N,sd15N)))
}
t2 <- ddply(t1, .(SizeClass), hh)

# ---------------------------------------------------------
# test rmultireg directly, by avoiding call to bayesMVN() which
# is in turn called by siber.hull.metrics()

with(t2, { 
     reps <- 10^3 # fewer samples for testing
     
     # split the data based on SizeClass
     M <- length(unique(SizeClass))
     spx <- split(d13C, SizeClass)
     spy <- split(d15N, SizeClass)

     # just test the first group for now
     x <- spx[[1]]
     y <- spy[[1]]
     
     # prepare information for call to rmultireg
     Y <- cbind(x, y) # response variable
     X <- matrix(1, length(x), 1) # intercept
     Bbar <- matrix(0, 1, 2) # prior means
     A <- matrix(10^-3, 1, 1) # prior precision
     nu <- 2 # d.f. for prior on Sigma
     V <- 2 * diag(2) # location for prior on Sigma
     
     # matrices to store posterior draws
     b <- matrix(double(reps * 2), ncol = 2)
     S <- matrix(double(reps * 4), ncol = 4)
     
     # loop rmultireg
     for (i in 1:reps) {
       out <- rmultireg(Y, X, Bbar, A, nu, V)
       b[i, ] <- out$B
       S[i, ] <- out$Sigma
     }
     
     
     
     # quick test that the means at least look ok
     par(mfrow=c(1,2))
     hist(b[,1], 100) # posterior of d13C
     abline(v = mean(x), col = "red", lwd = 2) 
     hist(b[,2], 100) # posterior of d15N
     abline(v = mean(y), col = "red", lwd = 2)
     
})