使用r和winbug的多变量测试

Question

如何使用R和WinBUGS14对多变量进行均值差(t检验)我有一个多变量结果y和分类变量x。我可以使用下面的代码从多元变量中获得MCMC采样值的平均值，但是如何通过变量x来测试平均值的差异

下面是R代码

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)

这是名为mult_normal.bug的WinBUGS14代码

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)
}

Answer 1

2个步骤：

1. 使用样本统计数据加载一个函数来运行t.test，而不是直接执行。

t.test2 <-函数(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 <- ( (s1^2/n1 + s2^2/n2)^2 )/( (s1^2/n1)^2/(n1-1) + (s2^2/n2)^2/(n2-1) )} else {#合并标准差，按样本大小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(“均值差”，“标准误差”，"t"，"p-value") return(dat) }

• 解析出我们想要对x进行测试的东西的平均值和标准差，然后将它们传递给函数。

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

平均标准差t p值的

差

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

当我将结果从我的屏幕复制到，所以很难让结果的标签适当地间隔开，我道歉。