[Fall 2017] Meeting 8 - Example Logistic Regression

Attached are the revised slides for more on Logistic Regression with an example... 
Catie

* all possible pair linearity test using Box-Cox.

* modified Box-Cox equation.
* (y ^ lambda -1)/(lambda * y_geomean ^(lambda-1)) .
* y_geomean * ln (y) .
* draper and smith ( 1998) page 280-282.

preserve.
set printback=off mprint=off.
define !linear (var =  !charend('/') /  width = !default(3) !tokens(1) /groupvar = !tokens(1) /maxloop = !default(100000) !tokens(1)).

set errors = none.
set mxloop = !maxloop.

matrix.

compute ngroup = !length(!groupvar).

do if (ngroup = 0).
+ get data /variables = !var /names = dnames /missing = omit /sysmis = omit.
+ compute totaln = nrow(data).
+ compute p = ncol(data).
+ compute g = make(totaln,1,1).
+ compute gname = 1.
end if.

do if (ngroup > 0).
+ get data /variables = !var !groupvar /names = dnames /missing = omit /sysmis = omit.
+ compute totaln = nrow(data).
+ compute p = ncol(data) - 1.
+ compute gvar = data(:,p+1).
+ compute g = design(gvar). 
end if.

compute totalp = p.
compute p = 2.

compute ng = ncol(g).

loop m = 1 to ng.

compute briefres = make(totalp,totalp,0).

loop pair1 = 1 to totalp.
loop pair2 = 1 to (pair1 - 1). 

compute xvar = data(:,pair1).
compute y = data(:,pair2).

compute x0 = make(totaln,1,1).
compute x = {x0, xvar}.

compute lmin = -1 * !width + 1.
compute lmax = !width + 1.
compute inc = (lmax - lmin)*4 + 1.

compute n = 0.
compute nx = make(totaln,p,0).
compute ny = make(totaln,1,0).

loop i = 1 to totaln.
do if (g(i,m) = 1).
+ compute n = n + 1.
do if (n=1 and ngroup > 0).
+ compute gname = gvar(i,1).
end if.
+ compute nx(n,1:p) = x(i,1:p).
+ compute ny(n,1) = y(i,1).
end if.
end loop.

compute nx = nx(1:n,1:p).
compute ny = ny(1:n,1).

compute j = 0.
compute lambda = make(inc,1,0).
compute rsq = make(inc,1,0).
compute sse = make(inc,1,0).
compute sstotal = make(inc,1,0).
comptue ssreg = make(inc,1,0).
compute logl = make(inc,1,0).
compute ind = make(inc,1,0).

do if (cmin(ny) <= 0).
+ compute ny = ny + (1 - cmin(ny)).
end if.

compute yl = 0.
loop k = 1 to n.
+ compute yl = yl + ln(ny(k)).
end loop.
compute yltemp = yl * (1/n).
compute ylmean = exp(yltemp).

loop i = lmin*4 to lmax*4.
+ compute j = j + 1.
+ compute lambda(j) = i/4.
+ do if (lambda(j) <> 0).
+ compute ylambda = (ny &** lambda(j) - 1) &/ (lambda(j) * ylmean ** (lambda(j)-1)).
+ end if.
+ do if (lambda(j) = 0).
+ compute ylambda = ylmean * ln(ny).
+ end if.
+ compute b = inv(t(nx)*nx)*t(nx)*ylambda.

+ compute e = ylambda - nx*b.
+ compute sse(j) = msum(t(e) * e).
+ compute ylm = msum(ylambda)/n.
+ compute yldif = ylambda - ylm.
+ compute sstotal(j) = t(yldif) * yldif.
+ compute ssreg(j) = sstotal(j) - sse(j).
+ compute rsq(j) = ssreg(j)/sstotal(j).
+ compute logl(j) = -0.5*n*ln(sse(j)/n).

end loop.

compute bestlogl = mmax(logl).
compute llogl = bestlogl - 1.920729.

compute lb = 0.
compute ub = 0.

loop j = 1 to inc.
+ do if (logl(j) >= llogl).
+ compute ind(j) = 2.
+ compute ub = lambda(j).
+ do if (lb = 0).
+ compute lb = lambda(j).
+ end if.
+ end if.
+ do if (logl(j) = bestlogl).
+ compute ind(j) = 1.
+ compute best = lambda(j).
+ end if.
end loop.

do if (lb > 1 or ub < -1).
+ compute briefres(pair2,pair1) = 1.
end if.

end loop. 

end loop. 

do if (ngroup > 0).
+ print gname /title '============================== Group =============================='.
end if.
print briefres /title 'Linearity Matrix - DV = row, IV = column' /rnames=dnames/cnames=dnames /format F8.0.
print /title '0 = linear, 1 = non-linear'.

end loop. 

end matrix.
!enddefine.
restore.