HPGENSELECT Restrict语句

Question

我正在尝试理解如何在PROC HPGENSELECT中使用RESTRICT语句来拟合多项式logit模型。

我有一个分类变量evt_type，它有3个值(0,1,2)。我有两个连续的回归变量x1和x2。我希望将evt_type=2的x1参数限制为0。

不受约束的拟合：

proc hpgenselect data=to_reg;
   class evt_type;
   model evt_type(ref='0') = x1 x2 / link=glogit;
run;

给了我：

                                     Parameter Estimates

                                                      Standard
   Parameter    evt_type    DF        Estimate           Error    Chi-Square    Pr > ChiSq

   Intercept    1            1       -0.092251        0.109479        0.7100        0.3994
   Intercept    2            1       -0.181718        0.119293        2.3204        0.1277
   x1           1            1        0.875623        0.097909       79.9807        <.0001
   x1           2            1       -0.017942        0.110512        0.0264        0.8710
   x2           1            1       -0.220358        0.113552        3.7659        0.0523
   x2           2            1        2.902520        0.199583      211.4960        <.0001

根据文档尝试使用restrict，它看起来像是：

proc hpgenselect data=to_reg;
   class evt_type;
   model evt_type(ref='0') = x1 x2 / link=glogit;
   restrict x1 0, x1 1  = 0;
run;

然而，这似乎限制了evt_type=1上的参数。我做错了什么：

                                         Parameter Estimates

                                                      Standard
   Parameter    evt_type    DF        Estimate           Error    Chi-Square    Pr > ChiSq

   Intercept    1            1        0.028338        0.100646        0.0793        0.7783
   Intercept    2            1       -0.202028        0.117543        2.9541        0.0857
   x1           1            0               0               .         .             .
   x1           2            1       -0.423127        0.099020       18.2597        <.0001
   x2           1            1       -0.159225        0.104655        2.3147        0.1282
   x2           2            1        2.914689        0.199375      213.7191        <.0001

使用以下命令生成示例数据：

%macro generate_sample(seed);
    proc iml;
    n=1000;
    call randseed(&seed);
    cov = {1 .0, .0 1};
    x = randnormal(n,{0,0},cov);

    c0 = {0, 0};
    c1 = {1, 0};
    c2 = {0, 3};

    vProb = exp(x*c0) || exp(x*c1) || exp(x*c2);

    vProb = vProb/ vProb[,+];

    NumTrials = 1;
    events = J(n,3,0);

    do i=1 to n;
        prob = vProb[i,];
        events[i,] = RandMultinomial(1,1,prob);
    end;

    out = events || x;

    create events from out[colname={'e0','e1','e2','x1', 'x2'}];
    append from out;
    close events;

    out = vProb || x;

    create true_prob from out[colname={"P_0","P_1","P_2",'x1', 'x2'}];
    append from out;
    close true_prob;
    quit;

    data to_reg;
    format evt_type $1.;
    set events;
    array e[3] e0-e2;

    do i=0 to 2;
        if e[i+1] then do;
            evt_type = put(i,1.);
            leave;
        end;
    end;
    drop i e0-e2;
    run;
%mend;
%generate_sample(2);

Answer 1

与SAS技术支持一起工作。似乎问题出在ref='0'语句上。将模型顺序从ref='0'更改为descending似乎可以让事情变得更好。

这将失败：

proc hpgenselect data=to_reg ;
   model evt_type(ref='0')= x1 x2 / link=glogit dist=multinomial;
   restrict x1 1 , x1 0 =0;
   restrict x2 0 , x2 1 =0;
run;

在此过程中：

proc hpgenselect data=to_reg ;
   model evt_type(descending)= x1 x2 / link=glogit dist=multinomial;
   restrict x1 1 , x1 0 =0;
   restrict x2 0 , x2 1 =0;
run;

Answer 2

这似乎与使用class语句时SAS内部排序的方式有关。如果您设置了类选项order = data，它将适当地应用限制。

proc hpgenselect data=to_reg;
   class evt_type / order=data;
   model evt_type(ref='0') = x1 x2 / link=glogit;
   restrict x1 0, x1 1 = 0;
run;