Многофакторного ТТЕСТ использования г и winbug

Question

Как можно сделать разницу в средствах (ТТЕСТ) для многофакторного с использованием 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) 
} 

Answer 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 Ошибка т р-значение

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

Когда я скопировал результаты от моего экрана, чтобы ТАК было трудно сделать ярлыки результатов должным образом разнесенные, мои извинения.