psnreg.src

来自「没有说明」· SRC 代码 · 共 718 行 · 第 1/2 页

/*
** psnreg.src - Poisson Regression
**
**
** (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.
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**> psnreg
**
**
**  Purpose:  To estimate the poisson regression model using a GAUSS data set.
**            The details of the model are presented in the remarks section.
**
**  Format:   { vnames,b,vc,ndtran,meanx,sdx,fit,df,tol }
**              = psnreg(dataset,depvar,indvars);
**
**  Input:    dataset   string, name of data file
**
**            depvar    string, name of dependent variable
**                         - or -
**                      scalar, index of dependent variable.
**
**                      The value of depvar will be truncated before analysis.
**                      Thus, 1.4 is treated as category 1.
**
**            indvars    Kx1 character vector, names of independent variables.
**                          - or -
**                       Kx1 numeric vector, indices of independent variables.
**
**                       The program adds one variable for the constant term.
**
**            Defaults are provided for the following global input
**            variables.  They can be ignored unless you need control
**            over the other options provided by this procedure.
**
**            WARNING:  If you change the defaults in a command file, the new
**            values will apply in the next program you run using PSNREG unless
**            you change them back.  This can be done by running the files
**            gauss.dec and quantal.dec that contains the declare statements
**            for the globals.  If you want to make permanent changes to the
**            globals, they should be changed in gauss.dec and quantal.dec.
**
**            __output    global scalar, default 1.
**
**                        if 1, sends results to the output device (including
**                        the screen).
**
**                        if 0, no information is sent to output.
**
**            __range     2*1 vector.  The range of record in data
**                        set used for analysis.  The first element is the
**                        starting row index, the second element is the
**                        endding row index.
**
**                        Default is the whole dataset.
**
**            __row    global scalar, default 0.
**
**                     if 0, the number of rows to read per iteration of the
**                     read loop is calculated by the program.
**
**                     if not 0, the specified number of rows will be read.
**
**          __miss    global scalar, default 1.
**
**                      if 0, there are no missing values (fastest).
**
**                      if 1, do listwise deletion, drop an observation
**                      if there are any missing values among the independent
**                      and dependent variables.
**
**          __tol    global scalar controlling the iterations.  __tol
**                   indicates the maximum difference between estimates
**                   of the coefficients in two adjacent iterations.
**
**          _qriter    global scalar, default 0.
**
**                     if 0, do not print information on iterations.
**
**                     if 1, send detailed information on iterations to
**                     the screen but not to the output device.
**
**                     if 2, send detailed information on iterations to
**                     the output device.
**
**          _qrstat    global scalar, default 0.
**
**                     if 0, do not print descriptive statistics.
**
**                     if 1, print descriptive statistics.
**
**           _qrmiter   maximum number of iterations, default = 1000.
**
**
**  Output:
**
**          vnames    a (K+2)x1 character vector containing the names of
**                    the variables in the model.  The order is:
**                    depvar|"CONSTANT"|indvars.
**
**          b    an NPARM=*(K+1) vector of parameter estimates in
**               the order:  intercept|var1|var2|...varK.
**
**               If errors are encountered a message will be sent to the
**               error log.  Also, b will contain a scalar error code.  This
**               code appears as missing unless it is translated with
**               the command scalerr(b).  The codes are defined as:
**
**                       1  data file not found
**                       2  found undefined variables
**                       30  system singular
**                       31  too few nonmissing observations.
**                       72  one of the outcome categories has no cases
**                       73  an independent variable has no variation
**                       74  can't open file for predicted values
**                       75  out of disk space
**                       77  all cases were deleted
**
**          vc    NPARMxNPARM variance covariance matrix for the parameters
**                in b.
**
**          ndtran      2x1 vector of observations.  Element 1 contains
**                      number of cases read from dataset; element 2
**                      contains number of cases left after deletion of
**                      missing cases controlled by __miss, it is the
**                      number of cases used in the analysis.
**
**          meanx    the means based on nused cases of the independent
**                   variables in the order in indvars.
**
**          sdx    the standard deviations based on nused cases of the
**                 independent variables in the order in indvars.
**
**          fit    3x1 vector of goodness of fit measures.  Element 1 is
**                 the likelihood ratio chi-square assessing the overall
**                 fit of the model against a null model with only an
**                 intercept; element 2 is -2 times the log
**                 likelihood function evaluated at the estimated values;
**                 element 3 is -2 times the log likelihood function
**                 evaluated with the slopes fixed to zero.
**
**          df    the degrees of freedom associated with lrx2; this is
**                equal to the number of parameter less 1 for the intercept.
**
**          tol    the tolerance reached.  If convergence was obtained,
**                 tol must be less than __tol.
**
**  Remarks:
**
**         POISSON estimates the Poisson regression model for count data.
**         In this model the depedent variable Y(i) is the number of
**         occurrences of some event during a fixed interval of time.
**         The probability of Yi being some value yi is defined as:
**
**                          -Li   yi
**           Prob(Yi=yi) = e    Li   / yi!  for yi=0,1,...; i=1,...N.
**
**         This probability is influenced by exogenous variables through
**         the effect of the exogenous variables on Li (read: Lambda sub i)
**         through the equation:
**
**           Li = exp(Xi'B)
**
**         where Xi is a Kx1 vector of values of the exogenous variables for
**         observation i; B is a Kx1 vector of parameters.  It is assumed that
**         the first element of Xi is 1 and the first value of B is the
**         intercept.
**
**         For additional details, see:  A. Colin Cameron and Pravin K. Trivedi,
**         1986, "Econometric models based on count data."  Journal of Applied
**         Econometrics, 1, 29-53; or, William H. Greene, 1990, Econometric
**         Analysis, Macmillan, pages 707-709.
**
*/

