#Generate 10,000 samples of size 10 for two different groups, each with the same variance
X1_values=matrix(rnorm(100000,mean=0,sd=1),ncol=10)
X2_values=matrix(rnorm(100000,mean=0,sd=1),ncol=10)
#Compute variances, F statistic and p-value
V1=apply(X1_values,1,var)
V2=apply(X2_values,1,var)
F=V1/V2
pval=1-pf(F,9,9)
hist(pval)
#Proportion of simulations for which the pvalue is less than alpha
sum(pval<0.05)/10000

#Repeat above simulation with sigma1=1.2
X1_values=matrix(rnorm(100000,mean=0,sd=1.2),ncol=10)
X2_values=matrix(rnorm(100000,mean=0,sd=1),ncol=10)
#Compute variances, F statistic and p-value
V1=apply(X1_values,1,var)
V2=apply(X2_values,1,var)
F=V1/V2
pval=1-pf(F,9,9)
hist(pval)
#Proportion of simulations for which the pvalue is less than alpha
sum(pval<0.05)/10000

#Repeat above simulation with sigma1=2
X1_values=matrix(rnorm(100000,mean=0,sd=2),ncol=10)
X2_values=matrix(rnorm(100000,mean=0,sd=1),ncol=10)
#Compute variances, F statistic and p-value
V1=apply(X1_values,1,var)
V2=apply(X2_values,1,var)
F=V1/V2
pval=1-pf(F,9,9)
hist(pval)
#Proportion of simulations for which the pvalue is less than alpha
sum(pval<0.05)/10000

#Repeat above simulation with sigma1=2 and n1=n2=20
X1_values=matrix(rnorm(200000,mean=0,sd=2),ncol=20)
X2_values=matrix(rnorm(200000,mean=0,sd=1),ncol=20)
#Compute variances, F statistic and p-value
V1=apply(X1_values,1,var)
V2=apply(X2_values,1,var)
F=V1/V2
pval=1-pf(F,9,9)
hist(pval)
#Proportion of simulations for which the pvalue is less than alpha
sum(pval<0.05)/10000

#Robustness
#Confirm that departures from normality are hard to detect for small samples
X_values=rchisq(10,5)
qqnorm(X_values)
qqline(X_values)
xvals=seq(0,15,by=0.1)
plot(xvals,dchisq(xvals,5),type="l")

#Generate 10,000 chi-squared samples of size 10 for two different groups, each with the same variance
X1_values=matrix(rchisq(100000,df=5),ncol=10)
X2_values=matrix(rchisq(100000,df=5),ncol=10)
#Compute variances, F statistic and p-value
V1=apply(X1_values,1,var)
V2=apply(X2_values,1,var)
F=V1/V2
pval=1-pf(F,9,9)
hist(pval)
#Proportion of simulations for which the pvalue is less than alpha
sum(pval<0.05)/10000

#Generate 10,000 chi-squared samples of size 10 for two different groups, but sigma1=1.2*sigma2
X1_values=matrix(rchisq(100000,df=7.2),ncol=10) #df=5*1.2^2
X2_values=matrix(rchisq(100000,df=5),ncol=10)
#Compute variances, F statistic and p-value
V1=apply(X1_values,1,var)
V2=apply(X2_values,1,var)
F=V1/V2
pval=1-pf(F,9,9)
hist(pval)
#Proportion of simulations for which the pvalue is less than alpha
sum(pval<0.05)/10000

#Generate 10,000 chi-squared samples of size 10 for two different groups, but sigma1=1.2*sigma2
X1_values=matrix(rchisq(100000,df=20),ncol=10) #df=5*2^2
X2_values=matrix(rchisq(100000,df=5),ncol=10)
#Compute variances, F statistic and p-value
V1=apply(X1_values,1,var)
V2=apply(X2_values,1,var)
F=V1/V2
pval=1-pf(F,9,9)
hist(pval)
#Proportion of simulations for which the pvalue is less than alpha
sum(pval<0.05)/10000

#Graphing power curves
mu0=100
sigma=4
n=10
alpha=0.05
mua=seq(90,110,by=0.01)
#plot(mua,pnorm(-qnorm(1-alpha)+sqrt(n)*(mu0-mua)/sigma),type="l",xlab=expression(mu),ylab="Power")
#plot(mua,1-pnorm(qnorm(1-alpha)+sqrt(n)*(mu0-mua)/sigma),type="l",xlab=expression(mu),ylab="Power")
plot(mua,pnorm(-qnorm(1-alpha/2)+sqrt(n)*(mu0-mua)/sigma)+1-pnorm(qnorm(1-alpha/2)+sqrt(n)*(mu0-mua)/sigma),type="l",xlab=expression(mu),ylab="Power")
abline(h=alpha,col="red",lty=3)
abline(v=100,col="red",lty=3)