ordered.src

没有说明

    endif;
    ncon = ncat-1;          /* number of constant terms */
    nparm = ncon+nvar;      /* Total number of parameters */
    obspred = zeros(ncat*ncat,1);           /* table of obs and predicted
                                            :: outcomes
                                            */
    if readdisk;
        minx = (1e+50)*ones(nvar,1);
        maxx = -minx;
        pct = zeros(ncat,1);
        clear xx,xy,meanx;
        call seekr(fin,start);
        count = counter;
        do while count < lastobs;
            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;
            if scalmiss(dta);
                continue;
            endif;
            if iswt;
                wt = dta[.,cols(dta)];
                dta = trimr(dta',0,1)';
            else;
                wt = ones(rows(dta),1);
            endif;
            y = maxindc(((ones(rows(dta),1)*ycat .== trunc(dta[.,1])))');
            x = trimr(dta',1,0)';
            wx = x.*wt;
            cx = ones(rows(dta),1)~x;
            xx = xx + (cx.*wt)'*cx;
            xy = xy + (cx.*wt)'*y;
            meanx = meanx+sumc(wx);
            minx = minc((minx~minc(x))');
            maxx = maxc((maxx~maxc(x))');
            pct = pct + countwts(y,seqa(1,1,ncat),wt);
        endo;
    else;
        if iswt;
            wt = dta[.,cols(dta)];
            dta = trimr(dta',0,1)';
        else;
            wt = ones(rows(dta),1);
        endif;
        ydum = trunc(dta[.,1]);
        y = maxindc(((ones(rows(dta),1)*ycat .== ydum))');
        x = trimr(dta',1,0)';
        wx = x.*wt;
        cx = ones(rows(dta),1)~x;
        xx = (cx.*wt)'*cx;
        xy = (cx.*wt)'*y;
        meanx = sumc(wx);
        minx = minc(x);
        maxx = maxc(x);
        pct = countwts(y,seqa(1,1,ncat),wt);
    endif;

    nb = pct;
    pct = pct/wnused;       /* % in categories 1 to ncon */
    meanx = meanx/wnused;

    /* start values */

    if ((_qrstart eq 0) .and (rows(_qrstart) .eq 1));
        oldtrap = trapchk(0xffff);
        trap 1;
        bstrt = solpd(xy,xx);
        trap oldtrap,0xffff;
        if scalerr(bstrt);
            errmsg = "ERROR:  X'X matrix is singular.";
            goto errout(bstrt);
        endif;
        bstrt = trimr(bstrt,1,0);
        bm = bstrt'*meanx;

/* Estimate thresholds based on means of x and marginal percentages */
        let cd = -6 -1.7 -1.3 -1 -.8 -.6 -.5 -.4 -.25 -.1 0 .1 .25 .4 .5
            .6 .8 1 1.3 1.7 6;
/* approximate the cummulative distribution at 5% points */
        pl = cd[ceil(20*pct[1])];
/*        pu = cd[ceil(20*sumc(pct))]; */
        pu = cd[20];
        astrt = seqa(pl,(pu-pl)/(ncon-1),ncon)+bm';
        abstrt = astrt|bstrt;

    else;           /* user supplied start values */
        if rows(_qrstart) == (ncon+nvar);
            abstrt = _qrstart;
        else;
            errmsg = "ERROR:  Wrong number of start values in _qrstart.";
            goto errout(error(79));
        endif;
    endif;

    xx = diag(xx);
    sdx = sqrt((xx[2:rows(xx)]./wnused)-(meanx.*meanx));
    if __output;
        call header("ORDERED",dataset,_qr_ver);
        print;
        print "CASES PROCESSED BY THIS PROGRAM:";
        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;
        print;
    endif;
    call formatcv("*.*lf"~10~8);
    call formatnv("*.*lf"~10~4);

    if __output;
        print "DEPENDENT CATEORIES ARE DESIGNATED AS:";
        print;
        i1 = 1;
        i2 = 0;
        do until i1 > ncat;
            print ("  " $+ ftocv(_qrycat[i1],1,0) $+ " - " $+ _qrtmp[i1,1]
                $+ "    ");;
            if i2 == 3;
                i2 = 0;
                print;
            else;
                i2 = i2 + 1;
            endif;
            i1 = i1 + 1;
        endo;
        print;
    endif;

    if __output and _qrstat;
        print ("DISTRIBUTION AMONG ORDINAL GROUPS FOR VARIABLE " $+ ylbl $+
            ":");
        print;
        call printfmt(" "~_qrtmp',zeros(1,rows(_qrtmp)+1));
        print;
        print "PROPORTION";;
        call printfmt(pct',1);
        print;
        if minc(pct) < .00000001;
            errmsg = "ERROR:  An outcome category has no cases.";
            goto errout(error(72));
        endif;
        print;
        print ("DESCRIPTIVE STATISTICS (N=" $+ ftos(nused,"*.*lf",1,0) $+ ""\
            "):");
        if iswt;
            print ("                       (N="$+ftos(wnused,"*.*lf",1,0)
                $+ " weighted):");
        endif;
        print;
        let lbl[1,5] = "  " "Mean" "Std Dev" "Minimum" "Maximum";
        call printfmt(lbl,0);
        print;
        mask = 0~1~1~1~1;
        omat = xlbl~meanx~sdx~minx~maxx;
        let fmt[5,3] = "*.*lf" 10 8 "*.*lf" 10 4 "*.*lf" 10 4 "*.*lf" 10 4
            "*.*lf" 10 4;
        call printfm(omat,mask,fmt);
        print;
    endif;

    if maxc(minx .eq maxx) == 1;
        errmsg = "ERROR:  An independent variable has no variation.";
        goto errout(error(73));
    endif;

    if wnused < nvar+1;
        errmsg = "ERROR:  Too few nonmissing observations.";
        goto errout(error(31));
    endif;

    abold = abstrt;
/* Start iterations */

    abchng = 1;
    llchng = 1;
    iter = 1;

    clear llnew,absqz,cdfd,g;
    llold = 1000;
    ky = upper(chrs(key)) $== "C";
    do until abchng < __tol or iter > _qrmiter or ky;
        clear m,s,llwrk,sqz;
        if readdisk;
            call seekr(fin,start);
            count = counter;
            do while count < lastobs;
                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;
                if scalmiss(dta);
                    continue;
                endif;
                y = maxindc(((ones(rows(dta),1)*ycat .== trunc(dta[.,1])))');
                if iswt;
                    wt = dta[.,cols(dta)];
                    dta = trimr(dta',0,1)';
                else;
                    wt = ones(rows(dta),1);
                endif;
                x = trimr(dta',1,0)';
                n = rows(y);
                mask = reshape(seqa(1,1,ncon+2)',n,ncon+2);
                gosub grad;
                m = m+moment(g,0);          /* computes g'g */
                s = sumc(g)+s;
                if cdfd > 0;
                    llwrk = -sumc(wt.*ln(cdfd))+llwrk;
                else;
                    llwrk = -sumc(wt.*ln(cdfd-((cdfd.<=0).*(cdfd-.0001))))+
                        llwrk;
    /* fix up zero's */
                endif;
            endo;
        else;
            n = rows(y);
            mask = reshape(seqa(1,1,ncon+2)',n,ncon+2);
            gosub grad;
            m = moment(g,0);        /* computes g'g */
            s = sumc(g);
            if cdfd > 0;
                llwrk = -sumc(wt.*ln(cdfd));
            else;
                llwrk = -sumc(wt.*ln(cdfd-((cdfd.<=0).*(cdfd-.0001))));
    /* fix up zero's */
            endif;
        endif;
        llnew = llwrk;
        llchng = (llnew-llold)/llold;       /* % change in log-likelihood  */

        if iter gt 1 and _qrsqz eq 0;       /* check if sqeeze is needed  */
            if llchng > -_qrsqtol;
                _qrnsqz0 = _qrnsqz1;
                _qrsqz = 1;
                if llchng gt 0;
                    abold = absqz;
                    continue;
                endif;
            endif;
        endif;
        oldtrap = trapchk(0xffff);
        trap 1;
        abadj = solpd(s,m);         /* modified method of scoring */
        trap oldtrap,0xffff;
        if scalerr(abadj);
            errmsg = "ERROR:  Singular matrix encountered.  Try new start v"\
                "alues.";
            goto errout(error(78));
        endif;

        if _qrnsqz0 > 0;
            gosub squeeze(abold,abadj);     /* adjust step length */
            pop abnew;
            llnew = llwrk;
            llchng = (llnew-llold)/llold;
        else;
            abnew = abold + abadj;
        endif;
        if __output;
            if _qriter == 0;
                if iter == 1;
                    print "ITERATION: ";;
                endif;
                print ftos(iter,"*.*lf ",1,0);;
            else;
                omat = ("LogLik"|(0$+"Alpha_"$+ftocv(seqa(1,1,ncon),1,0))|
                    xlbl);
                if _qriter == 1;
                    output off;
                    call orditer(iter,sqz,(llold|abold),(llnew|abnew),omat);
                    output on;
                elseif _qriter == 2;
                    call orditer(iter,sqz,(llold|abold),(llnew|abnew),omat);
                endif;
            endif;
        endif;
        llold = llnew;
        iter = iter + 1;
        absqz = abold;
        abchng = maxc(abs(abold-abnew));
        abold = abnew;
        ky = upper(chrs(key)) $== "C";
    endo;           /* End of iteration loop */
    abml = abold;

/* Print Results */
    vcml = invpd(m);
    seml = sqrt(diag(vcml));
    tml = abml./seml;
    df = wnused-nparm;
    p_ml = 2*cdftc(abs(tml),df);
    tol = maxc(abs(abchng));
    llfull = 2*llnew;       /* -2 log lik of full model */
    llrest = -2*sumc(nb.*ln(pct));          /* -2 log lik of full model  */
    lrx2 = llrest - llfull;
    df = nvar;
    lrxp = cdfchic(lrx2,df);

    if __output;
        if _qriter == 0;
            print;
            print;
        endif;
        if _qrlogit;
            print ("ESTIMATES FROM ORDINAL LOGIT ANALYSIS OF VARIABLE:  "
                $+ ylbl);
        else;
            print ("ESTIMATES FROM ORDINAL PROBIT ANALYSIS OF VARIABLE:  "
                $+ ylbl);
        endif;
        print;
        if iter > _qrmiter;
            print "   WARNING:  Maximum iterations - " ftos(iter-1,"*.*lf "
                ,1,0) "- exceeded.";
        elseif ky;
            print "   WARNING:  Convergence forced after " ftos(iter-1,"*.*"\
                "lf ",1,0) "iterations.";
        else;
            print "   Convergence after " ftos(iter-1,"*.*lf ",1,0) "iterat"\
                "ions.";
        endif;
        print ("   Tolerance of " $+ ftos(tol,"*.*lf",5,4) $+ " achieved in"\
            "" $+ ftos((hsec-sttime)/6000,"*.*lf",5,2) $+ " minutes.");
        print;
        if _qrlogit;
            print "             Logit        Std.";
        else;
            print "             Probit       Std.";
        endif;
        print "Variable    Estimate     Error      t-value     p>|t|";
        print "-----------------------------------------------------";
        mask = 0~1~1~1~1;
        let fmt[5,3] = "*.*lf" 8 8 "*.*lf " 12 5 "*.*lf " 10 4 "*.*lf " 10
            2 "*.*lf " 10 3;
        omat = xlbl~abml[ncon+1:nparm,1]~seml[ncon+1:nparm,1];
        omat = omat~tml[ncon+1:nparm,1]~p_ml[ncon+1:nparm,1];
        call printfm(omat,mask,fmt);
        print;
        print "Constant    Estimate   Std Error    t-value     p>|t|";
        print "-----------------------------------------------------";
        mask = 0~1~1~1~1;
        omat = 0 $+ "Alpha_"$+ftocv(seqa(1,1,ncon),1,0);
        omat = omat~abml[1:ncon,1]~seml[1:ncon,1];
        omat = omat~tml[1:ncon,1]~p_ml[1:ncon,1];
        call printfm(omat,mask,fmt);
        print;
        print "MEASURES OF FIT:";
        print;
        print ("  Likelihood Ratio Chi-square:            " $+ ftos(lrx2,"*"\
            ".*lf",10,4));
        print ("    with " $+ ftos(df,"*.*lf",1,0) $+ " d.f., prob= " $+
            ftos(lrxp,"*.*lf",3,4));
        print ("  -2 Log Likelihood for full model:       " $+ ftos(llfull,
            "*.*lf",10,4));
        print ("  -2 Log likelihood for restricted model: " $+ ftos(llrest,
            "*.*lf",10,4));
    elseif iter > _qrmiter;