logit.src
来自「没有说明」· SRC 代码 · 共 1,297 行 · 第 1/3 页
SRC
1,297 行
endif;
call formatcv("*.*lf "~10~8);
call formatnv("*.*lf "~10~4);
if __output;
print "\nDEPENDENT CATEGORIES ARE DESIGNATED AS: \n";
i1 = 1;
i2 = 0;
do until i1 > ncat;
if vtype ==1;
print (" " $+ ftocv(_qrycat[i1],1,0) $+ " - " $+
_qrtmp[i1,1] $+ " ");
else;
print (" " $+ ftocv(i1,1,0) $+ " - " $+ _qrtmp[i1,1] $+ " "\
" ");
endif;
if i2 == 3;
i2 = 0;
print;
else;
i2 = i2 + 1;
endif;
i1 = i1 + 1;
endo;
if (ncat/4) /= trunc(ncat/4);
print;
endif;
endif;
if __output and _qrstat == 1;
print ("\nDISTRIBUTION AMONG OUTCOME CATEGORIES FOR " $+ ylbl);
print;
call printfmt(" "~_qrtmp',zeros(1,rows(_qrtmp)+1));
print "\nPROPORTION ";;
call printfmt(pctall',1);
print;
print "\nDESCRIPTIVE STATISTICS (N="$+ftos(nused,"*.*lf",1,0)$+"):";
/*
if iswt;
print (" (N="$+ftos(wnused,"*.*lf ",1,0) $+
" weighted):");
endif;
*/
print;
let title[1,5] = " " "Mean" "Std Dev" "Minimum" "Maximum";
call printfmt(title,0);
print;
mask = 0~1~1~1~1;
omat = trimr(xlbl~meanx~sd~minx~maxx,1,0);
let fmt[5,3] = "*.*lf " 10 8 "*.*lf " 10 4 "*.*lf " 10 4 "*.*lf "
10 4 "*.*lf " 10 4;
call printfm(omat,mask,fmt);
if rows(omat) eq 1;
print;
endif;
endif;
cv = seqa(2,1,rows(minx)-1);
if maxc(minx[cv] .eq maxx[cv]) == 1;
errmsg = "ERROR: An independent variable has no variation.";
goto errout(error(73));
endif;
/* Compute OLS approximation */
oldtrap = trapchk(0xffff);
trap 1, 1;
bstrt = solpd(xy,xx);
trap oldtrap,0xffff;
if scalerr(bstrt);
errmsg = "ERROR: X'X matrix is singular.";
goto errout(bstrt);
endif;
/* Adjust the OLS start values to approximate logit coefficients */
vc = (eye(ncatm1)-pct).*pct';
dl = pct-vc*(ln(pct/pctlast));
vv = zeros(nivar,ncatm1);
vv[1,.] = ones(1,ncatm1).*(dl');
/* bmat = (bstrt-vv)*invpd(vc);*/
bmat = solpd((bstrt-vv)',vc)';
clear vv,dl;
/* ML iterations using the method of scoring */
if _qriter == 0 and __output;
print;
print "ITERATION: ";;
endif;
bchng = 1000;
iter = 0;
bvec = vec(bmat);
clear ypred;
ky = upper(chrs(key)) $== "C";
do until abs(bchng) < __tol or iter >= _qrmiter or ky;
m = zeros(nparms,nparms);
tm = m;
tzu = zeros(nparms,1);
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;
if iswt;
wt = dta[.,cols(dta)];
dta = trimr(dta',0,1)';
else;
wt = ones(rows(dta),1);
endif;
/* prepare data for logit analysis */
x = ones(rows(dta),1)~trimr(dta',1,0)';
if vtype;
ydum = trunc(dta[.,1]) .== ycat1;
else;
ydum = dta[.,1] .$== ycat1;
endif;
gosub compute;
endo;
clear dta;
else;
gosub compute;
endif;
/* Update estimates */
oldtrap = trapchk(1);
trap 1,1;
/* vc = invpd(tm);*/
bchng = solpd(tzu,tm);
trap oldtrap,1;
/* if scalerr(vc);*/
if scalerr(tm);
errmsg = "\nERROR: Iterations failed.";
goto errout(error(35));
endif;
/* bnew = bvec+vc*tzu;*/
/* bchng = bnew-bvec; */
bnew = bvec + bchng;
iter = iter + 1;
if __output;
if _qriter .eq 0;
print " " ftos(iter,"*.*lf ",1,0);;
else;
if ncat==2;
omat = reshape(xlbl,nparms,1);
else;
omat = reshape(xlbl~reshape(""|" ",nivar,ncat-2),
nparms,1);
endif;
if _qriter == 1;
output off;
call logiter(iter,vec(reshape(bvec,ncatm1,nivar)),
vec(reshape(bnew,ncatm1,nivar)),omat);
output on;
elseif _qriter == 2;
call logiter(iter,vec(reshape(bvec,ncatm1,nivar)),
vec(reshape(bnew,ncatm1,nivar)),omat);
endif;
endif;
endif;
bvec = bnew;
bmat = reshape(bvec,ncatm1,nivar)';
ky = upper(chrs(key)) $== "C";
endo; /* end of ML iterations loop */
goto convrged;
compute:
xb = x*bmat;
if xb < 500 and xb > -500;
fl = exp(xb);
ypred = fl./(1+sumc(fl')); /* sumc(fl') : add up all
:: categories
*/
resid = ydum-ypred;
ii = 1;
do until ii > ncatm1;
ij = ii;
do until ij > ncatm1;
if ii == ij;
s = ypred[.,ii].*(1-ypred[.,ii]);
cc = 1;
else;
s = ypred[.,ii].*ypred[.,ij];
cc = -1;
endif;
i2 = nivar*ii;
i1 = i2 + 1 - nivar;
j2 = nivar*ij;
j1 = j2 + 1 - nivar;
sx = sqrt(wt.*s).*x;
m[i1:i2,j1:j2] = cc*moment(sx,0);
ij = ij + 1;
endo;
ii = ii + 1;
endo;
tm = tm + m;
tzu = vec((wt.*x)'resid) + tzu;
return; /* end of compute subroutine */
else;
errmsg = "\nERROR: Iterations failed due to bad intermediate estim"\
"ates of coefficients.";
goto errout(error(36));
endif;
convrged:
tol = maxc(bchng); /* tolerance at convergence */
bnew = bmat;
bprt = vec(reshape(bvec,ncatm1,nivar)); /* b in order to print */
if _qrev;
bprt = -bprt;
endif;
vc = solpd(eye(rows(tm)),tm);
local vci;
vci = vec(reshape(seqa(1,1,rows(bprt)),ncatm1,nivar));
vcml = vc[vci,vci];
sdml = sqrt(diag(vcml));
tml = bprt./sdml;
pml = 2*cdfnc(abs(tml));
if __output;
print;
print;
print ("ESTIMATES FROM LOGIT ANALYSIS OF VARIABLE: " $+ ylbl);
print;
if iter >= _qrmiter;
print " WARNING: Maximum iterations - " ftos(iter,"*.*lf ",
1,0) "- exceeded.";
elseif ky;
print " WARNING: Convergence forced after " ftos(iter,"*.*lf"\
" ",1,0) "iterations.";
else;
print " Convergence after " ftos(iter,"*.*lf ",1,0) "iteratio"\
"ns.";
endif;
print (" Tolerance of " $+ ftos(tol,"*.*lf ",5,4) $+ "achieved af"\
"ter" $+ ftos((hsec-sttime)/6000,"*.*lf ",5,2) $+ "minutes.");
elseif iter >= _qrmiter;
errorlog "WARNING: Maximum iterations exceeded";
elseif ky;
errorlog "WARNING: convergence forced";
endif;
/* Compute LRX2 for each variable: Aldrich & Nelson: 59-60. */
omat2 = 0~0~0;
iv = 1;
do until iv == nivar+1;
ii = seqa(iv,nivar,ncatm1);
covar = submat(vc,ii,ii);
bj = submat(bvec,ii,1);
lrx2 = bj'*invpd(covar)*bj;
lrxp = cdfchic(lrx2,ncatm1);
omat2 = omat2|(lrx2~ncatm1~lrxp);
iv = iv+1;
endo;
if __output;
print "\n Logit Std "\
" 2-tailed Exp\n Variable Comparison Estimate "\
"Error t-value Prob Estimate\n -----------------"\
"------------------------------------------------------------"
;
mask = 0~0~1~1~1~1~1;
let fmt[7,3] = "*.*lf " 10 8 "*.*s" 8 8 "*.*lf " 14 5 "*.*lf " 10
4 "*.*lf " 10 2 "*.*lf " 8 3 "*.*lf " 12 4;
if maxc(bprt) > 10;
fmt[7,1] = "*.*le";
endif;
if ncat==2;
omat = reshape(xlbl,nparms,1);
else;
omat = reshape(xlbl~reshape(""|" ",nivar,ncat-2),nparms,1);
endif;
if vtype;
if _qrev;
omat = omat~reshape(ftocv(_qrycat[rows(_qrycat)],1,0) $+ "/"\
"" $+ ftocv(_qrycat[1:ncatm1],1,0),nparms,1);
else;
omat = omat~reshape(ftocv(_qrycat[1:ncatm1],1,0) $+ "/" $+
ftocv(_qrycat[rows(_qrycat)],1,0),nparms,1);
endif;
else;
if _qrev;
omat = omat~reshape(_qrycat[rows(_qrycat)] $+ "/" $+
_qrycat[1:ncatm1],nparms,1);
else;
omat = omat~reshape(_qrycat[1:ncatm1] $+ "/" $+
_qrycat[rows(_qrycat)],nparms,1);
endif;
endif;
omat = omat~bprt~sdml~tml~pml;
call printfm(omat~exp(omat[.,3]),mask,fmt);
endif;
/* Compute likelihood ratio statistics */
ll = 0;
if _qrpred;
/* Open file for predicted probabilities */
_qrpredn = "" + _qrpredn;
if inm or vtype == 0;
prednm = (0 $+ "P_" $+ _qrtmp)|ylbl;
else;
prednm = (0 $+ "P_" $+ ftocv(_qrycat,1,0))|ylbl;
endif;
create fpred = ^_qrpredn with ^prednm,ncat+1,8;
if fpred == -1;
errmsg = "ERROR: Can't open file for residuals: " $+ _qrpredn
$+ ".";
goto errout(error(74));
endif;
endif;
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;
if iswt;
wt = dta[.,cols(dta)];
dta = trimr(dta',0,1)';
else;
wt = ones(rows(dta),1);
endif;
/* prepare data for logit analysis */
ynum = dta[.,1];
if vtype;
ynum = trunc(ynum);
ydum = ynum .== ycat1;
else;
ydum = ynum .$== ycat1;
endif;
x = ones(rows(dta),1)~trimr(dta',1,0)';
gosub pred;
if _qrpred; /* write observed outcome and probabilities */
if writer(fpred,g~ynum) /= rows(ynum);
errmsg = "ERROR: Disk full, predicted file not complet"\
"ed.";
goto errout(error(75));
endif;
endif;
endo;
clear dta;
else;
gosub pred;
if _qrpred; /* write observed outcome and probabilities */
if writer(fpred,g~ynum) /= rows(ynum);
errmsg = "ERROR: Disk full, predicted file not written.";
goto errout(error(75));
endif;
endif;
endif;
fin = close(fin);
/* Compute goodness of fit measures */
nb = wnused*pctall;
/* log-likelihood for all parameters */
llfull = ll;
/* log-likelihood with only intercepts */
llrest = sumc(nb.*ln(pctall));
/* likelihood ratio test statistics - Maddalla, p.40 */
lrx2 = -2*(llrest - llfull);
df = (nivar-1)*(ncatm1);
lrxp = cdfchic(lrx2,df);
/* McFadden's pseudo-R2 - Madalla, p.40 */
mcfr2 = 1-(llfull/llrest);
/* Madalla's R2 - Madalla, p. 39 */
madr2 = 1-exp(-lrx2/wnused);
/* Cragg & Uhler's pseudo-R2 - Madalla, p.40 */
ii = 2/wnused;
cur2 = (exp(ii*llfull)-exp(ii*llrest))/(1 - exp(ii*llrest));
obspred = reshape(obspred',ncat,ncat)./wnused;
/* percent correct predictions - Maddall, p.70 */
success = 100*sumc(diag(obspred));
obspred = wnused*(obspred~sumc(obspred'));
if __output;
omat2[1,1] = lrx2;
omat2[1,2] = df;
omat2[1,3] = lrxp;
print "\nMEASURES OF FIT:\n";
if ncat == 2;
print (" Likelihood Ratio Chi-square: " $+
ftos(lrx2,"*.*lf ",10,4));
print (" with " $+ ftos(df,"*.*lf ",1,0) $+ " d.f., prob="
$+ ftos(lrxp,"*.*lf ",4,3));
else;
print " Test LRX2 df Prob";
print " --------------------------------------";
mask = 0~1~1~1;
omat = ("Overall"|xlbl)~omat2;
let fmt[4,3] = "*.*lf " 10 8 "*.*lf " 12 4 "*.*lf " 6 0 "*.*lf "\
"" 9 3;
call printfm(omat,mask,fmt);
print;
endif;
print (" -2 Log Likelihood for full model: " $+
ftos(-2*llfull,"*.*lf ",10,4));
print (" -2 Log likelihood for restricted model: " $+
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?