negbin.src

来自「没有说明」· SRC 代码 · 共 285 行

/*
** negbin.src - Negative Binomial Regression Model
**              (Truncation-at-zero and Variance Functions Optional)
**
** (C) Copyright 1988-1995  Aptech Systems, Inc.
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APTECH
** SYSTEMS, INC.    This File Header must accompany all files using
** any portion, in whole or in part, of this Source Code.   In
** addition, the right to create such files is strictly limited by
** Section 2.A. of the GAUSS Applications License Agreement
** accompanying this Software Product.
**
** If you wish to distribute any portion of the proprietary Source
** Code, in whole or in part, you must first obtain written
** permission from Aptech Systems.
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**  FORMAT:     { bg,vc,llik } = negbin(dataset,dep,ind1,ind2);
**
**  INPUT:
**      dataset = name of Gauss dataset or name of matrix in memory
**      dep     = dependent variable name or column number
**      ind1     = vector of independent variable names or column numbers
**      ind2    = 0 for negbin (scalar dispersion parameter), vector of var
**                names for variance function
**
**  OUTPUT:
**      bg    = vector of effect parameters that maximize the likelihood
**             on top of parameter(s) corresponding to ind2.
**             PARAMETERIZATION: bg=b|g;
**                 E(Y) = exp(ind1*b)
**                 V(Y) = E(Y)*(1+exp(g))       for ind2==0;
**                 V(Y) = E(Y)*(1+exp(ind2*g))  for ind2==vector
**      vc   = variance-covariance matrix of b
**      llik = value of the log-likelihood at the maximum
**
**  GLOBALS:
**    _cn_Inference    = MAXLIK for maximum likelihood estimates
**                    = BOOT for bootstrapped estimates
**                    = PROFILE for likelihood profile and profile t traces
**
**       _cn_ZeroTruncate  1 = regular model, 0 = truncated-at-zero model,
**                         (default 1)
**
**       _cn_Fix        0 = do nothing extra (default)
**                 symbol = include this variable from a dataset,
**                            constraining its coefficient to 1.0.
**                integer = include log of this column of input matrix,
**                            constraining its coefficient to 1.0.
**
**      _cn_Start   choose method of calculating starting values.
**                     0 = LS (default),
**                     1 = set to vector stored in _cn_StartValues,
**                     2 = rndu-0.5,
**                     3 = zeros, or set to vector
**
**      _cn_Dispersion   starting value for scalar dispersion parameter.
**                 Default = 3.
**
**      __output    1  =  print output to screen (default),
**                  0  =  do not print to screen
**
**  OTHER GLOBALS:
**      see MAXLIK.
**
**  EXAMPLE 1:                         EXAMPLE 2:
**    let dep = wars;                    let dep = wars;
**    let ind1 = age party unem;         let ind1 = age party unem;
**    dataset = "\\gauss\\prg\\sample";  let ind2 = coups us; @ var function @
**    call negbin(dataset,dep,ind1,0);   dataset = "sample";
**                                         call negbin(dataset,dep,ind1,ind2);
**
**  REFERENCE (for the basic negative binomial model):
**      Gary King. 1989. "Variance Specification in Event Count
**      Models:  From Restrictive Assumptions to a Generalized Estimator,"
**      AMERICAN JOURNAL OF POLITICAL SCIENCE, 33, 2 (forthcoming, August).
**
**  REFERENCE (negative binomial model with truncation and variance functions):
**      Gary King. 1989. "Event Count Models for International Relations:
**      Generalizations and Applications," INTERNATIONAL STUDIES QUARTERLY.
**      (forthcoming, June).
*/

#include count.ext
#include gauss.ext
#include maxlik.ext

proc _cn_svnb(dataset,dep,ind1,ind2);
    local res,b0,b1,b,pars;
    if ind1==0;
        pars = 1;
    else;
        pars = rows(ind1)+1;
    endif;
    If ind2==0;
        pars = pars+1;
    else;
        pars = pars+rows(ind2)+1;
    endif;
    if _cn_Start==0;
        if ind1==0;
            b0 = 0;
        else;
            b0 = lols(dataset,dep,ind1);
        endif;
        if ind2$/=0;        /* if variance function */
            b1 = lols(dataset,dep,ind2);
            res = b0|b1;
        else;       /* if scalar ancillary param */
            res = b0|_cn_Dispersion;
        endif;
    elseif _cn_Start==1;
        b = _cn_StartValues;
        res = b;
        if rows(b)/=pars;
            print "b is the wrong size for _cn_Start\g";
            end;
        endif;
    elseif _cn_Start==2;
        res = rndu(pars,1)-0.5;
    elseif _cn_Start==3;
        res = zeros(pars,1);
    else;
        res = _cn_Start;
        if rows(res)/=pars;
            errorlog "rows(_cn_Start) is wrong.\g";
            end;
        endif;
    endif;
    retp(res);
endp;

