Solutions for Stat 421 Homework 7: Two-sample t-test and Normality Testing - Prof. Friedri, Assignments of Statistics

Download Solutions for Stat 421 Homework 7: Two-sample t-test and Normality Testing - Prof. Friedri and more Assignments Statistics in PDF only on Docsity! Stat 421, Fall 2008 Fritz Scholz Homework 7 Solutions Problem 1: Write a function Satterthwaite(xdat,ydat){....} that takes an x- sample xdat and a y-sample ydat and computes the two-sample t-statistic txy (we called it t*(X,Y) on slide 107) without assumption of equal variance in both sampled populations, computes the degrees of freedom fhat by the Satterthwaite method, and the corresponding p-value , pval, corresponding to the observed value txy. The output should be a list with the values txy, fhat, and pval. Compare your results with those of t.test(xdat,ydat,var.equal=F) when the two samples are given as xdat=c(-0.7,-1.6,1.1,-0.2,0.4,-1.4,1.6,0.2,-0.3,0.3) ydat=c(-0.4,-2.7,2.6,-1.6,2.8,1.5,-2.2,3.2,-0.1,2.2,-1.9, 0.9) This means that you should not use t.test as a shortcut in your function Satterthwaite. Give your function and the results of your comparison. After having convinced yourself that t.test does produce the correct p-values (according to slide 108 in Stat421NormalPopulation.pdf), write another function test.Satterthwaite(Nsim=10000,sigx=1,sigy=2,m=10,n=20){...} that generates in a loop, Nsim times, independent samples of size m and n from normal populations with means 0 and standard deviations sigx and sigy, respectively, and uses t.test to compute the p-value for each such pair of samples (Use out=t.test(xdat,ydat,var.equal=F) and out$p.value to get the p-value). This function should then also produce a histogram of these Nsim p-values, alternating colors in the bars and using breaks=seq(0,1,.025) and probability=T. During the development of this function use Nsim=1000, but ultimately use the default arguments as shown above. Show the histogram (with appropriate annotations) and give the code. The histogram should look approximately uniform, since p-values (computed from continuous null distributions) should have a U(0,1) distribution when the data are generated according to the null distribution. Does this distribution of p-values suggest that these p-values will lead to a more or less correct type I error probability alpha, if we reject the hypothesis of equal means whenever an observed p-value is less than or equal to alpha, whatever alpha is chosen? Solution: The code Satterthwaite= function (xdat,ydat){ m=length(xdat) n=length(ydat) sx2=var(xdat) sy2=var(ydat) fhat=(sx2/m+sy2/n)^2/((sx2/m)^2/(m-1)+(sy2/n)^2/(n-1)) txy=(mean(ydat)-mean(xdat))/sqrt(sx2/m+sy2/n) pval=2*pt(-abs(txy),fhat) out=list(pval=pval,txy=txy,fhat=fhat) out } The comparison on the given samples yields > Satterthwaite(xdat,ydat) $pval [1] 0.5535037 $txy [1] 0.6050084 $fhat [1] 16.29990  t.test(ydat,xdat,var.equal=F)  Welch Two Sample t-test data: ydat and xdat t = 0.605, df = 16.3, p-value = 0.5535 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -1.045287 1.881954 sample estimates: mean of x mean of y 0.3583333 -0.0600000 The results are correct within given precision. The simulation function is given below, followed by the resulting plot test.Satterthwaite=function (Nsim=10000,sigx=1,sigy=2,m=10,n=20) {pval=rep(0,Nsim) for(i in 1:Nsim){ xdat=rnorm(m,0,sigx) ydat=rnorm(n,0,sigy) out=t.test(xdat,ydat,var.equal=F) pval[i]=out$p.value } hist(pval,xlab="P-value of Satterthwaite-Welch t-test", breaks=seq(0,1,.025), col=c("blue","orange"),probability=T,main="") segments(0,1,1,1,lty=2) } and indicate the proportion of values that are greater or equal to that critical value, indicating the power achieved for uniform samples. Now repeat this by writing a slightly modified function using the Anderson-Darling test of fit statistic (using ad.test instead of lillie.test). Here use breaks=seq(0,1.1*M,.1). Show your properly annotated plots (showing power, n, Nsim, and axis labels, indicating the test used) and interpret what you see in the compared histograms, why are they different? How do the results change for the two types of tests? Does it make sense (see slide 124 of Stat421Normal Population.pdf)? In understanding the different results make a plot that superimposes a uniform(0,1) CDF over a normal CDF with mean ½ and variance 1/12 (variance of U(0,1)). Show this plot as well. Give the code for the two functions. Reminder: For each new R session you have to invoke library(nortest) before you can use lillie.test or ad.test. Solution: The two sets of plots are as follows Sample size n = 100 , Nsim = 10000 KS-Statistic for Normal Samples F re q u e n cy 0.00 0.05 0.10 0.15 0.20 0 4 0 0 8 0 0 1 2 0 0  0.05 KS-Statistic for Uniform Samples F re q u e n cy 0.00 0.05 0.10 0.15 0.20 0 2 0 0 6 0 0 1 0 0 0  0.5795 Sample size n = 100 , Nsim = 10000 AD-Statistic for Normal Samples F re q u e n cy 0 1 2 3 4 0 5 0 0 1 5 0 0 2 5 0 0  0.05 AD-Statistic for Uniform Samples F re q u e n cy 0 1 2 3 4 0 2 0 0 4 0 0 6 0 0 8 0 0  0.9479 The top histogram in each pair shows the null distribution of the KS-discrepancy statistic or Anderson-Darling discrepancy statistic. The bottom histogram in each pair shows the distribution of the test statistic for normality when the samples do not come from a normal distribution. In fact, in that case the sampled distribution is uniform. The distribution of the test statistic then moves to the right, i.e., yielding larger values of the discrepancy metric, inducing us to reject the hypothesis of normality when it is false (desirable!). In the case of the KS-statistic there is still a fair amount of overlap between the null distribution and the distribution under the U(0,1) alternative, even for samples as large as n=100. Hence we still have .4205 estimated probability of type II error here when the significance level is set to .05, while in the case of the Anderson-Darling criterion the probability of type II error is .0521. Of course, these numbers are based on 10000 simulations and may vary to some degree from one simulation to the next. -0.5 0.0 0.5 1.0 1.5 0 .0 0 .2 0 .4 0 .6 0 .8 1 .0 x C D F Normal CDF Uniform CDF The fact that the AD criterion shows clearer separation is due to the fact that the normal and uniform distributions show strongly different tail behavior (see plot above) and discrepancies extend over a wider range of x values. The AD-test integrates discrepancies over the full x range and weights tail discrepancies more strongly (see slide 124). Thus, all makes sense. The code behind the first two plots follows lillie.test.simulation=function (Nsim=10000,n=100,alpha=.05) { sim.out.norm=rep(0,Nsim) sim.out.unif=rep(0,Nsim) for(i in 1:Nsim){ out.n=lillie.test(rnorm(n)) out.u=lillie.test(runif(n)) sim.out.norm[i]=out.n$statistic sim.out.unif[i]=out.u$statistic } m=min(sim.out.norm,sim.out.unif) M=max(sim.out.unif,sim.out.norm) par(mfrow=c(2,1)) out.hist=hist(sim.out.norm,breaks=seq(0,1.1*M,.005),col=c("blu e","orange"), xlab="KS-Statistic for Normal Samples", main=paste("Sample size n =",n," , Nsim =",Nsim)) crit.alpha=quantile(sim.out.norm,1-alpha) Sample size n = 100 , Nsim = 10000 KS-Statistic for Normal Samples F re q u e n cy 0.00 0.05 0.10 0.15 0.20 0 4 0 0 8 0 0 1 2 0 0  0.05 KS-Statistic for Kinked Normal Samples F re q u e n cy 0.00 0.05 0.10 0.15 0.20 0 2 0 0 4 0 0 6 0 0 8 0 0  0.6446 Sample size n = 100 , Nsim = 10000 AD-Statistic for Normal Samples F re q u e n cy 0.0 0.5 1.0 1.5 2.0 0 5 0 0 1 5 0 0 2 5 0 0  0.05 AD-Statistic for Kinked Normal Samples F re q u e n cy 0.0 0.5 1.0 1.5 2.0 0 5 0 0 1 5 0 0  0.139 Clearly the AD-Statistic shows much less power than the KS-statistics for such kinked normal samples. To understand why, we derive the CDF of this kinked normal random variable V, defined as V = Z when |Z| > .2 and V = 0 when |Z|  .2. Here Z is a standard normal random variable. Note that P( V  x )=P( Z  x )=(x) when x  -.2. For x  .2 we have P(V x)=P( V-.2 )+P( V = 0) + P( .2 V  x) = P(Z-.2)+P( -.2 Z  .2) + P( .2 Z  x) )=(x) For -.2  x < 0 we have P(V x)=P(V-.2) = (-.2) And for 0  x  .2 we have P(V x)=P(V-.2) + P(V = 0 ) = (Z -.2)+P( -.2 Z  .2) = (.2) The plot of this CDF is given below followed by the empirical CDF of a sample of size n=10000 from this CDF. The discrepancy between this kinked CDF and the standard normal CDF (superimposed) is local near zero (within the interval [-.2,.2]), not in the tails. Thus the KS-statistic should be sensitive to it while the AD-statistic will be less effective since discrepancies don’t extend over a wider range and thus do not contribute much in the integration, and discrepancies are not in the tails. -3 -2 -1 0 1 2 3 0 .0 0 .2 0 .4 0 .6 0 .8 1 .0 x P V  x -4 -2 0 2 4 0 .0 0 .2 0 .4 0 .6 0 .8 1 .0 x ki n ke d C D F