umai88
2012-May-29 06:06 UTC
[R] need help to find type I error rate for modified F statistic
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]]
Rui Barradas
2012-May-29 11:26 UTC
[R] need help to find type I error rate for modified F statistic
Hello, Without trying it, I very much doubt that your code runs: First you use 'asim' and create it after. So I've moved asim <- 1000 up some lines. And it still doesn't work. 1. Run your own code before posting it. 2. Learn to use indentation. 3. Learn to use spaces to make your code easier to read. Hope this helps, Rui Barradas P.S. At a glance, it seems you are dividing 'num' by J and only then subtracting 1. Missing parenthesis? Impossible to say. Follow the above first. R.B. Em 29-05-2012 07:06, umai88 escreveu:> 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]] > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.