#include gauss.ext
#include quantal.ext

proc (9) = psnreg(dataset,depvar,indvars);

    local fin,errmsg,nivar,nr,dta,y,x,minx,maxx,miny,maxy,xx,xy,nused,
        mask,omat,fmt,bstrt,bchng,iter,tm,tzu,w,vcml,bnew,sdml,tml,pml,tol,
        llrest,lrx2,df,lrxp,mcfr2,err,llfull,oldtrap,xlbl,ylbl,nm,readdisk,
        nd, ndtran,wx,sdy,meany,yy,meanlny,iswt, wlbl,windx,wt,meanx,f,sdx,
        wnused,olsst,xi,lam,nparm,bonly,wy,varindx,count,counter, lastobs,
        start,range,cv,ky;

    fin = -1;
/*
    if type(__weight) == 13;
        wtf = stof(__weight);
    else;
        wtf = __weight;
    endif;
    if wtf == 0;
        iswt = 0;
    else;
        wvar = wtf;
        iswt = 1;
        iv = iv|wvar;
    endif;
*/
    dataset = "" $+ dataset;
    open fin = ^dataset;
    if fin == -1;
        errmsg = "ERROR:  Can't open file " $+ dataset $+ ".";
        goto errout(error(1));
    endif;

    { nr,start,counter,lastobs } = _rngchk(dataset,__range);
    range = lastobs - start + 1;
    ndtran = zeros(2,1);
    ndtran[1] = range;
