(*:Mathematica:: V2.1 *) (*:Context:: "exlrsys`" *) (*Title:: Coupled lLogistic regression for partially observed categorical data *) (*References: Baker, S.G. 1992. A simple method for computing the observed information matrix when using the EM algorithm with categorical data, Journal of Computational and Graphical Statistics,1, example 4.3 (modified from loglinear to logistic regression*) (*Keywords: nonignorable nonresponse *) (*Requirements: "clm`" *) (*History: 1992 Stuart G. Baker *) Clear[exlrsys,createx1a,createx2a] y=Transpose[{{100,20,30,50,40,60}}]; mod1={"ind", "mcar"}; mod2={"ind", "mar"}; mod3={"sat", "mcar"}; mod4={"sat", "mar"}; mod5={"sat", "mni"}; mod={mod1,mod2,mod3,mod4,mod5}; exlrsys[y_,modelname_]:= Module[{q1,w1,g1,h1,z1,x1, q2,w2,g2,h2,z2,x2, nv,ncore,n,q,w,g,h,z,x, func1,x1a,m1,arg1, func2,x2a,m2,arg2, components,mstepfunc,msteparg, parametername,rationame, c,model,ratio,name}, (*model 1; margin model, i.e., pr(y|x) *) q1=Identity; w1=J[8,1,0]; g1=Id[8]; h1=#2 #1 - Log[1+Exp[#2]]&; z1=Dir[J[4,1],{{1},{0}}]; x1a=createx1a[modelname]; x1=Dir[J[2,1],x1a,J[2,1]]; func1=NLR; m1=Id[4] ~Hcat~ Id[4]; arg1={m1,x1a}; (*model 2: nonresponse model, i.e., pr(r|y) *) q2=Identity; w2=J[8,1,0]; g2=Id[8]; h2=#2 #1 - Log[1+Exp[#2]]&; z2=Dir[{{1},{0}},J[4,1]]; x2a=createx2a[modelname]; x2=J[2,1] ~Dir~ x2a; func2=NLR; m2left=Id[4] ~Dir~ {{1},{0}}; m2right=Id[4] ~Dir~ {{0},{1}}; m2=m2left ~Hcat~ m2right; arg2={m2,x2a}; (*combination*) ncore={{1,1,0,0,1,0},{0,0,1,1,0,1}} . y; nv=Dir[J[2,1],ncore,J[2,1]]; n={nv,ncore}; q={q1,q1}; w={w1,w2}; g={g1,g2}; h={h1,h2}; z={z1,z2}; x={x1,x2}; c=BlockDiag[Id[4],J[1,2],J[1,2]]; components={n,q,w,g,h,z,x}; mstepfunc={func1,func2}; msteparg={arg1,arg2}; model={components,mstepfunc,msteparg}; rmata=J[1,2] ~Dir~ Id[4]; rmatb=J[4,8]; ratio={rmata,rmatb}; ratio="None"; parametername=Automatic; rationame={"p11","p12","p21","p22"}; rationame=Automatic; name={parametername,rationame}; Return[{c,model,ratio,name}]]/;MemberQ[mod,modelname] (*creates design matrix*) createx1a[modelname_]:= Switch[modelname[[1]], "ind", {{1},{1}}, "sat", {{1,0},{1,1}}] createx2a[modelname_]:= Switch[modelname[[2]], "mcar", J[4,1], "mar", J[4,1] ~Hcat~ {{0},{0},{1},{1}}, "mni", J[4,1] ~Hcat~ {{0},{1},{0},{1}}]