fit_model_mixed_lltm.m

The BNL toolbox is a set of Matlab functions for defining and estimating the parameters of a Bayesi

%main file for estimation prob graphical models
%mixed lltm (rasch) 
%see the examples in the example map for other specific models
%or define your own model
%to be done for each model
    %1.load data 
    %2.define bayes net (model structure)
    %3.define equivalence sets of nodes, and link function and design
        %matrix for each equivalent set
    %4. estimation, automatic after specification of 1-3
    %5. postestimation: compute standard errors, save results 

%%%%%
%1.load data
%%%%%
    %obs_variables: the observed variables that are nodes in the
    %graphical model
        %obs_variables is structure :
            %obs_variables.names: cell array of names (does not have to be
        %defined)
            %obs_variables.datamatrix: 'wide' format datamatrix, there is a
                %separate row for each independent case (e.g. two repeated
                %measurements defines two separate variables)
                %missings are coded by the value -1
        
        %matlab import wizard can be used to import data from e.g. excel
        obs_variables.datamatrix=load('spo1_it_b.txt');
    
    %covariates: the covariates, variables that vary over cases
        %but that are themselves not included as separate nodes in the network
        %covariates.names: cell array of names
        %covariates.datamatrix: again 'wide' format.
         %no case covariates
        

        
%%%%%       
%2.define bayes net (model structure)
%%%%%

    %specify bayes net using one of the construct_bnet functions or use your own

    %specifications of options for the construct_bnet function that is used. 
    S=3;% number of latent states 
    S_item=2;%number of response categories for items
    nrquadpoints=20;%number of locations for gaussian quadrature
    [N,I]=size(obs_variables.datamatrix); %number of cases
    %bnet construction
    [bnet,evidence_nodes,partial_evidence_nodes,terminal_merged_nodes,hid_nodes,gausskwadnodes,...
    names,onames,order,inv_o]=...
    construct_bnet_mixlltm(S,nrquadpoints,I,S_item,obs_variables.datamatrix);
    %convert bnet object from BNT (Murphy) into what we need 
    [parents,child,node_sizes,postorder,preorder,cliquetable,septable,clq_ass_to_node,pot_to_CPT]= franks_from_BNT(bnet,N);
    

%%%%
%3.define equivalence sets of nodes, link function,  design
    %matrix and restrictions on parameters for each equivalent (parameter) set
%%%%
    %equivalent nodes are nodes with the same design matrix, the same link and governed by the same set of parameters 
    %param_equivalent nodes are nodes with a differnetn design matrix, but
    %the same link and governed by the same set of parameters
    equiv_class=zeros(nr_nodes,1);
    equiv_class(strmatch('class', onames))=1;
    equiv_class(strmatch('theta', onames))=2;
    equiv_class(terminal_merged_nodes.nodenrs)=3;
     
    param_equiv_class=equiv_class;

    nr_equiv=max(equiv_class);
    for i=1:nr_equiv
        link{i}='multinomial';
    end
    
    %design matrix for all equiv classes
        %first, define pred_mat which is a cell array. each cell contains
        %the design matrix of the
        %variable(s) in a node without modelling the parents. e.g. 
        %for a four-category variable, a 3*3 design matrix
        %case covariates are not yet included here 
        %default is an identity matrix
        pred_mat=construct_predmat(equiv_class,evidence_nodes,partial_evidence_nodes,terminal_merged_nodes,hid_nodes,gausskwadnodes);
        %node number for obs variables is 3, we have design matrix for them
        des=load('desLLTMspo1.txt');
        pred_mat{3}=des(:,1:end-1);%rank is only six, skip last column of des (ones)
        %second, define the structure of the linear predictor for each equiv class (no restrictions, only main effects for parents etc)    

        lin_pred_struct=cell(max(equiv_class),1);
        %first two are theta and class
        lin_pred_struct(1)={[0 ]};
        lin_pred_struct(2)={[0 ]};
        %check order of theta and class
        if gausskwadnodes.nodenrs<hid_nodes.nodenrs
            a={[0 2] ;1 };
        else
            a={[0 1];2  };
        end
        lin_pred_struct(3)={a};

   
        %third, construct design matrix for each node based on its pred_mat and  linear
        %predictor

        design=construct_design_mats(parents,node_sizes,equiv_class,lin_pred_struct,pred_mat,gausskwadnodes);
        
        
        %fourth, inclusion of case covariates, makes design three dimensional, last
        %dimension represents cases
        %no covariates 
        
            
        %fifth, further changes to design matrices can be done manually or
        %with code provided by user
            
    % define restrictions on parameters (gauss kwad nodes are
    %restricted automatically in gen_rand_start function below)
    restparms=cell(max(param_equiv_class),1);
       
    %default =0, no restrictions (1 =restricted at fixed value)
    %restrictions for gauss quad nodes are defined in gen_random_start
    for i=1:max(param_equiv_class);
        node_nr=find(param_equiv_class==i,1);
        restparms{i}=zeros(size(design{equiv_class(node_nr)},2),1);
    end
   
%%%%
%4. estimation, automatic after specification of 1-3
%%%%
    %estimation with EM
    %can be altered: number of runs, starting values, convergence criteria

    %number of runs
    nr=1; 
    for run=1:nr
        loglikelihood=-1e10;
        v=Inf;
        it=1;
        %starting values   (here: random start)
        [parms,restparms]=gen_random_start(design,parents,equiv_class,param_equiv_class,gausskwadnodes,link);
        loglikelihood=-1/eps;
        logl_diff=1;
        logl_avg=1;
        %EM loop
        while (v>.00000001 && logl_diff/logl_avg>.00000001 )
            %1 iteration
            parms_old=parms;
            loglikelihood_old=loglikelihood;
            [parms, restparms, loglikelihood ]=EM_iteration(link,parms,restparms,design,parents,node_sizes...
            ,equiv_class,param_equiv_class,evidence_nodes,partial_evidence_nodes,terminal_merged_nodes,hid_nodes,...
            gausskwadnodes,preorder, postorder ,cliquetable, septable,pot_to_CPT,N);    
            %check convergence
            v=0;
            for i=1:max(param_equiv_class)
                vvec=abs(parms_old{i}-parms{i});
                v=[v max(vvec)];
                v=max(v);
            end
            v
            it=it+1;
            logl_diff=loglikelihood-loglikelihood_old;
            if logl_diff<-.0001
                disp('decrease in loglikelihood. press any key to continue');
                pause;
            end
            logl_avg=(abs(loglikelihood)+abs(loglikelihood_old))/2
        end
%%%%
%5. postestimation: compute standard errors, save results ?
%%%%

        %compute final CPTs?
    	lin_pred=construct_lin_pred(parms,design,parents,node_sizes,equiv_class,param_equiv_class,terminal_merged_nodes,N);
        equiv_class_CPTs=construct_equiv_class_CPT(link,lin_pred,N);
    
        
        %compute standard errors?
        info=num_infomatrix_anal_score(link,parms,restparms,design,parents,node_sizes,equiv_class,...
        postorder,septable, cliquetable,preorder,param_equiv_class,clq_ass_to_node,evidence_nodes,...
        partial_evidence_nodes,terminal_merged_nodes,hid_nodes,gausskwadnodes,N,.000001);
        serr=sqrt(diag(inv(info)));
    end