非参数回归和相关统计检验

Generalized additive model （简称GAM）在1986年由Hastie & Tibshirani发表以后，一直被应用于许多领域，如环境科学、医学、生物统计，在方法上也得到发展，如Amato等人（2002）提出的基于傅立叶变换最优逼近模式估计，Sardy（2003）基于非线性小波估计方法，范剑青2004年提出的基于核函数的广义可加模型估计以及基于GLR的统计推断，对可加模型的估计和统计推断做了进一步的发展，而且在应用上该方法也被拓展到经济金融领域，Wensui Liu在2007年提出了基于GAM的信用评分模型改进无疑是该方法在金融实践中的一个成功的尝试，在不久的将来也许非参数计量将会是下一代计量理论体系的核心补丁。
下面说说GAM方法的优势，首先它没有参数方法对先验分布的限制前提（如normalassumption、homoscedasticity assumption），其次在处理非线性模式的能力要远远强于参数模型，当然新的模式出现也会给我们的应用增添一些疑惑，比如在可加成分是否显著或如何判断自变量与因变量间是线性还是非线性模式，基于这个问题，理论学者和应用大牛们也都分别给出了相应的方法和尝试，比如范剑青的GLR检验、profile least-squares estimation，Hastie的Back-Fitting and Local Scoring Algorithms，他们的方法本质上都是基于非参数部分的局部分解以及整体模式残差最小的思想，如果因变量是二分变量估计方法相应替换成似然方法便是。
接下来给出相应的估计和检验方法，之前博客里已经给出范剑青的核方法的半变系数非参数模型的估计和检验程序（理论可以参考Profile Likelihood inferences on Semiparametric Varying-coeffcient Partially Liner Models），接下来介绍基于SAS集成的程序（更灵活简单）。
先介绍下GAM 基本的概念。一般线性模式可如下所示：
E(Y)=b0+b1*X1+b2*X2+…+bp*Xp
但在 GAM 下，模式相应改为：
E(Y)=s0+s1(X1)+s2(X2)+…+sp(Xp)
其中 si(Xi) 为smooth function（平滑函数）。值得注意的是这个方程式没有确定的模式，而以 non-parametric 的形态去估计。SAS的过程布中采用「backfitting algorithm」。由于 GAMscore equation 没有close form，所以必须依赖backfitting algorithm回溯拟合或构建估出si(Xi)的局
部核函数。
接下来采用Dong Xiang 在 SUGI 26发表的文章里的案例简要说明Proc GAM的使用。Xiang 引用了 Bell et al. 在1989 年的一份有关推背的相关数据，其反应变量是个二分类变量（有(1)和无(0)）。三個预测变量分别是 Age, StartVert和 NumVert。
sas估计程序：
PROC GAM data=kyphosis;
model kyphosis = spline(NumVert,df=3)
spline(Age,df=3)
spline(StartVert,df=3)
/dist=logist;
output out=estimate pred uclm lclm;
run;
非参数logistic回归模型

其中spline是平滑方式，可用loess和spline2代替，三者的解释为：spline的平滑方式为cubic spline，即三次样条逼近，loess为局部估计（即为核估计），spline2是bivariate thin-plate smoothing；语法分别为：SPLINE(variable<,df=number>)，loess(variable<,df=number>)，SPLINE2(variable,variable<,df=number>)，其中df为平滑统计量的自由度，缺省为4，可以自己指定，在option添加GCV（广义交叉确认）时不需在平滑函数中指定自由度，可以自动确定各变量的最优自由度。如果有的变量与因变量为线性关系，也可以通过param（variable（s））指定参数变量，如果一个回归模型中既有param又有loess部分的话，模型即为半参数模型，proc gam过程目前还不支持系数为非参数函数的模型。当因变量为二分类变量时，模型变是logistic additive model，需要通过制定link function进行拟合，利用dist option指定依赖变量的分布模式，即设定dist=logist,如果link function为dist=log时，则模式即为poission model。文中kyphosis为二分类变量，故而指定logist为连接函数。
注：dist可以是GAUSSIAN, BINOMIAL, BINARY, GAMMA, or POISSON.
另外OUTPUT statement 可以将参数估计、非参数平滑预测值、目标变量的预测概率和其 95% 置信区间分别给出（pred uclm lclm）并产生新变量。估计的数据可以用来做图。
非参数logistic回归模型

