fit_ordered_logistic.m

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

function [parms,restparms,fischer,se]=fit_ordered_logistic( link,parms,restparms,design, freq,total,prec)
%fits ordered logistic regression model (cumulative or adjacent categories)


%frequencies: row vector or matrix n by q (q=number of categories)
%total= total counts: column vector of length n
%design is design matrix n*(q-1) by #pars, first all rows for first observation
%(row of freq)
%returns MLE parms (column vector) and optional: fischer info matrix fischer, and standard errors se


%take out parms that are restricted
f=find(restparms~=1);

if ~isempty(f)
    [n,q]=size(freq);

    rel_freq=freq./(total*ones(1,q));
    rel_freq(:,1)=[];  %discard baseline freqs
    rel_freq=rel_freq';
    rel_freq=rel_freq(:);% one vector: first all categories for first observation, etc.
    v=Inf;
    
    
    while v>prec
    %fit logistic model using iteratively reweighted least squares
    %compute working observations and weigths
        
        lin_pred=design*parms;
        lin_pred=reshape(lin_pred,q-1,n);
        if strmatch(link,'cumulative')
        mu=cum_logistic([zeros(n,1) lin_pred']); %add column of zeros 
        
        %mu=multinom_logistic([zeros(n,1) lin_pred']);
        D=deriv_cum_logist(lin_pred');
        %D=deriv_multinom_logist(lin_pred');
        elseif strmatch(link,'adjacent')
                mu=adj_logistic([zeros(n,1) lin_pred']); %add column of zeros 
                D=deriv_adj_logist(lin_pred');
            else error('no valid link function');
            end
        mu(:,1)=[];%remove first column again
        
        m=mu';
        m=m(:);% one vector: first all categories for first observation, etc.
        s=numel(m);
              
        W=[];
        DD=[];
        V=[];

        for i=1:size(D,3)
            varfnctn=(diag(mu(i,:))-mu(i,:)'*mu(i,:))/total(i);
            W=blkdiag(W,squeeze(D(:,:,i))*inv(varfnctn)*squeeze(D(:,:,i))');
            DD=blkdiag(DD,sparse(D(:,:,i)));
            V=blkdiag(V,inv(varfnctn));
        end
        %fischer scoring 
        
        fischer=design'*W*design;
        score=design'*DD*V*(rel_freq-m);
        parms_old=parms;

        f=find(restparms~=1);
        parms(f)=parms_old(f)+fischer(f,f)\score(f);
            
        
        
        %nr_it=nr_it+1;
       % a=find(abs(parms)>50);
        
        %parms(a)=sign(parms(a))*50;
        
        %restparms(a)=1;
        v=max(abs(parms_old-parms));
        %if ~isempty(a) 
        %   v=1;
        %end
    end
end
%parms
if nargout>2
se=sqrt(diag(inv(fischer)));
end