#------------------------------------------------------------------------# # Modified on : July 29, 2013 # Created on : April 1, 2012 # AUTHOR : Hyokyoung Grace Hong jointly with Xuming He and Lan Wang # AFFILIATION : Michigan State University # EMAIL : hhong@stt.msu.edu # Simulate Example 1 case (1b) in He, Wang, and Hong (2013, Annals of Statistics, # QUANTILE-ADAPTIVE MODEL-FREE VARIABLE SCREENING FOR HIGH-DIMENSIONAL HETEROGENEOUS DATA) # This program computes the minimum model size(denoted by R), that is, # the smallest number of covariates needed to ensure that all the active # variables are selected. # The proportion of active variables (denoted by S) can be simply obtained using # the threshold [n/log(n)]. # This program also computes L2-based nonlinear screening and SIRS for comparison. # R: 2.15.1 #-------------------------------------------------------------------------# ### load packages library(quantreg) library(VGAM) library(MASS) #------------------------------------------------------------------------# ### Compute Quantile Adaptive sure Independence Screening (QaSIS) QaSIS<-function(Y,X,tau) { p=dim(X)[2] n=length(Y) s<-matrix(0,nrow=p,ncol=1) bfit<-matrix(0,nrow=n,ncol=p) Y<-Y-quantile(Y,tau) #centered y X<-apply(X,2, function(x) (x-mean(x))/sd(x)) for(j in 1:p){ x0<-X[,j] knots=quantile(x0,c(1/3,2/3)) a=bs(x0, knots=knots,degree=1) b=rq(Y~a,tau=tau) s[j,] <-sum((b$fitted)^2)/n #avg of f^2 for each j bfit[,j]<-b$fitted #pi(x)%*%bhat } list(rank=(p+1)-rank(s),fit=bfit,w.k=s) } ### Compute L2-based Nonparametric Independence Screening (NIS) NIS<-function(Y,X) { p=dim(X)[2] n=length(Y) s<-matrix(0,nrow=p,ncol=1) bfit<-matrix(0,nrow=n,ncol=p) Y<-Y-mean(Y) #centered y X<-apply(X,2,function(x) (x-min(x))/(max(x)-min(x))) for(j in 1:p){ x0<-X[,j] knots=quantile(x0,c(1/3,2/3)) a=bs(x0, knots=knots,degree=1) b=lm(Y~a) s[j,] <-sum((b$fitted)^2)/n #avg of f^2 for each j bfit[,j]<-b$fitted #pi(x)%*%bhat } list(rank=(p+1)-rank(s),fit=bfit,w.k=s) } ### Compute Sure Independent Ranking and Screening (SIRS) by Zhu et al. (2011) SIRS<-function(X,Y){ p=dim(X)[2] n=length(Y) w<-matrix(0, nrow=p, ncol=1) X<-apply(X,2, function(x) (x-mean(x))/sd(x)) for(k in 1:p){ w.k.j<-NULL for(j in 1:n){ s<-(t(X[,k])%*%(1*(Y