10001036 : Study 2

rm(list=ls());alpha=.05
Prob=.001;lambda=1.01;
p.max=20;##p is number of IVs including intercept
q.max=8; ##q is number of auxiliary variables

ProbUR <- function(N,p,q,lambda,alpha){
## Probability of that scale ratio of unrestricted model confidence region to
##   restricted model confidence region is larger than lambda
## Df of unrestricted model is N-p-q while df of restricted model is N-p.
## Alpha is the type I error of the confidence region.
  F=(1/(lambda)^2)*(qf(1-alpha,p,N-p-q)/qf(1-alpha,p,N-p))*(N-p)/q - (N-p-q)/q;
  pf(F,q,N-p-q)
}

N.min <- function(Prob,p,q,lambda,alpha){
  N <- p+q+10:500;
  min(N[ProbUR(N,p,q,lambda,alpha)<Prob])
}

N.matrix <- function (p.max=20,q.max=8,Prob=.01,lambda=1.01,alpha=.05){
  N.set<-matrix(nrow=p.max,ncol=q.max);
  rownames(N.set)=paste('p=',1:p.max,sep='');
  colnames(N.set)=paste('q=',1:q.max,sep='');
  N.seq<-N.set;
  for (p in 1:p.max)
    for (q in 1:q.max){
      N.set[p,q]<- N.min(Prob,p,q,lambda,alpha);
      N.seq[p,q]<- N.min(Prob/q,p+q-1,1,lambda^(1/q),alpha);
  ##  N.seq[p,q]<- N.min(Prob/q,p,1,lambda^(1/q),alpha);
  }
  plot(c(row(N.seq),row(N.set)),c(N.seq,N.set)
    ,pch=as.character(c(col(N.seq),col(N.set)))
    ,col=c(rep('grey',p.max*q.max),rep('black',p.max*q.max))
    ,xlab='Number of IVs in restricted model'
    ,ylab='Sample Size'
    ,main=paste(sep='\n','Grey:Sequential Enter  Dark:Set Enter'
      ,sprintf('Prob:%g  Lambda:%g alpha:%g',Prob,lambda,alpha)
     )
    ,ylim=c(0,min(500,max(N.seq,N.set)))
  )

pmin(N.seq,N.set)
  ##require(R2HTML);require(Rpad);.HTML.file = "";
  ##HTML(pmin(N.seq,N.set)
  ##  ,caption=sprintf('Prob:%g  Lambda:%g alpha:%g',Prob,lambda,alpha));

}

N.matrix(p.max,q.max,Prob,lambda,alpha);
N.matrix(p.max,q.max,0.001,1.05,alpha);

10001036

Study 2

Categories

Meta

Blog Search

Recent Comments

链接表

Recent Posts

Archives

Pages