% Program logregB.m This program uses the function logregA.m
% and the MATLAB function fminsearch to find maximum likelihood
% estimates for a logistic regression model of metapopulation 
% dynamics. 
% As written below, this function is tailored to the particular
% example used in Morris and Doak Chapter 11 (Delphinium
% uliginosum data from the study of Harrison et al. 2000), but
% it contains the code to include effects, such as site sizes,
% that are not used in this example, and also can be modified 
% to use more years of data or more explanatory variables.
clear all;
global FF1 FF2 AA sitenum Dists AAmx NNL randtime ...
 mintime bestparams bestNLL; % define global variables 

%*****************User-defined Parameters*********************
% The program is set up to use two years of data.
% Load the stored information on each site.
% Each row of data should be one site.
% The columns of variables in the data file should be: 
%   X-coordinate of site,
%   Y-coordinate,
%   occupancy state in first year(0=unoccupied, 1=occupied)
%   occupancy state in second year
%   independent variable columns (e.g., size of site)
load delugA.txt;                 
siteinfo = delugA; % assign data to the siteinfo variable       
randtime = 1; % the number of replicate random simulations 
% used to calculate the likelihood of initial occupancy pattern
mintime = 2; % number of initial years to discard in 
% simulations to get rid of transitory dynamics

% There are two arrays of initial parameter values: 
% 1) eeparams stores the values of the extinction parameters
%    ae, betae, se, alpha, and b, used in the expression  
%    Ei=ae + se*Ai + betae*sum{(exp(-alpha*Dij)*pj*Aj^b}
%    where Ai is the area of patch i, Dij is the distance from
%    patch i to patch j, pj indexes occupancy of patch j 
%    (1=occupied, 2=unoccupied), and the sum is taken over 
%    all patches other than patch i.
% 2) ccparams stores the values of the colonization parameters
%    ac, betac, and sc in the expression: 
%	Ci=ac + sc*Ai + betac*sum{(exp(-alpha*Dij)*pj*Aj^b}
% Note that both expressions use the same values for the 
% parameters alpha and b. Also, se, sc, and b (governing site 
% size effects) are not used for Delphinium uliginosum.
eeparams = [0.233   0.000000  -2.862   5.85   0.000000];    
ccparams = [-1.108     0.000000   0.653];  
% IMPORTANT NOTE: starting values can have drastic effects on
%  the final parameter values to which a search routine 
% converges. It is very important to try different starting 
% values in eeparams and ccparams to verify your results. 
%*************************************************************

rand('state',sum(100*clock)); randn('state',sum(100*clock));
% Making variables to use in calculations:
XX = siteinfo(:,1); 	% latitude in km
YY = siteinfo(:,2); 	% longitude in km
FF1 = siteinfo(:,3);	% occupancy state in year 1 (0,1)
FF2 = siteinfo(:,4); 	% occupancy state in year 2 (0,1)
AA = siteinfo(:,5); 	% size of site
sitenum = length(XX);	% number of sites
x1 = XX(:,[ones(1,sitenum)]); 	% making matrices 
x2 = x1';				% of x 
y1 = YY(:,[ones(1,sitenum)]);	% and y values
y2 = y1';
Dists = sqrt((x1-x2).^2 + (y1-y2).^2); % dist’s between sites 
AAmx = AA(:,[ones(1,sitenum)])'; 	% matrix of sizes 
% combine parameters into one array to pass to fminsearch:
params = [eeparams,ccparams]; 
bestparams=params; bestNLL=10000; % these help to track results
% Now, call the MATLAB minimization routine fminsearch to find
% the best values for logregA, a function that calculates the
% negative log likelihood of the data given the parameters. 
fminsearch('logregA',params);  % call fminsearch
disp('Locally optimum parameter values:');
disp('extinction parameters:');
disp([params(1:5)']);
disp('Locally optimum parameter values:');
disp([params(6:8)']);