/*
    { nm, varindx } = indices(dataset, dv|iv);
    if scalerr(nm);
        goto errout(nm);
    endif;

    ylbl = nm[1];
    xlbl = "CONSTANT"|trimr(nm,1,0);
    if iswt;
        xlbl = trimr(xlbl,0,1);
    endif;
    nivar = rows(nm)-iswt;
*/
    { nm, varindx } = indices(dataset, depvar|indvars);
    if scalerr(nm);
        goto errout(nm);
    endif;
    ylbl = nm[1];
    xlbl = "CONSTANT"|trimr(nm,1,0);
    nivar = rows(nm);
    if __weight $== 0;
        wlbl = "";
        windx = 0;
        iswt = 0;
    else;           /* INCLUDES __WEIGHT IF DEFINED */
        { wlbl, windx } = indices(dataset, __weight);
        if scalerr(wlbl);
            goto errout(wlbl);
        endif;
        iswt = 1;
        nm = nm|wlbl;
        varindx = varindx|windx;
    endif;

    minx = (1E+50)*ones(nivar,1);
    miny = (1E+50);
    maxy = -miny;
    maxx = -minx;
    clear readdisk,nused;
    /* OLS Approximation and Descriptive Statistics */
    clear xx,xy,meanx,wnused,meany,yy,meanlny;
    nd = 0;
    call seekr(fin,start);
    count = counter;
    do while count < lastobs;
        readdisk = readdisk+1;
        dta = readr(fin,nr);
        count = count + rows(dta);
        if count > lastobs;
            dta = trimr(dta,0,count-lastobs);
        endif;
        if __miss == 1;
            dta = packr(dta[.,varindx]);
        else;
            dta = dta[.,varindx];
        endif;
        nd = rows(dta);
        if nd == 0;
            continue;       /* no observations returned */
        endif;
        ndtran[2] = ndtran[2] + nd;
        y = trunc(dta[.,1]);        /* truncate for poisson */
        x = ones(rows(dta),1)~trimr(dta',1,0)';
        if iswt;
            wt = dta[.,cols(dta)];
            x = trimr(x',0,1)';
        else;
            wt = ones(rows(dta),1);
        endif;
        wnused = wnused+sumc(wt);
        wy = y.*wt;
        meany = meany+sumc(wy);
        yy = yy + wy'*y;
        wy = ln(y+.5).*wt;
        meanlny = meanlny+sumc(wy);
        wx = x.*wt;
        xx = xx + wx'*x;
        xy = xy + wx'*ln(y+.5);
        meanx = meanx+sumc(wx);
        minx = minc((minx~minc(x))');
        miny = minc((miny~minc(y))');
        maxx = maxc((maxx~maxc(x))');
        maxy = maxc((maxy~maxc(y))');
    endo;
    nused = ndtran[2];
    if nused == 0;
        errmsg = "ERROR:  No observations left after deleting missing values.";
        goto errout(error(77));
    endif;
    readdisk = readdisk>1;
    clear dta;
    if not iswt;
        wnused = nused;
    endif;
    meanx = meanx/wnused;
    meany = meany/wnused;
    meanlny = meanlny/wnused;
    sdx = sqrt((diag(xx)./(wnused))-(meanx.*meanx));
    sdy = sqrt((diag(yy)./(wnused))-(meany.*meany));

    if __output;
        call header("PSNREG",dataset,_qr_ver);
        print;
        print "CASES PROCESSED BY PSNREG:";
        print;
        print "  " ftos(nused,"%-*.*lf",1,0) " cases were kept out of "
            ftos(range,"%-*.*lf",1,0) " in file.";
        if iswt;
            print;
            print "WEIGHTING IS IN EFFECT:  ";
            print;
            print ("  Computations have used the weight variable:  " $+ wlbl);
        endif;
        call formatcv("*.*lf"~10~8);
        call formatnv("*.*lf"~10~4);
        if _qrstat;
            print;
            print ("DESCRIPTIVE STATISTICS (N="$+ftos(ndtran[2],"*.*lf",1,
                0) $+"):");
            if iswt;
                print ("                       (N="$+ftos(wnused,"*.*lf",1,
                    0) $+" weighted):");
            endif;
            print;
            print "                Mean   Std Dev   Minimum   Maximum";
            mask = 0~1~1~1~1;
            omat = trimr(xlbl~meanx ~sdx~minx~maxx,1,0);
            omat = (ylbl~meany~sdy~miny~maxy)|omat;
            let fmt[5,3] = "*.*lf" 10 8 "*.*lf" 10 4 "*.*lf" 10 4 "*.*lf"
                10 4 "*.*lf" 10 4;
            call printfm(omat,mask,fmt);
        endif;