先看下输出结果，第一部分是summary statistics，提供一些参数估计过程的信息，如下：
非参数logistic回归模型

目标变量为二项分布（binomial），连接函数为logist，局部评分迭代9次等。接下来看下第二部分的输出结果：
非参数logistic回归模型

这部分包含三张报表。第一个表示模式的线性部分的参数估计值（包含非参数的线性分解部分，具体会在后面介绍）。第二个是用 GAM 跑出的 smoother 即窗宽值。第三个表比较特别。此处利用 deviance analysis检验比较full model（包含全部 smoothing function）和 reduced model（少掉某一个smoothing function）那個统计量的检验结果（这正是GLR的核心思想，范剑青在2005的一篇文献中证明，这个似然统计量遵从wilks现象，即服从卡方分布，在具体应用中可根据卡方分布拟合及Bootstrap具体判断参数的显著推断）。所以零假设 H0: reduce model vs H1: full model。如果的 p-value小于0.05 則表示拒绝 reduced model。反之则表示不拒绝 reduce model。以 Age 和 StartVert 来看，两者的 p-value 分別是 0.0009 和 0.0350，皆小于0.05。所以可以同时分別少了这两个变量的 reduced model（对比参数模型中b=0的假设）。而 NumVert 的 p-value 为0.3311 > 0.05，表示不拒绝少了 NumVert 的 smoothing function 的reduced model，也就是说f(NumVert) 在模式中不显著。
很多人可能会疑问，程序中设定的变量均为非参数形式为何输出的结果会有回归结果的parameter估计值和非参数部分的smoothing parameter，这得从gam过程的估计原理说起：当模型设定为下式时
model y = param(z) spline(x1) spline(x2)
则表示估计的模型为y = a0 + a1*z + s1(x1) + s2(x2)，而GAM过程实际估计的是y = a0 + a1*z + b1*x1+f1(x1) + b2*x2 + f2(x2)，其中s1(x1) = b1*x1 + f1(x1), GAM 分离出线性趋势项和非参数的平滑函数f1(x1),其中线性部分决定自变量的主要趋势，非参数部分为调整项。
所以上述输出结果中 “Regression Model Analysis” with “Parameter Estimates” 就是a0 (labeled Intercept)的估计值,，a1 (labeled z),，b1 (labeled L_X1 for the linear part of x1，例子中的age)，and b2 (labeled L_X2，例子中的startvert)等等。在程序的output设定中的Pred会在结果数据集中自动输出非参数的估计部分：
P_Y = a0 + a1*z + b1*x1+f1(x1) + b2*x2 + f2(x2)
P_X1 = f1(X1)
P_X2 = f1(X2)
故而，非参数部分的full partial prediction即为：b1*x1+p_X1às(x1)，b2*x2+p_x2às(x2)；总体的预测必须在data步完成，或在proc步中添加score语句。

另附半参数回归估计和基于bootstrap参数检验程序：

