logit.src

来自「没有说明」· SRC 代码 · 共 1,297 行 · 第 1/3 页

            ftos(-2*llrest,"*.*lf ",10,4));
        print ("  Percent Correctly Predicted:            " $+
            ftos(success,"*.*lf ",10,4));
        if _qrfit;
            print ("  Madalla's pseudo R-square:              " $+
                ftos(madr2,"*.*lf ",10,4));
            print ("  McFadden's pseudo R-square:             " $+
                ftos(mcfr2,"*.*lf ",10,4));
            print ("  Cragg and Uhler's pseudo R-square:      " $+
                ftos(cur2,"*.*lf ",10,4));
            print "\nOBSERVED AND PREDICTED OUTCOMES\n";
            dash = "-------------------------------------------------------"\
                "----------------------";
            spc = "                                                        "\
                "                     ";
            spc = strsect(spc,1,(4.5*ncat-4));
            print "           | " spc "Predicted";
            dash = strsect(dash,1,22+9*ncat);
            print "  Observed |";;
            call formatcv("*.*lf "~8~8);
            call formatnv("*.*lf "~8~0);
            call printfmt(_qrtmp'~"Total",0);
            print;
            print ("  " $+ dash);
            ii = 1;
            call formatcv("*.*lf "~10~8);
            do until ii>ncat;
                call printfmt(_qrtmp[ii],0);
                print "|";;
                call printfmt(obspred[ii,.],1);
                print;
                ii = ii+1;
            endo;
            print;
            print ("  " $+ dash);
            print "     Total |";;
            call printfmt(sumc(obspred)',1);
            print;
        endif;
        print;

    endif;
    if _qrpred;
        call formatcv("*.*lf "~10~8);
        nout = rowsf(fpred);
        fpred = close(fpred);
        if fpred .eq 0 and __output;
            print;
            print "PREDICTED VALUES SUCCESSFULLY WRITTEN TO DISK:";
            print;
            print ("  The file " $+ _qrpredn $+ " was written with " $+
                ftos(nout,"*.*lf ",1,0) $+ " cases.");
            print;
            print "  The following variables are in the file:";
            print;
            print ("    Prob Y=i for i=1 to " $+ ftos(ncat,"*.*lf ",1,0) $+
                ":   ");;
            call printfmt(prednm[1:ncat,1]',0);
            print;
            print;
            print ("    Dependent variable (Y):  " $+ ylbl);
        endif;
    endif;

    goto done;

pred:

    /* Compute predicted probabilities */

    fl = exp(x*bmat);
    ypred = fl./(1+sumc(fl'));      /* prob for all but last category  */
    g = ypred~(ones(rows(x),1) - sumc(ypred'));
    if vtype;
        ydum = ydum~(ynum.==_qrycat[ncat]);
    else;
        ydum = ydum~(ynum.$==_qrycat[ncat]);
    endif;
    ll = ll+sumc(wt.*sumc((ydum.*ln(g))'));

/* compute obs/predicted table */

    ii = ((vec((seqa(1,1,ncat)*ones(1,ncat))')*10)+reshape(seqa(1,1,ncat),
        ncat*ncat,1));
    obspred = obspred+countwts(10*maxindc(ydum') + maxindc(g'),ii,wt);
    return;

done:

    call formatcv(oldcfmt);
    call formatnv(oldnfmt);
    ndpclex;
    lrx2 = lrx2|-2*llfull|-2*llrest|success;
    retp(ylbl|xlbl,bprt,vcml,ndtran,pctall,meanx,sd,lrx2,df,tol);

errout:

    pop err;
    cls;
    if fin > 0;
        fin = close(fin);
    endif;
    call formatcv(oldcfmt);
    call formatnv(oldnfmt);
    ndpclex;
    if not trapchk(1);
        errorlog errmsg;
        end;
    endif;
    retp(0,err,0,0,0,0,0,0,0,99999);

endp;
/* End of LOGIT*/

proc(10) = logitprt(vnam,b,vc,nused,pct,mn,sd,lrx2,df,tol);
    local mask,i,fmt,fmt2,fmt3,iv,ij,ind,p,xx,il,ik,im,in,unstd,stand,vv,
        pp,ivv,jj,iw,kk,zp,ppv,tt,pp1,uu,ww,n, l1,l1end,ddf, omat, i1,l2,
        fmt4,ncat,catnm,vnames,dum,lx,j,ii,fm,yoff,x, mxt,maxx,minx,xtick,
        axis,wrk, strng,ioff,irow,icol,oldnfmt;
    if _qrev;
        b = -b;
    endif;
    msym "-----";
    dum = 0;
    ncat = 1+rows(b)/(rows(vnam)-1);
    if rows(_qrycat) /= ncat;
        _qrycat = seqa(1,1,ncat);
    endif;
    if __vtype[1] == 1 or vartype(vnam[1]);
        _qrycat = ftocv(_qrycat,1,0);
    endif;
    if rows(_qrtmp) /= ncat;
        catnm = 0 $+ "Cat_" $+ _qrycat;
    else;
        catnm = _qrtmp;
    endif;
    vnames = trimr(vnam,1,0);
    vnames[1] = "CONSTANT";
    fmt4 = { "*.*s" 8 8, "#*.*lG" 12 4, "*.*lf" 10 2, "*.*lf" 10 2, "*.*lf"
        10 3, "*.*lf" 10 2, "*.*lf" 10 3 };
    if (rows(dum) == 1) and (dum == 0);
        dum = ones(rows(vnames)-1,1);
    endif;
    print "LOGIT EFFECT COEFFICIENTS";
    print;
    print "NOTE:  Percent change coefficients indicate the percent change i"\
        "n the odds";
    print "       for a change in the independent variable.";
    print;
    print "Categories are designated as:";
    print;
    i = 1;
    j = 0;
    do until i>ncat;
        if __vtype[1] == 1;
            print ("  " $+ _qrycat[i] $+ " - " $+ catnm[i,1] $+ "    ");
        else;
            print ("  " $+ ftocv(i,1,0) $+ " - " $+ catnm[i,1] $+ "    ");
        endif;
        if j==2;
            j = 0;
            print;
        else;
            j = j+1;
        endif;
        i = i+1;
    endo;
    print;
/* prepare table output */
    mask = 0~ones(1,ncat);
    i = 1;
    let fmt[1,3] = "*.*ld" 10 8;
    fmt2 = fmt;
    fmt3 = fmt;
    do until i>ncat;
        fmt = fmt|("*.*lf "~10~4);
        fmt2 = fmt2|("*.*lf "~10~2);
        fmt3 = fmt3|("*.*lf "~10~3);
        i = i+1;
    endo;
    vv = sqrt(diag(vc));    /* get starting standard error values */
    iv = rows(sd);
    if _qrcon==0;           /* skip over constant or not */
        ij = 2;
        ind = 1+(ncat-1);
    else;
        ij = 1;
        ind = 1;
    endif;
    do until ij>iv;
    /* Compute Logit Coefficients */
        p = b[ind:(ind+(ncat-2)),1];
        xx = zeros(ncat,(ncat-1));
        p = p|0;
        xx = xx~p;
        il = 1;
        do until il==ncat;
            ik = il+1;
            do until ik==ncat;
                xx[il,ik] = xx[il,ncat] - xx[ik,ncat];
                ik = ik+1;
            endo;
            il = il+1;
        endo;
        im = 1;
        do until im==ncat;
            in = im+1;
            do until in>ncat;
                xx[in,im] = -xx[im,in];
                in = in+1;
            endo;
            im = im+1;
        endo;
    /* standardized and unst effect coefficients */
        unstd = (exp(xx)-1)*100;
        stand = (exp(xx.*sd[ij,.])-1)*100;
    /* figuring standard errors */
        pp = vv[ind:ind+(ncat-2),1];
        pp = pp|0;
        zp = zeros(ncat,ncat-1);
        pp = zp~pp;
        ivv = 1;
        jj = 0;
        do until ivv==ncat-1;
            iw = ivv+1;
            kk = 1;
            do until iw>ncat-1;
                uu = ind+jj;
                ww = ind+jj+kk;
                pp[ivv,iw] = vc[uu,uu]+vc[ww,ww]-2*vc[uu,ww];
                pp[ivv,iw] = sqrt(pp[ivv,iw]);
                kk = kk+1;
                iw = iw+1;
            endo;
            jj = jj+1;
            ivv = ivv+1;
        endo;
        im = 1;
        do until im==ncat;
            in = im+1;
            do until in>ncat;
                pp[in,im] = pp[im,in];
                in = in+1;
            endo;
            im = im+1;
        endo;
    /* t statistics and probabilities */
        pp1 = diagrv(pp,ones(ncat,1));
        tt = xx./pp1;
        tt = diagrv(tt,zeros(ncat,1));
        ddf = nused-((ncat-1)*iv);
        ppv = 2*cdftc(abs(tt),ddf[1,1]);
    /* printing results */
        omat = ones(1,7);
        n = upper(vnames[ij]);
        l1 = 1;
        l1end = ncat;
        i1 = l1;
        do until i1 > l1end;
            l2 = 1;
            do until l2 > ncat;
                if i1 .ne l2;
                    omat = omat|(0 $+ _qrycat[i1] $+ " vs " $+ _qrycat[l2]
                        $+ ":  ")~xx[i1,l2]~ unstd[i1,l2]~stand[i1,l2] ~
                        pp1[i1,l2]~ tt[i1,l2] ~ppv[i1,l2];
                endif;
                l2 = l2+1;
            endo;
            i1 = i1+1;
        endo;
        let mask = 0 1 1 1 1 1 1;
        print "              Logit    Unstd %     Std %      Std";
        print ("" $+ n $+ chrs(32*ones(1,10-strlen(n))) $+ "     Coef     C"\
            "hange    Change     Error   t-value     p>|t|");
        print "------------------------------------------------------------"\
            "------------";
        omat = trimr(omat,1,0);
        if ij == 1;         /* constant */
            omat[.,3] = miss(1,1)*ones(rows(omat),1);
            omat[.,4] = miss(1,1)*ones(rows(omat),1);
        endif;
        call printfm(omat,mask',fmt4);
        print;
        ij = ij+1;
        ind = ind+(ncat-1);
    endo;

/* PLOT EFFECT COEFFICIENTS */
    if _qrplot;
        local _same, updown,myt,xwrk;
        print;
        print "PLOTS OF EFFECT COEFFICIENTS:";
        print;
        _same = 1;          /* 1 if same scale for std and ustd */
                            /* sd = sd[2:iv+1,1]; */
    /* y offsets for plots depending on # of categories */
        if ncat eq 2;
            let updown = 0 0;
        elseif ncat eq 3;
            let updown = 0 -1 0;
        elseif ncat eq 4;
            let updown = 0 -1 1 0;
        elseif ncat eq 5;
            let updown = 0 -1 1 -1 0;
        elseif ncat eq 6;
            let updown = 1 -1 0 1 -1 0;
        elseif ncat eq 7;
            let updown = 2 1 -1 0 1 -1 0;
        elseif ncat eq 8;
            let updown = 2 1 -1 0 1 -1 0 -2;
        elseif ncat eq 9;
            let updown = 2 1 -1 3 0 1 -1 0 -2;
        elseif ncat eq 10;
            let updown = 2 1 -1 3 0 1 -1 -3 0 -2;
        endif;
        lx = reshape(b,iv,ncat-1);
        lx = lx~zeros(rows(lx),1);
        lx = trimr(lx,1,0);
        sd = trimr(sd,1,0);
        vnames = trimr(vnames,1,0);
        iv = rows(vnames) + 1;
    /* Loop over unstandardized and standardized plot */
        ii = 1;     /* type of plot */
        do until ii>2;
            print;
            if ii == 1;     /* fm is coefficients to plot */
                print "Unstandardized % Change:";
                fm = exp(lx);
            else;
                print "Standardized % Change:";
                fm = exp(lx.*sd);
            endif;
            print;
            yoff = 0;
            i = 1;
            do until i > iv-1;
                x = (fm[i,1:ncat]|seqa(1,1,ncat)')';
                x = sortc(x,1)~updown;
                x = sortc(x,2);
                yoff = yoff|x[1:ncat,3];
                i = i+1;
            endo;
            yoff = yoff[2:rows(yoff),1];

/* PLOT LOG OF EFFECTS */
            mxt = ln(fm);
            myt = reshape(yoff,iv,ncat);
            if _same==1;    /* if same metric for std and ustd plot  */
                maxx = maxc(maxc(((lx.*sd)|lx)));
                minx = minc(minc(((lx.*sd)|lx)));
            else;
                maxx = maxc(maxc(mxt));
                minx = minc(minc(mxt));
            endif;
            xtick = ( seqa(0,1,7) * ((maxx-minx)/6) ) + minx;
            xwrk = round((mxt-minx)./((maxx-minx)/60))+1;
            i = 1;
            axis = "+---------+---------+---------+---------+---------+----"\
                "-----+";
            axis = "            " $+ axis;
            print "                                  Percent Change Scale  "\
                "     ";
            print "     ";;
            oldnfmt = formatnv("*.*lf"~10~2);
            call printfmt(100*(exp(xtick)-1)',1);
            print;
            do until i > iv-1;
                print axis;
                wrk = xwrk[i,.]|myt[i,.]|seqa(1,1,ncat)';
                wrk = wrk|(wrk[2,.].*1000+wrk[1,.]);
                wrk = sortc(wrk',4);
                strng = strsect(""$+vnames[i,1]$+":              ",1,12);
                ioff = minc(updown);
                irow = 1;
                do until ioff>maxc(updown);
                    icol = 1;
                    if wrk[irow,2] == ioff;
                        do until ((icol>61) or (irow>ncat));
                            if wrk[irow,1]==icol;
                                strng = strng $+ _qrycat[wrk[irow,3]];
                                irow = irow+1;
                            else;
                                strng = strng$+" ";
                            endif;
                            icol = icol+1;
                        endo;
                    endif;
                    print strng;
                    strng = "" $+"            ";
                    ioff = ioff+1;
                endo;
                i = i+1;
            endo;
            print axis;
            print "     ";;
            call printfmt(xtick',1);
            print;
            print "                              Logit Coefficients Scale";
            print;
            ii = ii+1;
        endo;
        call formatnv(oldnfmt);
    endif;
    msym ".";
    retp(vnam,b,vc,nused,pct,mn,sd,lrx2,df,tol);
endp;       /* LOGITPRT */

proc(0) = logiter(iter,st,en,lbl);
    local omat,fmt,mask,st0,i;
    print;
    print "ITERATION " ftos(iter,"*.*lf",1,0) ":";
    print;
    print "  Parameter    Start Value        End Value    % Change";
    print "  -----------------------------------------------------";
    mask = 0~1~1~1;
    let fmt[4,3] = "*.*lf " 10 8 "*.*lf " 16 8 "*.*lf " 16 8 "*.*lf " 11 5;
    i = 0;
    st0 = st;
    do until i == rows(st);
        i = i+1;
        if st[i]/=0;
            st0[i] = 100*(en[i]-st[i])/st[i];
        else;
            st0[i] = -999;
        endif;
    endo;
    omat = lbl~st~en~st0;
    call printfm(omat,mask,fmt);
endp;       /* logiter */