4
$\begingroup$

A traditional t-test to compare two means assumes independent random variables, $X_1,...,X_n \sim n(\mu_1,\sigma_1 ^2)$ and $Y_1, ...,Y_m \sim n(\mu_2,\sigma_2 ^2)$. One usually assumes normal distributions with $\sigma_1$ = $\sigma_2$. Then the statistics $$\frac{\bar{X}-\bar{Y}}{S\sqrt{1/n+1/m}}$$ is used and it's observed value is compared to $t_{n+m-2,1-\alpha/2}$ where $S^2$ is the pooled estimate of variance, valid under both the null and alternative, given the assumption of equal variances.

If the variances are known one can use $$\frac{\bar{X}-\bar{Y}}{\sqrt{\sigma_1^2/n+\sigma_2^2/m}}$$ and compare with a $z_{1-\alpha/2}$ test.

If we assume $n=m$ and for example $\sigma_1^2=1$ and $\sigma_2^2=2$ and investigate using simulation which test is better. How would we do that. It seems wise to check level and powers but how would one do that?

$\endgroup$
  • 1
    $\begingroup$ Are the $\sigma$'s known to the person doing the test or not? $\endgroup$ – Glen_b Nov 16 '13 at 17:26
  • $\begingroup$ Are you asking how to conduct a simulation study to investigate this? $\endgroup$ – gung Nov 16 '13 at 17:29
  • $\begingroup$ Well I'm asking for how to check the levels and powers of the test. And yes $\sigma's$ are known, $\sigma_1^2=1$ and $\sigma_2^2=2$. $\endgroup$ – Raxel Nov 16 '13 at 17:39
  • $\begingroup$ Then why would the person doing the test use a t-test - where you estimate something that's already known - and even then why would you use one that assumes what you know with certainty to be false? Wouldn't you at least do something like a Welch-Satterthwaite t-test? $\endgroup$ – Glen_b Nov 16 '13 at 19:17
  • $\begingroup$ What language would you want to use for the simulations, R? $\endgroup$ – gung Nov 16 '13 at 19:20
3
$\begingroup$

Without loss of generality, you can assume that the mean of the first group, $\mu_1$ is zero, and the mean of the second group is $\mu_2 = \mu_1+\delta = \delta$.

You'll need to assume some particular sample sizes (specific m and n), or do this over a grid of such values.

Let's take level: here $\delta = 0$.

So you know the parameters of both distributions. Basically:

Repeat nsim times:
    simulate sample 1
    simulate sample 2
    compute z statistic 
    if z rejects, increment the count of z rejections
    compute t statistic
    if t rejects, increment the count of t rejections

Compute rejection rates for t and z

The case for power is the same, expect you do it for some set of values of $\delta$.

A common thing to do at the end is plot the rejection rates against delta, so you can get an idea what the power curve looks like.


Edit: Okay, you work in R; that would have been handy to know.

In R you probably want to use replicate rather than a loop.

Since you would do this for a set of delta values, you'd probably wrap that in a call to sapply.

[and then, if you then did that over various values of n, you might even put that in a loop]


Here's an example of some R code for estimating level and power at a fixed n

nsim <- 1000
sigxsq <- 1
sigysq <- 2
n <- 10
m <- 10
sigx <- sqrt(sigxsq)
sigy <- sqrt(sigysq)

testasample <- function(n,m,delta=0,mux=0,muy=mux+delta,sigx=1,sigy=1,alpha=0.10,
                 zs=sqrt(sigx^2/n+sigy^2/m)){
  x <- rnorm(n,mux,sigx)
  y <- rnorm(m,muy,sigy)
  tt <- t.test(x,y,var.equal=TRUE)
  tp <- tt$p.value
  zp <- as.numeric(2*pnorm(-abs(diff(tt$estimate)/zs)))
  c(tp,zp) < alpha
}

del <- c(0,0.2,0.5,0.7,1,1.5,2)
res <- sapply(del,function(d){
         rowMeans(replicate(nsim,testasample(n,m,delta=d,sigx=sigx,sigy=sigy)))})

res <- cbind(del,t(res))
colnames(res)<-c("delta","trej","zrej")

It does a little unnecessary calculation, but that will hardly change the overall speed so I didn't bother with it.

Here's a plot that I generated after running the above code:

 plot(trej~delta,res,pch="t",col=4,ylim=c(0,1),ylab="rejection rate")
 points(zrej~delta,res,pch="z",col=2)
 abline(h=c(0,1,.1),col=8,lty=c(1,1,3))

power curves

Note that you can compute standard errors (and intervals) for these estimated rejection rates.

$\endgroup$
  • $\begingroup$ Since $m=n$ I won't need to do this over a grid just over a simple for loop. So after computing the rejection rates I would prefer the one with lower rejection rate? And then I would run another loop with extra for loop outside where I test different appropriate values of $\delta$ and look at the rejection rate again and plot it. I'm basically thinking how I would interpret the results? $\endgroup$ – Raxel Nov 16 '13 at 19:51
  • $\begingroup$ For level you want one with its rejection rate as close as possible to the specified significance level (or possibly 'closest without going over', depending on what you're trying to achieve), while for power you want the highest possible rejection rate, given the same actual (rather than nominal) significance level. You should change your question to incorporate the new information. Before I go into much detail, I ask again -- is this for some subject? $\endgroup$ – Glen_b Nov 16 '13 at 20:15
  • $\begingroup$ This is just selfstudy. I picked $\alpha =0.1$ and I want to compare these values. I will post my code for the level very soon $\endgroup$ – Raxel Nov 16 '13 at 20:20
  • $\begingroup$ Thanks, I think my code below is right for the first part. $\endgroup$ – Raxel Nov 16 '13 at 20:42
  • $\begingroup$ Hard to deduce from all this which test is better overall although I expected z-test to become increasingly better when the m and n became bigger. Any thoughts on that? $\endgroup$ – Raxel Nov 16 '13 at 20:53
0
$\begingroup$
n <- 1000
h <- 30
rt <- rep(0,h)
rz <- rep(0,h)
for (j in 2:h) {
    for (i in 1:n) {
        p <- j
        tstat <- 0
        zstat <- 0
        Xvec <- rnorm(p, mean = 0, sd = 1)
        Yvec <- rnorm(p,mean = 0, sd = sqrt(2))
        xsum <- 0
        ysum <- 0
    for (k in 1:p) { 
        xsum <- (Xvec[k]-mean(Xvec))^2+xsum 
        ysum <- (Yvec[k]-mean(Yvec))^2+ysum 
    }
        Ssq <- (xsum+ysum)/(2*p-2)
        tstat <- (mean(Xvec)-mean(Yvec))/(sqrt(Ssq)*sqrt(2/p))
        if (abs(tstat) > qt(0.95,2*p-2)) {
            rt[p] <- rt[p]+1
        }
        zstat <- (mean(Xvec)-mean(Yvec))/(sqrt(3/p))
        if (abs(zstat) > qnorm(0.95, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)) { 
            rz[p] <- rz[p]+1
        }
    }
}
rt <- rt/n
rz <- rz/n
$\endgroup$
  • $\begingroup$ Answers without a single word of explanation, a bare code doing what only God knows, are hardly acceptable. $\endgroup$ – ttnphns Nov 17 '13 at 6:19
  • $\begingroup$ Well it was just mine attempt to create a code for the level problem. $\endgroup$ – Raxel Nov 17 '13 at 15:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.