%%**程序说明：SAS IML的半变系数部分线性模型的非参数估计其中，x1,x2为可变部分，z1-z3为线性部分。Reference:
Profile Likelihood Inferences on Semiparametric Varying-Coefficient Partially Linear Models
by Jianqing Fan****
%let x=%bquote(x1 x2);
%let z=%bquote(z1 z2 z3);
%let dataset=simul_dataset;
%let varnum=2;
proc iml;
use &dataset;
read all into u var{u};
read all into x var{&x};
read all into z var{&z};
read all into y var{y};
close &dataset;
sort &dataset out=data_set0 by u;
use data_set0;
read all into u1 var{u};
n=nrow(x);
x1=j(n,1,1);
x=x1||x;
%let vary_var=%eval_r(&varnum+1);
define kernel function
start kde(o);
kernel=0.75*(1-o##2)#(abs(o<=1));
return(kernel);
finish;
* Varying partly linear model estimation*
start S_metric(_y,u0,u10,_n,x0,z0,vary_var0);

mean=u0[:,];
stderr=sqrt((u0-mean)[##,]);
h=2.34#stderr#_n##(-0.2);
I=i(_n);
do j=1 to _n;
k0=(u0-u10[j,])/h;
w0=(1/h)#kde(k0);
w=diag(w0);
hx=(u0-u10[j,])/h;
hhx=hx#x0;
du=x0||hhx;
l=(x0[j,]||j(1,vary_var0,0))*inv(t(du)*w*du)*t(du)*w;
s=s//l[1,];
end;
return(s);
finish S_metric;
*define the estimation bhat*
start Para(_y,u0,u10,_n,x0,z0,vary_var);
I=i(_n);
S=S_metric(_y,u0,u10,_n,x0,z0,vary_var);
*linear part estimation*
bhat=inv(t(z0)t(I-s)(I-s)z0)*t(z0)*t(I-s)(I-s)*_y;
return(bhat);
finish Para;
Define the varying-coefficient named E*
start vary_para(_y,u0,u10,_n,x0,z0,vary_var);
bhat=para(_y,u0,u10,_n,x0,z0,vary_var);
mean=u0[:,];
stderr=sqrt((u0-mean)[##,]);
h=2.34#stderr#_n##(-0.2);
do i=1 to _n;
k1=(u0-u10[i,])/h;
w0=(1/h)#kde(k1);
w=diag(w0);
hx=(u0-u10[i,])/h;
hhx=hx#x0;
du=x0||hhx;
ahat=inv(t(du)w*du)*t(du)*w(_y-z0*bhat);
* 按照有无截距项和可变系数个数决定并修改*
if i=1 then e=ahat[1,]||ahat[2,]||ahat[3,];
else e=e//(ahat[1,]||ahat[2,]||ahat[3,]);
end;
return(e);
finish vary_para;
*Get resid of vary_semi estimation
start resid_boot(_y,u0,u10,_n,x0,z0,vary_var);
bhat=para(_y,u0,u10,_n,x0,z0,vary_var);
S=S_metric(_y,u0,u10,_n,x0,z0,vary_var);
resid_boot=_y-s*(_y-z0*bhat)-z0*bhat;
return(resid_boot);
finish resid_boot;
*Get the PLR statistic*
******Hypothsis testing for Parametric part,Statistic:PLR
firstly,Estimate Parametric part********
* where vary_var is not the same only*
start test_H(_y,u0,u10,_n,x0,z0,vary_var0,vary_var1);
bhat=para(_y,u0,u10,_n,x0,z0,vary_var0);
e=vary_para(_y,u0,u10,_n,x0,z0,vary_var0);
S=S_metric(_y,u0,u10,_n,x0,z0,vary_var0);
m_1=s*(_y-z0*bhat);
raw_rss1=_y-m_1-z0*bhat;
rss_1=raw_rss1[##,];
dim_z=ncol(z0);
do p=1 to dim_z;
* Take metric Z with pth observation deleted *;
nobs_bhat=nrow(bhat);
if p=1 then do;
_BHATdel_=bhat[2:nobs_bhat,]; ***
Zdel=z0[,2:ncol(z0)];
end;
else if p=nobs_bhat then do;
* _BHATdel_=bhat[1:(nobs_bhat-1),];*
Zdel=z0[,1:ncol(z0)-1];
end;
else do;
8*_BHATdel_=bhat[1:(p-1),]//bhat[(p+1):nobs_bhat,];**
Zdel=z0[,1:p-1]||z0[,(p+1):ncol(z0)];
end;
bhat0=para(y,u0,u10,_n,x0,_zdel,vary_var0);
e0=vary_para(y,u0,u10,_n,x0,_zdel,vary_var0);
S0=S_metric(y,u0,u10,_n,x0,_zdel,vary_var0);
m0=s0*(_y-_Zdel_*bhat0);
raw_rss0=_y-m0-_Zdel_*bhat0;
rss0=raw_rss0[##,];
delt=(1/_n)#rss0;
del_fin=del_fin||delt;
PLR=(_n/2)#log(rss0/rss_1);
PLR_PAR=PlR_PAR||PLR;
end;
*Get the GLR statistic***
*Hypothsis testing for varying Parametric part,Statistic:GLR*****
* create bhat and varying para******
dim_vary=ncol(x0)-1;
e_need=e[,2:dim_vary+1];*
do p=1 to dim_vary;
* Take metric X with pth observation deleted and append to metric Z*;
nobs=ncol(X0)-1;
if p=1 then do;
* _varydel_=e_need[,2:nobs];*
Xdel=X0[,1]||X0[,3:nobs+1];
*vary_X=_varydel_#_Xdel_[,+];***
_Zagu=X0[,2]||z0;
end;
else if p=nobs then do;
**_varydel_=e_need[,1:(nobs-1)];*
Xdel=X0[,1]||X0[,2:nobs];
*vary_X=_varydel_#_Xdel_[,+];**
_Zagu=X0[,nobs+1]||Z0;
end;
else do;
_varydel_=e_need[,1:(p-1)]||e_need[,(p+1):nobs];
Xdel=X0[,1]||X0[,2:p]||X0[,p+2:nobs+1];
vary_X=_varydel_#_Xdel_[,+];*
_Zagu=X0[,p+1]||Z0;
end;
bhat00=para(y,u0,u10,_n,_Xdel,_zagu,vary_var1);
e00=vary_para(y,u0,u10,_n,_Xdel,_zagu,vary_var1);
S00=S_metric(y,u0,u10,_n,_Xdel,_zagu,vary_var1);
m00=s00*(_y-_Zagu*bhat00);
raw_rss00=_y-m00-_Zagu*bhat00;
rss00=raw_rss00[##,];
delt=(1/_n)#rss00;
GLR=(_n/2)#log(rss00/rss_1);
GLR_NP=GLR_NP||GLR;
end;
* bb=bhat//j(nrow(e)-nrow(bhat),1,0);*
return(PLR_PAR||GLR_NP);
finish Test_H;
*************************************************************;
* Module BOOTSTRAP For Estimation of Null Distribution *;
*************************************************************;
start bootstrap(Y,_X,_Z,_U,_U1,_n,nboot,vary_var0,vary_var1);
resid_boot=resid_boot(y,_U,_U1,_n,_x,_z,vary_var0);
* Standardized observations (iid) **
* Generate bootstrap samples and the corresponding estimates *
do i=1 to nboot;
nobs=nrow(resid_boot);
noobs=2#nobs;
* Generate a bootstrap sample *;
do k=1 to nobs;
obs=int(_nobs_*ranuni(0))+1;
if k=1 then residboot=resid_boot[obs,];
else residboot=residboot//resid_boot[obs,];
end;
%let pi=3.1415926;
Yboot=3#sin(6#&pi#_u)#_x[,1]+2#sin(2#&pi#_u)#_x[,2]+2#_z[,1]+_z[,2]+0.5#_z[,3]+_residboot;
* Estimate parameters PLR and GLR*;
GLR=test_H(_Yboot,_U,_U1,_n,_x,_z,vary_var0,vary_var1);
get p_value
mean0=_u[:,];
stderr0=sqrt((_u-mean0)[##,]);
h0=2.34#stderr0#_n##(-0.2);
%let rk=2.1153;
%let ck=0.45;
glr_t=&rk#glr;
degree=&ck/h0;
p1=1-probchi(glr_t[,1],degree);
p2=1-probchi(glr_t[,2],degree);
p=p1||p2;
if i=1 then p_value=p;
else p_value=p_value//p;
* Store location estimates with sample identification *
if i=1 then GLR_boot=GLR;
else GlR_boot=GLR_boot//GLR;
end;
return(p_value);
finish bootstrap;
GLR_ORI=test_H(y,u,u1,n,x,z,3,2);
GLR_boot=bootstrap(y,x,z,u,u1,n,50,3,2);

GLR=GLR_ORI//GLR_boot;
do i=1 to nrow(glr);
if i>1 then do;
GLR_P=glr[i,]-glr[1,];
GLR_dif=GLR_dif//GLR_P;
end;
end;
glrpp=ranktie(glr);
* get solid parametric z1 z2 z3 and varying-coefficient x1 x2*
bhat_ori=Para(y,u,u1,n,x,z,3);
vary_ori=vary_para(y,u,u1,n,x,z,3);
create bhat_ori var{Bhat};
append from bhat_ori;
create vary_ori var{Inter_para para_x1 para_x2};
append from vary_ori;
create GLR_boot var{PLR_Z1 PLR_Z2 PLR_Z3 GLR_X1 GLR_X2};
append from GLR_boot;
create GLR_boot from glr_boot;
append from GLR_boot;
create GLR_ORI var{PLR_Z1 PLR_Z2 PLR_Z3 GLR_X1 GLR_X2};
append from GLR_ORI;
create glr_all from glr;
append from glr;
create glr_dif from GLR_dif;
append from GLR_dif;
*create p_value_glr from p_value;
*append from p_value;
quit;
data p_value;
set glr_dif;
array dif{*} col1-col5;
do i=1 to dim(dif);
if dif(i)<0 then dif(i)=0;
else dif(i)=1;
end;
run;
%macro varlist(type=NQ);
%do i=1 %to &nobs;
%if %upcase(&type)=Q %then sum(col&i=1)/count(col&i) as p_&i;
%else col&i;
%if &i ne &nobs %then %str(,);
%end;
%mend;
%let nobs=5;
proc sql;
create table p_cri as
select %varlist(type=Q)
from p_value;
quit;

非参数回归和相关统计检验相关推荐

在线进行 PCoA 分析和相关统计检验
PCoA分析一文学会PCA/PCoA相关统计检验(PERMANOVA)和可视化详细论述了PERMANOVA 检验(也包括最基本的方差检验基础),PERMANOVA检验的问题,并提供了代码生成 PC ...
一文学会PCA/PCoA相关统计检验(PERMANOVA)和可视化
方差分析中的元和因素试验中要考察的指标称为试验指标,影响试验指标的条件称为因素,因素所处的状态称为水平 (通常用于3个或更多水平时:如果只有2个水平考虑T-test);若试验中只有一个因素改变则称为 ...
对层级聚类树进行模块分割，定位基因在哪个模块中
拷贝数据到 ImageGP (http://www.ehbio.com/Cloud_Platform/front/#/analysis?page=b%27Ng%3D%3D%27),并设置参数. ``` ...
iMeta高被引论文｜陈同/刘永鑫等高颜值绘图网站imageGP被引500次(截止22/12/13)
https://doi.org/10.1002/imt2.5 2022年2月21日,iMeta 期刊在线发表了中国中医科学院黄璐琦院士.陈同副研究员和中国科学院遗传与发育生物学研究所刘永鑫高级工程师合 ...
iMeta：高颜值绘图网站imageGP+视频教程合集，已被引360次(220625更新)
https://doi.org/10.1002/imt2.5 2022年2月21日,iMeta 期刊在线发表了中国中医科学院黄璐琦院士.陈同副研究员和中国科学院遗传与发育生物学研究所刘永鑫高级工程师合 ...
你想要的宏基因组-微生物组知识全在这(2023.01)
欢迎点击上方蓝色"宏基因组"关注我们! 宏基因组/微生物组是当今世界科研最热门的研究领域之一,为加强宏基因组学技术和成果交流传播,推动全球华人微生物组领域发展,中科院青年科研人员创 ...
轻松在线差异基因/物种分析和可视化
7.1 limma用于转录组reads-count数据差异基因分析做转录组差异分析,我们通常是基于reads-count数据,这个数据可以来源于Salmon的结果.STAR的结果或从网上下载的数据. ...
fNIRS研究行文指南
近40年来,功能近红外成像技术(fNIRS)越来越多地应用于神经科学研究,利用各种实验范式研究不同人群.随着研究方法的快速发展与多样化,研究方法的不同给解释.重复研究带来很大困难.因此,fNIRS学会 ...
你想要的宏基因组-微生物组知识全在这(2022.12)
欢迎点击上方蓝色"宏基因组"关注我们! 宏基因组/微生物组是当今世界科研最热门的研究领域之一,为加强宏基因组学技术和成果交流传播,推动全球华人微生物组领域发展,中科院青年科研人员创 ...

非参数回归和相关统计检验

非参数回归和相关统计检验相关推荐

最新文章

热门文章