get_reml_ests.m

function ests = get_reml_ests( data )
%  This code is adapted from the code on the ESA archives from 
%  Staples et al 2004, Ecology, 85(4): 923–929
% 
%  muR, SEmuR, PvR, SvR are trend, standard error of trend,
%  process variance, non-process variance estimates
%  from Staples et al 2004 method
%
%  if either PvR or SvR are equal to 1.0000e-005, these are boundary
%    estimates for the parameters;  see Staples et al 2004
%
%  data is a row vector of population estimates!!
%  data is the raw data not log transformed
%  no missing values allowed

    data = data'; %transform to a column vector
    [years,rdata]=size(data);

    % w is vector of log differeces in population size
    w = log(data(2:years)) - log(data(1:years-1));
    [wr wc]= size(w);
    
    % construct initial estimates of V - cov matrix for w
    %  using method of moments
    sig1hat = var(w);
    covw = cov(w(2:wr),w(1:wr-1));
    sig2hat = covw(1,2);
    sigs = [sig1hat;sig2hat];
    V_hat = diag(sig1hat*ones(wr,1)) + diag(sig2hat*ones(wr-1,1),1)...
        + diag(sig2hat*ones(wr-1,1),-1);
    
    
    %  generate U matrix such that U'*X = 0 
    %  This annhilates fixed effect (mean) for REML estimation
    U = diag(-ones(wr,1)) + diag(ones(wr-1,1),-1);
    U = U(:,1:wr-1);
    
    % create time-series (z) with fixed mean annhilated
    z = U'*w;
    [zr zc]=size(z);
    
    %  conditions for optimization routing -- see MATLAB fmincon function
    %  necessary to restrict the parameter space such that V is invertable
    A = eye(2);
    A(1,1) = -1;
    A(3,:)=[-1 -2];
    A(4,:)=[-1 2];
    b = [0;0;-.00001;-.00001];
    
    %  REML procedure
    options = optimset('LargeScale','off','Display','off',...
        'MaxFunEvals',200,'MaxIter',200);
    x = fmincon('optimize',sigs,A,b,[],[],[0;-2*sig1hat],[2*sig1hat;0],[],options,z,U);
    
    %  model matrix
    X = ones(wr,1);
        
    %  REML estimates for variance paramters
    sig1 = x(1);
    sig2 = x(2);
    
    %  estimate of process variance
    k=[1 2];
    PvR = k*x;
    
    %  estimate of sampling variance
    SvR = -sig2;  
    
    %  calculate Variance for trend estimate
    c2 = sig2*ones(wr-1,1);
    V = diag(sig1*ones(wr,1)) + diag(c2,1) + diag(c2,-1);
    iV = inv(V);
    phi = inv(X'*iV*X);  % var(trend estimate)
   
    %   estimate of trend and SE trend
    muR = phi*X'*iV*w ;  
    SEmuR = sqrt(phi);

	ests.s2p=PvR;
	ests.mu=muR;
	ests.s2s=SvR;
	ests.semuR = SEmuR; %stand error of mu estimate