CODEWITHSUNDEEP
rbayesianwinbugsmultivariate-testing
Keniajin
Keniajin

Reputation: 1659

Multivariate ttest using r and winbug

How can I do difference in means (ttest) for a multivariate using R and WinBUGS14 I have a multivariate outcome y and the categorical variable x. I am able to get the means of the MCMC sampled values from the multivariate using the code below, but how can I test for the difference in means by variable x?

Here is the R code

library(R2WinBUGS)
library(MASS)  # need to mvrnorm
library(MCMCpack) # need for rwish
# Generate synthetic data
N <- 500
#we use this to simulate the data
S <- matrix(c(1,.2,.2,5),nrow=2)

#Produces one or more samples from the specified multivariate normal distribution.
#produces 2 variables with the given distribution
y <- mvrnorm(n=N,mu=c(1,3),Sigma=S)
x <- rbinom(500, 1, 0.5)

# Set up for WinBUGS
#set up of the mu0 values
mu0 <-  as.vector(c(0,0))


#covariance matrices
# the precisions
S2 <- matrix(c(1,0,0,1),nrow=2)/1000 #precision for unkown mu
# precison matrix to be passes to the wishart distribution for the tau
S3 <-  matrix(c(1,0,0,1),nrow=2)/10000 

#the data for the winbug code
data <- list("y","N","S2","S3","mu0")
inits <- function(){
                    list( mu=mvrnorm(1,mu0,matrix(c(10,0,0,10),nrow=2) ), 
                      tau <-  rwish(3,matrix(c(.02,0,0,.04),nrow=2)) )
                  }

# Run WinBUGS
bug_file <- paste0(getwd(), "/codes/mult_normal.bug")
multi_norm.sim <- bugs(data,inits,model.file=bug_file,
parameters=c("mu","tau"),n.chains = 2,n.iter=4010,n.burnin=10,n.thin=1,
bugs.directory="../WinBUGS14/",codaPkg=F)

print(multi_norm.sim,digits=3)

and this is the WinBUGS14 code called mult_normal.bug

model{
for(i in 1:N)
{
y[i,1:2] ~ dmnorm(mu[],tau[,])
}
mu[1:2] ~ dmnorm(mu0[],S2[,])
#parameters of a wishart
tau[1:2,1:2] ~  dwish(S3[,],3)
}

Upvotes: 1

Views: 354

Answers (1)

Hack-R
Hack-R

Reputation: 23210

2 Steps:

  1. Load a function to run the t.test using sample statistics instead of doing it directly.

    t.test2 <- function(m1,m2,s1,s2,n1,n2,m0=0,equal.variance=FALSE)
{
   if( equal.variance==FALSE ) 
   {
     se <- sqrt( (s1^2/n1) + (s2^2/n2) )
     # welch-satterthwaite df
     df <- ( (s1^2/n1 + s2^2/n2)^2 )/( (s1^2/n1)^2/(n1-1) + (s2^2/n2)^2/(n2-1) )
   } else
   {
     # pooled standard deviation, scaled by the sample sizes
     se <- sqrt( (1/n1 + 1/n2) * ((n1-1)*s1^2 + (n2-1)*s2^2)/(n1+n2-2) ) 
     df <- n1+n2-2
   }      
   t <- (m1-m2-m0)/se 
   dat <- c(m1-m2, se, t, 2*pt(-abs(t),df))    
   names(dat) <- c("Difference of means", "Std Error", "t",      "p-value")
  return(dat) 
 }

  2. Parse out the mean and standard deviation of the things we want to test against x, then pass them to the function.

    mu1   <- as.data.frame(multi_norm.sim$mean)$mu[1]
sdmu1 <- multi_norm.sim$sd$mu[1]
t.test2( mean(x), as.numeric(mu1), s1 = sd(x), s2 = sdmu1, 500, 500)

Difference of means Std Error t p-value 

  -4.950656e-01        2.246905e-02       -2.203323e+01        5.862968e-76

When I copied the results from my screen to SO it was hard to make the labels of the results properly spaced apart, my apologies.

Upvotes: 1

Related Questions