proc _cn_linb(b,dta);
    local res,xb1,l,t,lth,lthy,b0,b1,n,y,x0,x1,cx0,cx1;

    y = dta[.,1];
    n = rows(y);
    x0 = ones(n,1);
    if _cn_c1/=0;
        x0 = x0~dta[.,_cn_c1];
    endif;
    cx0 = cols(x0);
    x1 = ones(n,1);
    if _cn_c2/=0;
        x1 = x1~dta[.,_cn_c2];
    endif;
    cx1 = cols(x1);
    if _cn_ZeroTruncate==0;
        if sumc(y.==0)>0;
            print "Model not admissable.  Use _cn_ZeroTruncate==0 only if y"\
                "is truncated so that zeros do not appear in the data set.\g"\
                "";
            end;
        endif;
    endif;
    b0 = b[1:cx0];
    b1 = b[cx0+1:cx0+cx1];

    if _cn_Fix==0;
        l = exp(x0*b0);
    else;
        l = exp(x0*b0).*dta[.,cols(dta)];
    endif;
    xb1 = x1*b1;
    t = exp(xb1);
    lth = l./t;
    lthy = lth+y;
    res = _cn_lng(lthy)-_cn_lng(lth)+(y.*xb1)-lthy.*ln(1+t);
    if _cn_ZeroTruncate==0;
        res = res-ln(1-((1+t).^(-lth)));
    endif;
    retp(res);
endp;

proc 3 = negbin(dataset,dep,ind1,ind2);

    local vars,b,logl,vc,st,ret,nv;

    clearg _cn_c1,_cn_c2;
    _max_CovPar = 3;
    _cn_fn = dataset;       /* never used, should delete */
    if dep$==0;
        errorlog "DEP must = variable name or number";
        end;
    endif;
    if (type(dataset)==13) and (type(_cn_Fix)==13);
        _cn_Censor = indcv(_cn_Fix,getname(dataset));
    endif;
    if ((type(dataset)/=13) and (_cn_Fix>cols(dataset)));
        errorlog "If dataset=matrix, _cn_Fix must= 0 or a col of dataset\g";
        end;
    endif;
    if ((type(dataset)/=13) and ((maxc(ind1)>cols(dataset)) or
        (dep>cols(dataset))) );
        errorlog "If DATASET=matrix, DEP and ind1 must be column numbers of"\
            " the input matrix.\g";
        end;
    endif;
    vars = dep;
    if ind1==0;
        _cn_c1 = 0;
    else;
        _cn_c1 = seqa(2,1,rows(ind1));
        vars = vars|ind1;
    endif;
    if ind2==0;
        _cn_c2 = 0;
    else;
        _cn_c2 = seqa(rows(vars)+1,1,rows(ind2));
        vars = vars|ind2;
    endif;
    if _cn_Fix/=0;
        vars = vars|_cn_Fix;
    endif;
    st = _cn_svnb(dataset,dep,ind1,ind2);
    if __title $== "";
        if _cn_ZeroTruncate==0;
            __title = "Truncated ";
        else;
            __title = "";
        endif;
        __title = __title $+ "Negative Binomial Regression Model";
    endif;
    local infm,inf0,lcInf;
    infm = { MAXLIK, BOOT };
    inf0 = { 1, 2 };
    LcInf = _ml_check(_cn_Inference,1,infm,inf0,1);
    if LcInf == 1;
        { b,logl,vars,vc,ret } = maxlik(dataset,vars,&_cn_linb,st);
    elseif LcInf == 2;
        { b,logl,vars,vc,ret } = maxboot(dataset,vars,&_cn_linb,st);
    endif;
    if ret /= 0;
        errorlog "ERROR: Model estimation failed.";
        end;
    endif;

    vars = "beta0";

    if type(dataset)==13 and cols(dataset) >= 2;
        nv = getname(dataset);
        nv = rows(nv);
    else;
        nv = 1e+15;
    endif;

    if ind1/=0;
        if round(ind1) == ind1 and ind1 >= 1 and ind1 < 131072;
            if ind1 > nv;
                errorlog "ERROR: index of variable out of range: " $+
                    ftos(ind1,"%*.*lf",1,0);
                end;
            endif;
            vars = vars|
            ((0 $+ "Col." $+ zeros(rows(ind1),1))$+_cn_ftosm(ind1,2));
        else;
            vars = vars|ind1;
        endif;
    endif;
    vars = vars|"gamma0";
    if ind2/=0;
        if round(ind2) == ind2 and ind2 >= 1 and ind2 < 131072;
            if ind2 > nv;
                errorlog "ERROR: index of variable out of range: " $+
                    ftos(ind2,"%*.*lf",1,0);
                end;
            endif;
            vars = vars|
            ((0 $+ "Col." $+ zeros(rows(ind2),1))$+_cn_ftosm(ind2,2));
        else;
            vars = vars|ind2;
        endif;
    endif;

    _cn_vr = vars;
    _cn_dp = dep;
    ndpclex;
    retp(b,vc,logl*_max_NumObs);
endp;