(*:Mathematica:: V2.1 *) (*:Contex: "exmh`" *) (*Title: square matrix model for ordered marginals and cells 11,12,21,22 with missing data *) (*References: Baker, S.G. (1993) Composite Linear Models for Incomplete Multinomial Data, submitted to Statistics in Medicine *) (*Keywords: marginal homogeniety, ordered categorical data, missing-data mechanism *) (*Requirements: "clm`" *) (*History: 1992 Stuart G. Baker *) Clear[exmh, createx1,createx2, createg1,createm1,createm2, createc]; y=Transpose[{{29,11,8, 5,4,4, 6,9,3, 3,5,0, 2,0,1}}]; mod1={"mh", "ignorable"}; mod2={"mh+shift", "ignorable"}; mod3={"polytomous", "ignorable"}; mod4={"mh", "nonignorable"}; mod5={"mh+shift", "nonignorable"}; mod6={"polytomous", "nonignorable"}; mod={mod1,mod2,mod3,mod4,mod5,mod6}; exmh[y_,modelname_]:= Module[{q1,w1,g1,h1,z1,x1, func1,m1,components1,arg1, q2,w2,g2,h2,z2,x2, func2,m2,components2,arg2, nv,ncore,n,q,w,g,h,z,x, rmata,rmatb, components,mstepfunc,msteparg, parametername,rationame, c,model,ratio,name}, (*sub-model 1 response to A and B*) q1=Log; w1=J[3,1] ~Dir~ (J[8,1,0] ~Vcat~ {{1}}); g1=createg1[modelname]; h1= Exp[#2]/(1+Exp[#2])&; z1=J[8,1]; x1=createx1[modelname]; func1=NGLR; m1=createm1[modelname]; components1={{1},{q1},{w1},{g1},{h1},{z1},{x1}}; arg1={m1,components1}; (*sub-model 2: missing-data mechanism*) q2=Log; w2={{1},{0},{0}} ~Dir~ J[9,1]; g2={{-1,-1},{1,0},{0,1}} ~Dir~ Id[9]; h2= Exp[#2]/(1+Exp[#2])&; z2=J[18,1,0]; x2=createx2[modelname]; func2=NGLR; m2=createm2[]; components2={{1},{q2},{w2},{g2},{h2},{z2},{x2}}; arg2={m2,components2}; (*combination*) ncore=J[1,15].y; nv=ncore ~Dir~ J[27,1]; n={nv,ncore}; q={q1,q2}; w={w1,w2}; g={g1,g2}; h={h1,h2}; z={z1,z2}; x={x1,x2}; c=createc[]; components={n,q,w,g,h,z,x}; mstepfunc={func1,func2}; msteparg={arg1,arg2}; model={components,mstepfunc,msteparg}; (*cell & marginal probabilities*) rmata=Vcat[J[1,3] ~Dir~ Id[9], Dir[J[1,3],Id[3],J[1,3]], J[1,9] ~Dir~ Id[3]]; rmatb=J[15,27]; ratio={rmata,rmatb}; parametername=createparametername1[modelname] ~Join~ createparametername2[modelname]; rationame= {"pi11","pi12","pi13","pi21","pi22","pi23","pi31","pi32","pi33", "pi1+","pi2+","pi3+","pi+1","pi+2","pi+3"}; name={parametername,rationame}; Return[{c,model,ratio,name}]]/;MemberQ[mod,modelname] createc[]:= Module[{cab,ca,cb}, cab=Id[9]; ca=Id[3] ~Dir~ J[1,3]; cb=J[1,3] ~Dir~ Id[3]; c=BlockDiag[cab,ca,cb]; Return[c]]; (*design matrix for sub-model 1*) createx1[modelname_]:= Module[{margins}, margins=Switch[modelname[[1]], "mh", {{1,0},{0,1},{1,0},{0,1}}, "mh+shift", {{1,0,0},{0,1,0},{1,0,1},{0,1,1}}, "polytomous", Id[4]]; x1=BlockDiag[margins,Id[4]]; Return[x1]] (*design matrix for sub-model 2*) createx2[modelname_]:= Module[{des}, des={{1,0},{1,0},{0,1}}; x2=Switch[modelname[[2]], "ignorable", Id[2] ~Dir~ J[9,1], "nonignorable", Vcat[J[3,1] ~Dir~ des, des ~Dir~ J[3,1]]]; Return[x2]] createg1[modelname_]:= Module[{gpi,gim,gt,gmcm,gpcm}, gpi=Id[8] ~Vcat~ J[1,8,-1]; gim={{0,0, 0,0, 1,0,0,0}, {0,0, 0,0, 0,1,0,0}, {1,0, 0,0, -1,-1,0,0}, {0,0, 0,0, 0,0,1,0}, {0,0, 0,0, 0,0,0,1}, {0,1, 0,0, 0,0,-1,-1}, {0,0, 1,0, -1,0,-1,0}, {0,0, 0,1, 0,-1,0,-1}}; gt=Inverse[{{1,0},{1,1}}]; gmcm=Switch[modelname[[1]], "mh", BlockDiag[gt,gt,Id[4]], "mh+shift", BlockDiag[gt,gt,Id[4]], "polytomous", Id[8]]; gpcm= gpi . gim. gmcm; g1=J[3,1] ~Dir~ gpcm; Return[g1]] createm1[modelname_]:= Module[{mto,mtp}, (*creates pairs: {proportion and 1-proportion} for cumulative margins and interior cells*) mto={{1,1,1,0,0,0,0,0,0}, {0,0,0,1,1,1,1,1,1}, {1,1,1,1,1,1,0,0,0}, {0,0,0,0,0,0,1,1,1}, {1,0,0,1,0,0,1,0,0}, {0,1,1,0,1,1,0,1,1}, {1,1,0,1,1,0,1,1,0}, {0,0,1,0,0,1,0,0,1}, {1,0,0,0,0,0,0,0,0}, {0,1,1,1,1,1,1,1,1}, {0,1,0,0,0,0,0,0,0}, {1,0,1,1,1,1,1,1,1}, {0,0,0,1,0,0,0,0,0}, {1,1,1,0,1,1,1,1,1}, {0,0,0,0,1,0,0,0,0}, {1,1,1,1,0,1,1,1,1}}; mtp={{1,1,1,0,0,0,0,0,0}, {0,0,0,1,1,1,1,1,1}, {0,0,0,1,1,1,0,0,0}, {1,1,1,0,0,0,1,1,1}, {1,0,0,1,0,0,1,0,0}, {0,1,1,0,1,1,0,1,1}, {0,1,0,0,1,0,0,1,0}, {1,0,1,1,0,1,1,0,1}, {1,0,0,0,0,0,0,0,0}, {0,1,1,1,1,1,1,1,1}, {0,1,0,0,0,0,0,0,0}, {1,0,1,1,1,1,1,1,1}, {0,0,0,1,0,0,0,0,0}, {1,1,1,0,1,1,1,1,1}, {0,0,0,0,1,0,0,0,0}, {1,1,1,1,0,1,1,1,1}}; m1=J[1,3] ~Dir~ Switch[modelname[[1]], "mh", mto, "mh+shift", mto, "polytomous", mtp]; Return[m1]] createm2[]:= Module[{m2a,m2b,m2ab,ms}, (*creates pairs: {proportion and 1-proportion} for probability of observing only A, only B*) m2a=J[9,1] ~Dir~ {{0,1,0},{1,0,1}}; m2b=J[9,1] ~Dir~ {{0,0,1},{1,1,0}}; m2ab=Vcat[m2a,m2b]; (*position pair for each of 9 responses*) (*old: ms=J[4,1] ~Dir~ Id[9];*) ms=Dir[J[2,1],Id[9],J[2,1]]; m2=Hdir[m2ab,ms]; Return[m2]] (*parameter names for sub-model 1: cell model*) createparametername1[modelname_]:= Module[{margin,interior}, margin=Switch[modelname[[1]], "mh", {"theta1","theta2"}, "mh+shift", {"theta1","theta2","theta3"}, "polytomous", {"theta1","theta2","theta3","theta4"}]; interior={"theta11","theta12","theta21","theta22"}; res=Join[margin,interior]; Return[res]]; (*parameter names for sub-model 2*) createparametername2[modelname_]:= Module[{res}, res=Switch[modelname[[2]], "ignorable", {"thetaA","thetaB"}, "nonignorable", {"theta-s/i", "theta-r"}]; Return[res]]