[R] fixed trimmed mean for j-group
Hello...i want to find the empirical rate for type 1 error using fixed
trimmed mean. To make it easy, i'm referring to journal given by this
website
http://www.academicjournals.org/ajmcsr/PDF/pdf2011/Yusof%20et%20al.pdf.
I already run the programme and there is no error in it but i got zero for
the empirical rate of type 1 error. The empirical rate for the type 1 error
given in the journal should lied between 0.025 and 0.07. Below is my
programme. Hope anyone can help me. Thank you.
## use for data transformation
g=0
h=0
w-a*exp(h*a^2/2)
x-b*exp(h*b^2/2)
y-c*exp(h*c^2/2)
z-d*exp(h*d^2/2)
g=0.5 and h=0
g=0.5 and h=0.5
w-(exp(g*a)-1)/g*(exp((h*a^2)/2))
x-(exp(g*b)-1)/g*(exp((h*b^2)/2))
y-(exp(g*c)-1)/g*(exp((h*c^2)/2))
z-(exp(g*d)-1)/g*(exp((h*d^2)/2))
FIXED SYMMETRIC TRIMMED MEAN
###
asim-5000
pv-rep(NA, asim)
for(j in 1:asim)
{
print(j)
set.seed(j)
n1=15
n2=15
n3=15
n4=15
miu=0
sd1=1
sd2=1
sd3=1
sd4=1
a=rnorm(n1,miu,sd1)
b=rnorm(n2,miu,sd2)
c=rnorm(n3,miu,sd3)
d=rnorm(n4,miu,sd4)
## data transformation
g=0
h=0
w-a*exp(h*a^2/2)
x-b*exp(h*b^2/2)
y-c*exp(h*c^2/2)
z-d*exp(h*d^2/2)
mat1-sort(w)
mat2-sort(x)
mat3-sort(y)
mat4-sort(z)
alpha=0.15
k1=floor(alpha*n1)+1
k2=floor(alpha*n2)+1
k3=floor(alpha*n3)+1
k4=floor(alpha*n4)+1
r1=k1-(alpha*n1)
r2=k2-(alpha*n2)
r3=k3-(alpha*n3)
r4=k4-(alpha*n4)
## j-group trimmed mean
e1=k1+1
f1=n1-k1
e2=k2+1
f2=n2-k2
e3=k3+1
f3=n3-k3
e4=k4+1
f4=n4-k4
trim1=1/((1-2*alpha)*n1)*(sum(mat1[e1:f1]) + r1*(mat1[k1]+mat1[n1-k1+1]))
trim2=1/((1-2*alpha)*n2)*(sum(mat2[e2:f2]) + r2*(mat2[k2]+mat2[n2-k2+1]))
trim3=1/((1-2*alpha)*n3)*(sum(mat3[e3:f3]) + r3*(mat3[k3]+mat3[n3-k3+1]))
trim4=1/((1-2*alpha)*n4)*(sum(mat4[e4:f4]) + r4*(mat4[k4]+mat4[n4-k4+1]))
## sample winsorized mean
x1=(1-r1)*mat1[k1+1]+r1*mat1[k1]
x2=(1-r2)*mat2[k2+1]+r2*mat2[k2]
x3=(1-r3)*mat3[k3+1]+r3*mat3[k3]
x4=(1-r4)*mat4[k4+1]+r4*mat4[k4]
u1=(1-r1)*mat1[n1-k1]+r1*mat1[n1-k1+1]
u2=(1-r2)*mat2[n2-k2]+r2*mat2[n2-k2+1]
u3=(1-r3)*mat3[n3-k3]+r3*mat3[n3-k3+1]
u4=(1-r4)*mat4[n4-k4]+r4*mat4[n4-k4+1]
win1=1/n1* (sum(mat1[e1:f1])+k1*(x1+u1))
win2=1/n2* (sum(mat2[e2:f2])+k2*(x2+u2))
win3=1/n3* (sum(mat3[e3:f3])+k3*(x3+u3))
win4=1/n4* (sum(mat4[e4:f4])+k4*(x4+u4))
## g-winsorized sum of squared deviations
ssd1=sum((mat1[e1:f1]-win1)^2) + k1*((mat1[k1+1]-win1)^2 +
(mat1[n1-k1]-win1)^2)
ssd2=sum((mat2[e2:f2]-win2)^2) + k2*((mat2[k2+1]-win2)^2 +
(mat2[n2-k2]-win2)^2)
ssd3=sum((mat3[e3:f3]-win3)^2) + k3*((mat3[k3+1]-win3)^2 +
(mat3[n3-k3]-win3)^2)
ssd4=sum((mat4[e4:f4]-win4)^2) + k4*((mat4[k4+1]-win4)^2 +
(mat4[n4-k4]-win4)^2)
## calculate f statistic
J=4
h1=n1-2*floor(alpha*n1)
h2=n2-2*floor(alpha*n2)
h3=n3-2*floor(alpha*n3)
h4=n4-2*floor(alpha*n4)
H=h1+h2+h3+h4
xt=(h1*trim1+h2*trim2+h3*trim3+h4*trim4)/H
num= ((trim1-xt)^2+(trim2-xt)^2+(trim3-xt)^2+(trim4-xt)^2)/(J-1)
denom= (ssd1+ssd2+ssd3+ssd4)/(H-J)
ft=num/denom
pv[j]=1-pf(ft,(J-1),(H-J))
}
mean(pv0.05)
##
Graduate student of Universiti Putra Malaysia
--
View this message in context:
http://r.789695.n4.nabble.com/fixed-trimmed-mean-for-j-group-tp4633327.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@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.
[R] obtain type I error from modified F statistic
i want to find type I error rates using modified F statistic..below is my R
code..
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=4
trim1-rep(NA,asim)
trim2-rep(NA,asim)
trim3-rep(NA,asim)
trim4-rep(NA,asim)
pw-rep(NA,asim)
qw-rep(NA,asim)
px-rep(NA,asim)
qx-rep(NA,asim)
py-rep(NA,asim)
qy-rep(NA,asim)
pz-rep(NA,asim)
qz-rep(NA,asim)
ssd1-rep(NA,asim)
ssd2-rep(NA,asim)
ssd3-rep(NA,asim)
ssd4-rep(NA,asim)
madN1-rep(NA,asim)
madN2-rep(NA,asim)
madN3-rep(NA,asim)
madN4-rep(NA,asim)
lo.w-rep(NA,asim)
up.w-rep(NA,asim)
lo.x-rep(NA,asim)
up.x-rep(NA,asim)
lo.y-rep(NA,asim)
up.y-rep(NA,asim)
lo.z-rep(NA,asim)
up.z-rep(NA,asim)
hw-rep(NA,asim)
hx-rep(NA,asim)
hy-rep(NA,asim)
hz-rep(NA,asim)
H-rep(NA,asim)
xbar.t-rep(NA,asim)
num-rep(NA,asim)
df2-rep(NA,asim)
denom-rep(NA,asim)
F-rep(NA,asim)
F.table-rep(NA,asim)
asim-10
for(j in 1:asim)
{
print(j)
set.seed(j)
a-rnorm(15,0,1)
b-rnorm(15,0,1)
c-rnorm(15,0,1)
d-rnorm(15,0,1)
w-a*exp(h*a^2/2)
x-b*exp(h*b^2/2)
y-c*exp(h*c^2/2)
z-d*exp(h*d^2/2)
matw-sort(w)
matx-sort(x)
maty-sort(y)
matz-sort(z)
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]))
trim3[j]=1/((1-2*alpha)*n)*(sum(maty[i:m])+r*(maty[k]+maty[n-k+1]))
trim4[j]=1/((1-2*alpha)*n)*(sum(matz[i:m])+r*(matz[k]+matz[n-k+1]))
madN1[j]-mad(w)
madN2[j]-mad(x)
madN3[j]-mad(y)
madN4[j]-mad(z)
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])))
lo.y[j]-length(which(y-median(y)(-2.24*madN3[j])))
up.y[j]-length(which(y-median(y)(2.24*madN3[j])))
lo.z[j]-length(which(z-median(z)(-2.24*madN4[j])))
up.z[j]-length(which(z-median(z)(2.24*madN4[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
py[j]-up.y[j]+2
qy[j]-n-up.y[j]-1
pz[j]-up.z[j]+2
qz[j]-n-up.z[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
ssd3[j]-lo.y[j]*(maty[lo.y[j]+1]-trim3[j])^2 +
(maty[py[j]:qy[j]]-trim3[j])^2 + up.y[j]*(maty[n-up.y]-trim3[j])^2 -
((lo.y)*(matw[lo.y+1]-trim3[j])+ (up.y)*(maty[n-up.y]-trim3[j]))^2
ssd4[j]-lo.z[j]*(matz[lo.z[j]+1]-trim4[j])^2 +
(matz[pz[j]:qz[j]]-trim4[j])^2 + up.z[j]*(matz[n-up.z]-trim4[j])^2 -
((lo.z)*(matw[lo.z+1]-trim4[j])+ (up.z)*(matz[n-up.z]-trim4[j]))^2
hw[j]-n-lo.w[j]-up.w[j]
hx[j]-n-lo.x[j]-up.x[j]
hy[j]-n-lo.y[j]-up.y[j]
hz[j]-n-lo.z[j]-up.z[j]
H[j]=hw[j]+hx[j]+hy[j]+hz[j]
xbar.t[j]-(hw[j]*trim1[j]+hx[j]*trim2[j]+hy[j]*trim3[j]+hz[j]*trim4[j])/H[j]
num[j]-( (trim1[j]-xbar.t[j])^2 + (trim2[j]-xbar.t[j])^2 +
(trim3[j]-xbar.t[j])^2+ (trim4[j]-xbar.t[j])^2 )/J-1
df2[j]-H[j]-J
denom[j]-(ssd1[j]+ssd2[j]+ssd3[j]+ssd4[j])/df2[j]
F[j]-num[j]/denom[j]
}
--
View this message in context:
http://r.789695.n4.nabble.com/obtain-type-I-error-from-modified-F-statistic-tp4631660.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@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.
[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]]
__
R-help@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.
[R] code to iterate function apply to matrix
I got this code below and i want to repeat the loop for 100 times..
x-rnorm(60)
mat1-matrix(x,nrow=15,ncol=4)
trim-numeric(ncol(mat1))
win-numeric(ncol(mat1))
ssd-numeric(ncol(mat1))
for(j in 1:ncol(mat1))
{
n=length(mat1[,j])
alpha=0.1
k=floor(alpha*n)+1
r=k-(alpha*n)
i=k+1
m=n-k
y1-sort(mat1[,j])
y-y1[i:m]
x.low=(1-r)*y1[k+1]+r*y1[k]
x.upp=(1-r)*y1[n-k]+r*y1[n-k+1]
trim[j] =1/((1-2*alpha)*n)*(sum(y)+r*(y1[k]+y1[n-k+1]))
win[j]=1/n*(sum(y)+k*(x.low+x.upp))
ssd[j]-sum((y-win[j])**2)+k*( (y1[k+1]-win[j])**2 + (y1[n-k]-win[j])**2 )
}
trim.mean-matrix(trim, nrow=1)
win.mean-matrix(win, nrow=1)
sum.sq.dev-matrix(ssd, nrow=1)
--
View this message in context:
http://r.789695.n4.nabble.com/code-to-iterate-function-apply-to-matrix-tp4630221.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@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.
