(*:Mathematica:: V2.2.1 *) (*:CLM Version:: 1.1 *) (*:Context:: "clmdesig`" *) (*:Title:: Design and transformation matrices for CLM *) (*:Summary:: This package contains the functions for creating special design matrices (X) and transformation matrices (G) for Composite Linear Models. CategoricalDesignMatrix: creates design matrix for loglinear model HierarchicalInterations: creates a list of hierarchical interactions FailureMatrix: creates failure transformation matrix CensoringMatrix creates censoring transformation matrix *) (*:Keywords:: design matrix, loglinear model, IPF *) (*:Input For:: clmem`, clmnipf` and some example functions*) (*:Requirements:: matrixlx: matrix functions *) (*:Author:: 1992 Stuart G. Baker *) BeginPackage["clm`clmdesig`", "clm`matrixlx`"] CategoricalDesignMatrix::usage= "CategoricalDesignMatrix[dim,configset,options] creates a design matrix for the indicated loglinear model. For example, CategoricalDesignMatrix[{2,3,2},{{1,2},{2,3}}] creates a design matrix for three variables 1 (2 levels), 2 (3 levels) and, 3 (2 levels), where the highest order interactions are between variables 1 and 2 and between variables 2 and 3." (*options*) Parametrization::usage= "Parametrization is an option of CategoricalDesignMatrix. Parametrization-> \"ZeroOne\" specifies a design matrix of 0's and 1's. Parametrization-> \"ZeroSum\" specifies a design matrix of -1's and 1's which constrain sums of parameters of each order to equal 0, as in ANOVA." Baseline::usage= "Baseline is an option of CategoricalDesignMatrix, which specifies which level of each variable is a baseline. For example, Baseline-> {2,1,2} specifies baselines for level 2 of variable 1, level 1 of variable 2, and level 2 of variable 3. The defaut is Baseline-> \"First\",which specifies that the first level of each variable is the baseline. The baseline under the \"ZeroOne\" parametrization is a row of zeros; the baseline under the \"ZeroSum\" parametrization is a row of -1's." ModelType::usage= "ModelType is an option of CategoricalDesignMatrix. If ModelType-> \"Hierarchical\" higher order interactions imply lower-order interactions; otherwise, ModelType->\"NonHierarchical\"." VariablesVary::usage= "VariablesVary is an option of CategoricalDesignMatrix. VariablesVary-> \"SlowtoFast\" specifies that variable 1 varies slowest then 2, etc. This is the default option which should be used with NIPF or NIPFCount. VariablesVary->\"FasttoSlow\" specifies that that varible 1 varies fastest,then 2, etc." Options[CategoricalDesignMatrix] = {Baseline->"First", Parametrization->"ZeroOne", ModelType->"Hierarchical", VariablesVary->"SlowtoFast"} clmdesig::baddim= "The dimensions should not contain a 1." clmdesig::badconf= "The configuration is not consistent with the the dimensions." clmdesig::badpara= "The Parametrization option is not \"ZeroOne\" or \"ZeroSum\"." clmdesig::badvv= "The VariablesVary option is not \"FasttoSlow\" or \"SlowtoFast\"." clmdesig::badtype="The ModelType option is not \"Hierarchical\" or \"NonHierarchical\". " clmdesig::badbase= "The baseline is not consistent with the dimensions." (*---------------HierarchicalInteractions-------------------------*) HierarchicalInteractions::usage= "HierarchicalInteractions[config] creates a list of hierarchical interactions. For example, Hierarchical Interactions[{{1,2},{2,3}] returns {{1},{2},{3},{1,2},{2,3}}." (*--------------------FailureMatrix---------------------------*) FailureMatrix::usage= "FailureMatrix[number of times, list of failure times, list of censoring times] returns matrix for failures. See exhaz`." (*--------------------CensoringMatrix-------------------------*) CensoringMatrix::usage= "CensoringMatrix[number of times, list of failure times, list of censoring times] retruns matrix for censoring. See exhaz`" Clear[CategoricalDesignMatrix,HierarchicalInteractions, FailureMatrix,CensoringMatrix]; Begin["Private`"] (*------------------CategoricalDesignMatrix--------------------*) CategoricalDesignMatrix[dim:{___Integer?Positive},configset_:Automatic, options___Rule]:= Module[{baseline,param,modeltype,configset1,baseline1,matrix}, {baseline,param,modeltype,vv} = {Baseline,Parametrization,ModelType,VariablesVary}/. {options}/.Options[CategoricalDesignMatrix]; If[configset===Automatic, configset1={Range[Length[dim]]}, configset1=configset]; If[baseline==="First", baseline1=Flatten[J[1,Length[dim]]], baseline1=baseline]; matrix=pCategoricalDesignMatrix[dim,configset1,modeltype, param,baseline1,vv]; Return[matrix]]/; Which[MemberQ[dim,1], Message[clmdesig::baddim];False, True, True] pCategoricalDesignMatrix[dim_,configset_,modeltype_,para_,base_,vv_] := Module[{mainset,configset1,configset2, x,x1,x2,i,j}, (* create sets of main effect columns *) mainset=Table[pMainEffectColumn[dim,i,para,base,vv], {i,1,Length[dim]}]; (* write out configset *) If[configset==="Saturated", configset1={Range[Length[dim]]}, configset1=configset]; If[modeltype==="Hierarchical", configset2=HierarchicalInteractions[configset1], configset2=configset1]; (* create design matrix from configset and sets of main effects*) x=Table[pInteractionColumn[mainset,configset2[[j]]], {j,1,Length[configset2]}]; x1=Join[{J[Times @@ dim, 1]}, x]; x2=Hcat[Sequence @@ x1]; Return[x2]]/; Which[!(pDimConfigsetQ[dim,configset]), Message[clmdesig::badconf];False, !(pParametrizationQ[para]), Message[clmdesig::badpara];False, !(pVariablesVaryQ[vv]), Message[clmdesig::badvv];False, !(pModelTypeQ[modeltype]), Message[clmdesig::badtype];False, !(pDimBaselineQ[dim,base]), Message[clmdesig::badbase];False, True, True]; (*test if configuration and dimensions are consistent*) pDimConfigsetQ[dim_,configset_]:= And @@ (MemberQ[Range[Length[dim]],#]& /@ Flatten[configset]) (*test if valid parametrization*) pParametrizationQ[para_]:= (para === "ZeroOne") || (para === "ZeroSum") (*test if valid modeltype*) pModelTypeQ[modeltype_]:= (modeltype === "Hierarchical") || (modeltype === "NonHierarchical") (*test if valid variables vary*) pVariablesVaryQ[vv_]:= (vv === "FasttoSlow") || (vv === "SlowtoFast") pDimBaselineQ[dim_,baseline_]:= Module[{i}, And @@ Table[MemberQ[Range[dim[[i]]],baseline[[i]]], {i,1,Length[dim]}] ] (* create core unit of parametrization from number of levels pMainEffectCore[number of levels,parametrization,baseline], for example, pMainEffectCore[3,"ZeroOne",1] returns {{0,0},{1,0},{0,1}} and pMainEffectCore[3,"ZeroSum",1] returns {{-1,-1},{1,0},{0,1}}.*) pMainEffectCore[s_,para_,base_] := MoveRow[Join[ J[1,s-1,0], Id[s-1]],1,base]/; (para==="ZeroOne") pMainEffectCore[s_,para_,base_] := MoveRow[Join[ J[1,s-1,-1], Id[s-1]],1,base]/; (para==="ZeroSum") (*main effect based on dimension of variables; pMainEffectColumn[dimensions,variable,parametrization,baseline] for example, pMainEffectColumn[{2,3},1,"ZeroOne",{1,1}] returns {{0,0},{1,0},{0,1},{0,0},{1,0},{0,1}}.*) pMainEffectColumn[dim_,s_,para_,base_,vv_] := Module[{a,b,c,d,dimx,n,i,j}, dimx=Join[{1},dim,{1}]; n=Length[dimx]; a=Product[dimx[[i]],{i,s+2,n}]; b=dimx[[s+1]]; c=base[[s]]; d=Product[dimx[[j]],{j,1,s}]; If[vv=="FasttoSlow", y=Dir[ J[a,1], pMainEffectCore[b,para,c], J[d,1] ], y=Dir[ J[d,1], pMainEffectCore[b,para,c], J[a,1] ]]; Return[y]]; (* horizontal direct product of main effects indicated by interactions *) pInteractionColumn[mainset_,inter_]:= Hdir[Sequence @@ mainset[[inter]]] (*---Hierarchical Interactions-------------------------------*) HierarchicalInteractions[configset:{__List}] := Module[{n,all,i}, n=Length[configset]; all=Table[pSubsets[configset[[i]]],{i,1,n}]; interactions=Union[Flatten[all,1]]; Return[interactions]] (*create list subsets excluding null set*) (*from Stan Wagon,Mathematica Journal Spring 1991, page 46 *) pSubsets[x_List]:= Drop[Sort[Flatten /@ Distribute[{{},{#}}& /@ x,List]],1] (*--------------FailureMatrix---------------*) FailureMatrix[n_,v1_List,v2_List]:= Join[pFailureMatrixF[n,v1],pFailureMatrixC[n,v2]] (*matrix for failure hazards for failure at times v*) pFailureMatrixF[n_,v_]:= Module[{a,b,i,j,tab}, a=Length[v]; b=2 n; tab=Table[If[(EvenQ[j] && j < 2 v[[i]]) || j==2 v[[i]] - 1, 1, 0], {i,1,a},{j,1,b}]; Return[tab]] /; n >=Max[v] (*matrix for failure hazards for censoring at times v*) pFailureMatrixC[n_,v_]:= Module[{a,b,i,j,tab}, a=Length[v]; b=2 n; tab=Table[If[EvenQ[j] && (j <= 2 v[[i]]), 1, 0], {i,1,a},{j,1,b}]; Return[tab]] /; n >=Max[v] (*----------CensoringMatrix----------------------*) CensoringMatrix[n_,v1_List,v2_List]:= Join[pCensoringMatrixF[n,v1],pCensoringMatrixC[n,v2]] (*matrix for censoring hazards for failures at times v*) pCensoringMatrixF[n_,v_]:= Module[{a,b,i,j,tab}, a=Length[v]; b=2 n; tab=Table[If[EvenQ[j] && (j < 2 v[[i]]), 1, 0], {i,1,a},{j,1,b}]; Return[tab]] /; (n+1) >=Max[v] (*matrix for censoring hazards for censorings at times v*) pCensoringMatrixC[n_,v_]:= Module[{a,b,i,j,tab}, a=Length[v]; b=2 n; tab=Table[If[(EvenQ[j] && j < 2 v[[i]]) || j==2 v[[i]] - 1, 1, 0], {i,1,a},{j,1,b}]; Return[tab]] /; (n+1) >=Max[v] End[] (*Protect[CategoricalDesignMatrix,HierarchicalInteractions, FailureMatrix,CensoringMatrix];*) EndPackage[]