作者:

标题:

R-code：

#R-code for Zhong, PS and Chen, SX (2011) Tests for High #Dimensional Regression Coefficients with Factorial Designs. Journal of #the American Statistical Association, 106, 260-274. 
TestCoef<-function(p,n,rhop,eta,B)
{
m<-2*p-1
count<-1
non0<-sqrt(eta/5)
beta<-c(rep(0,p-5),rep(non0,5))
set.seed(10000)
mu<-runif(p,2,3)
rho1<-runif(rhop,0,1)
rho<-c(rho1,rep(0,p-rhop))
Gamma<-cbind(diag(rep(rho[1],p)),matrix(0,p,m-p))
for (i in 2:(p-1))
{Gamma<-Gamma+cbind(matrix(0,p,i-1),diag(rep(rho[i],p)),matrix(0,p,m-p-i+1))}
Gamma<-Gamma+cbind(matrix(0,p,m-p),diag(rep(rho[p],p)))
Tnvec<-NULL
repeat
{   set.seed(Sys.time())
    z<-matrix(rnorm(n*m,0,1),n,m)
    X<-z%*%t(Gamma)+matrix(rep(mu,n),n,p,byrow=T)
    Y<-X%*%beta+rnorm(n,0,2)
    TnStat1<-function(com)
    {
     k<-com[1]
     l<-com[2]
     s<-com[3]
     t<-com[4]
     Tn1<-sum((X[k,]-X[l,])*(X[s,]-X[t,]))*(Y[k]-Y[l])*(Y[s]-Y[t])/4
     Tn2<-sum((X[k,]-X[s,])*(X[l,]-X[t,]))*(Y[k]-Y[s])*(Y[l]-Y[t])/4
     Tn3<-sum((X[k,]-X[t,])*(X[s,]-X[l,]))*(Y[k]-Y[t])*(Y[s]-Y[l])/4
     Tn<-(Tn1+Tn2+Tn3)/3
     return(Tn)
     }
  XXn<-X%*%t(X)
  Y1<-sum(diag(t(X)%*%X))/n
  Y3<-(sum(XXn)-Y1*n)/(n*(n-1))
  Y2<-(sum(XXn^2)-sum(diag(XXn^2)))/(n*(n-1))
  XXn2<-XXn%*%XXn
  diagXX<-rep(diag(XXn),each=n-1)
  offdiagXX<-as.logical(lower.tri(XXn)+upper.tri(XXn))
  VecOffdiagXX<-matrix(XXn,n^2,1)[offdiagXX]
  Y4<-(sum(XXn2)-sum(diag(XXn2))-2*sum(diagXX*VecOffdiagXX))/(n*(n-1)*(n-2))
  Y5<-((n*(n-1)*Y3)^2-2*n*(n-1)*Y2-4*n*(n-1)*(n-2)*Y4)/(n*(n-1)*(n-2)*(n-3))
  trSigma2<-Y2-2*Y4+Y5
  Tn<-sum(combn(n,4,TnStat1))/choose(n,4)
  sigma2<-var(Y)
  varTn<-2*trSigma2*sigma2^2/(n*(n-1))
  Tnvec<-c(Tnvec,Tn/sqrt(varTn))
  count<-count+1
  if (count>B) break;
}
  rej2<-sum(Tnvec>qnorm(0.95,0,1))/B
  return(list(Ttest=rej2,eta=eta,Tn=Tnvec))
}

 p<-480
 n<-50
 rhop<-10
 result11try<-TestCoef(p,n,rhop,0,1)