1

我如何使用多元变量的平均值（ttest）进行差异，R并且WinBUGS14 我有一个多元结果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)

这是WinBUGS14代码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)
}
4

1 回答 1

1

2个步骤：

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

    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. 解析出我们要测试的事物的均值和标准差x，然后将它们传递给函数。

    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)

均值差 Std Error t p 值

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

当我将结果从我的屏幕复制到 SO 时，很难使结果的标签正确间隔开，我很抱歉。

于 2016-06-23T13:45:20.723 回答