Hello everyone, I want to calculate type I error rate for modified F
statistic for one way robust anova. I need to find the j group trimmed
mean and winsorized sum of squared deviations. Here I attached my code for
j=2 to make it simple. Originally I have j=4. Hope you can help. I need to
run it for 1000 times
My problem is:
i) the value of F-test obtain from my simulation below is in negative
value..There might be something wrong in my coding
ii) after obtain F value, how i can proceed to obtain type I error rate?
h=0
g=0
n=15
alpha=0.1
k=floor(alpha*n)+1
r=k-(alpha*n)
i=k+1
m=n-k
J=2
trim1<-rep(NA,asim)
trim2<-rep(NA,asim)
pw<-rep(NA,asim)
qw<-rep(NA,asim)
px<-rep(NA,asim)
qx<-rep(NA,asim)
ssd1<-rep(NA,asim)
ssd2<-rep(NA,asim)
madN1<-rep(NA,asim)
madN2<-rep(NA,asim)
lo.w<-rep(NA,asim)
up.w<-rep(NA,asim)
lo.x<-rep(NA,asim)
up.x<-rep(NA,asim)
hw<-rep(NA,asim)
hx<-rep(NA,asim)
xbar.t<-rep(NA,asim)
num<-rep(NA,asim)
df2<-rep(NA,asim)
denom<-rep(NA,asim)
F<-rep(NA,asim)
asim<-1000
for(j in 1:asim)
{
print(j)
set.seed(j)
a<-rnorm(15,0,1)
b<-rnorm(15,0,1)
w<-a*exp(h*a^2/2)
x<-b*exp(h*b^2/2)
matw<-sort(w)
matx<-sort(x)
trim1[j]=1/((1-2*alpha)*n)*(sum(matw[i:m])+r*(matw[k]+matw[n-k+1]))
trim2[j]=1/((1-2*alpha)*n)*(sum(matx[i:m])+r*(matx[k]+matx[n-k+1]))
madN1[j]<-mad(w)
madN2[j]<-mad(x)
lo.w[j]<-length(which(w-median(w)<(-2.24*madN1[j])))
up.w[j]<-length(which(w-median(w)>(2.24*madN1[j])))
lo.x[j]<-length(which(x-median(x)<(-2.24*madN2[j])))
up.x[j]<-length(which(x-median(x)>(2.24*madN2[j])))
pw[j]<-up.w[j]+2
qw[j]<-n-up.w[j]-1
px[j]<-up.x[j]+2
qx[j]<-n-up.x[j]-1
ssd1[j]<-lo.w[j]*(matw[lo.w[j]+1]-trim1[j])^2 +
(matw[pw[j]:qw[j]]-trim1[j])^2 + up.w[j]*(matw[n-up.w]-trim1[j])^2 -
((lo.w)*(matw[lo.w+1]-trim1[j])+ (up.w)*(matw[n-up.w]-trim1[j]))^2
ssd2[j]<-lo.x[j]*(matx[lo.x[j]+1]-trim2[j])^2 +
(matx[px[j]:qx[j]]-trim2[j])^2 + up.x[j]*(matx[n-up.x]-trim2[j])^2 -
((lo.x)*(matw[lo.x+1]-trim2[j])+ (up.x)*(matx[n-up.x]-trim2[j]))^2
hw[j]<-n-lo.w[j]-up.w[j]
hx[j]<-n-lo.x[j]-up.x[j]
H[j]=hw[j]+hx[j]
xbar.t[j]<-(hw[j]*trim1[j]+hx[j]*trim2[j])/H[j]
num[j]<-( (trim1[j]-xbar.t[j])^2 + (trim2[j]-xbar.t[j])^2)^2 )/J-1
df2[j]<-H[j]-J
denom[j]<-(ssd1[j]+ssd2[j])/df2[j]
F[j]<-num[j]/denom[j]
}
Graduate student,
Department of Mathematics& Statistics, UPM
--
View this message in context: http://r.789695.n4.nabble.com/need-help-to-find-type-I-error-rate-for-modified-F-statistic-tp4631661.html
Sent from the R help mailing list archive at Nabble.com.
[[alternative HTML version deleted]]