multivariée TTEST utilisant r et winbug

Question

Comment puis-je faire la différence dans les moyens (TTEST) pour une à plusieurs variables en utilisant R et WinBUGS14 J'ai un résultat à plusieurs variables y et la x variable catégorique. Je suis en mesure d'obtenir les moyennes des valeurs échantillonnées MCMC de la multivariable en utilisant le code ci-dessous, mais comment puis-je tester la différence de moyennes par variable x?

Voici le code 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) 

et c'est le code WinBUGS14 appelé 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 étapes:

1. charge une fonction pour exécuter le t.test en utilisant des statistiques d'échantillon au lieu de le faire directement.

   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 la moyenne et l'écart-type des choses que nous voulons tester contre x, puis les passer à la fonction.

   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) 
   

différence de moyens Std erreur t p-valeur

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

Lorsque j'ai copié les résultats de mon écran de sorte qu'il était difficile de faire les étiquettes des résultats correctement espacées, mes excuses.