%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The following code performs a bootstrap routine in order to resample gut
%content data, and subsequently estimate beta distribution fits to the
%proportions of each diet item that constitute the diet of a predator. The
%intent of this code is to develop beta distributions (defined by
%alpha and beta parameters) to be used as prior information in the Bayesian
%stable isotope mixing model MixSIR.
% To run the model you must enter in the gut content data, and specify the
% number of mixture data points you are using in the MixSIR program (e.g.
% we used isotope signatures from 8 salmonids in an alaska stream as our
% mix data)
%The fitting function used in this code (dirichlet_fit) is part of the
%Fitfit Matlab toolbox, which also requires the installation of the
%Lightspeed Matlab toolbox. The Dirichlet distribution is a generalized multinomial
%beta distribution and is commonly used to fit proportional data with K
%compartments. I've used it here to fit beta distributions given the cross
%dependencies of the compartments, and subsequently identify the marginals
%of the Dirichlet fit to characterize the individual beta distribs for
%priors.
% Written by Brice X. Semmens (brice.semmens(at)noaa.gov), improvements welcomed! --
% -- This program was last updated 12/03/07 --
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
% Define the number of mixtures used in the MixSIR program
mix = 8;
% Enter the data (e.g. dry-weight of each prey category (sources) as columns,
% and individual fish as rows)
data=[
0 0.021261454 0 0 4.722974128
0.031306164 0.002858017 0 0 0
0 0.003350291 2.470923945 0 0
0.004826265 0.000574232 0.153623732 0 0
0.002066613 0.011165416 0.135291275 0 0
0 0 0 0 10.63018052
0.011336159 0.087127041 0 0 0
0.011185263 0 0 0 0
0.037439171 0.331571928 0 0 0
0.021060915 0.090982954 0 0 0
0.004192756 0 4.451546489 0 0
0.001638277 0.028764563 0.901835813 0 1.80044733
0.040676992 0.137431711 0 0 0
0.00232928 0.001608522 0 0 0
0 0 2.33377449 0 0
0.022326426 0.004433752 0 0 0
0.060150711 0.098914199 0 0 0
0.003542604 0 1.650002248 0 0
0.011273913 0.010123391 0 0 0
0.036903155 0.033277541 0 0 6.033051545
0 0.074786399 0 0 0
0 0.011857462 0 0 0
0 0 0 0 6.804500938
0 0.168523728 0 0 0
0 0.042089001 0 0 0
0 0 0 0 10.81297058
0 0.008905664 0 0 8.508100536
0 0 0 0 7.402102982
0 0.020596039 0 0 0
0.059033578 0.002478583 0 0 0.750489357
0 0.01785594 7.906632943 0 0
0 0 2.33798286 0 0
0.363768347 0.000918603 1.350935166 0 0
0.003456071 0 17.77294249 0 0
0.002104457 0 2.326798266 0 0
0.002009662 0.024726633 0.132863768 0 0
0.005971094 0.000556831 0 0 0
0 0.000961441 1.214004779 0 3.399775542
0 0.046562947 0 0 0
0 0 2.09119578 0 0
0 0.002796883 0 0 11.0416595
0 0 4.186321372 0 0
0.05349736 0 0.737282152 0 0
0 0.062935779 0 0 0
0.020218802 0.028231489 0 0 0
1.40714234 0.004064415 0.684729313 0 0
0.044843571 0.060185078 0 0 4.235892874
0.022307914 0.015662486 0 0 0
0 0.056533075 0 0 0
0 0 0 0 7.158257663
0.015029105 0.013936096 0 0 0.283948236
0 0 4.243093334 0 0
0 0 2.678277499 0 0
0.002384868 0 2.110408821 0 0
0.014237514 0 1.324461356 0 0
0.026748274 0.025335647 0 0 0
0.06154712 0.048488112 4.373640349 0 0
0.077477477 0 2.113513514 0 0
0.005512714 0.010520583 0 0 0
0.004714596 0.046710541 0 0 0
0.093157405 0.28246243 0 0 0
0.014188588 0.01483713 1.902047376 0 0
0.260465116 0 0 0 0
0.185850835 0 0 0 0.540883696
0 0 2.293072454 0 0
0.008318417 0 0 0 0
0 0 3.013168704 0 0
0.074670388 0.024612897 0 0 0
0 0 0 0 11.43177768
0.09059151 0.014422459 0 0 1.596266517
0.029524107 0 6.866286504 0 0
0.028636651 0.002046536 0 0 0
0.008203052 0.016713676 0 0 0
0.009664828 0.002963138 1.468098906 0 0
0.03073075 0.018463393 0 0 0
0.299077046 0.057768337 0 0 0.549974178
0 0.000911918 0 0 0
0.080543788 0 7.24724819 0 0
0 0.009969292 2.553460386 0 0
0 0 2.682773516 0 0
0.050011991 0 1.168029182 0 0
0 0 0.353056563 0 0
0.001040266 0 0 0 11.07699953
0.530769938 0.763338751 0 0 0
0.16969697 0 0 0 0
0.000369975 0 3.701660693 0 0
0 0 0 0 8.130053068
0.031155245 0 0.205449709 0 0
0.009794663 0.000351123 0.09393553 0 0
0 0 0.92244051 0 0
0.000761182 0.032857584 1.724371001 0 0
0.129309755 0 1.50364891 0 0
0.007525988 0.00993422 0 0 0.430571889
0.097731424 0.043752366 0 0 0
0.082807421 0.061973794 0 0 0
0.005828037 0 2.132980293 0 0
0.070442465 0 0 0 0.549479563
0.01705237 0.079935701 0 0 0.212160887
0.002399755 0 4.283140338 0 0
0.018299931 0.099878289 0 0 0
0 0.018004025 1.860181212 0 0
0.003090723 0.153744722 0 0 0
0 0.017801994 2.627581945 0 0
0.079919394 0 1.927487633 0 0
0 0 2.130571257 0 0
0.057928145 0.017189764 0 0 0
0.003982837 0.000674413 0 0 0
0 0 1.589365759 0 0
0.07316844 0.094563755 0 0 0
0.00220855 0.005148077 0 0 0
0 1.794242029 1.899142227 0 0
0.038646956 0.007158339 0.109747691 0 0
0.035172933 0 1.708712804 0 0
0 0 5.159077221 0 0
0.012584794 0.001220246 0 0 0
0.02480315 0.042519685 0 0 0
0.062963196 0.019221005 0 32.4717608 0.249240008
0 0.011718305 0 3.809457652 0
0.008739668 0.001111932 2.134549839 0 0
0.315231788 0.61589404 0 0 0
0 1.327017544 0 0 0
0.032982456 0.027017544 0 0 0
0.016878895 0.000275555 0 0 0
0.038160349 0 0 0 0
0.021016949 0 0.037966102 0 0
0 2.870787246 1.302268955 0 0
0 0 3.663852269 0 0
0.021516139 0 1.343143425 0 0
0 1.999214145 2.787594948 0 0
0.048983957 0 0.767272727 0 0
0.251140915 4.767962107 1.308409039 0 0.160984901
0 0.174234453 1.784621087 0 0
0 0 1.298755705 0 0
0.060581824 0 0.692154562 0 0
0.012926521 0.002971369 0 0 0
0 1.926449336 0.89816576 0 0
0.002084808 2.440169159 3.68976204 0 0
0.001051003 0.738089663 0.446424005 0 0
0.001563572 0.004038319 0 0 0
0.018426591 0.102061178 1.558345023 0 0
0.007372762 0 0.341389198 0 0
0 0.421158903 0 0 0
0 0.088477486 0 0 0
0.000114608 0.000463768 1.501948632 0 0
0.012740644 0.003028726 0.030965195 0 0
0.002253915 0.00617768 0 0 0
0.0225 0 0 0 0
0.002267448 0 0 0 0
0.009028102 0.009081272 1.085723314 0 0
0.024378882 0.198136646 0 0 0
0 0.011016134 0 0 0
0 0.066221461 0 0 0
0 0.001628415 0 0 0
0 0 0.162493969 0 0
0.004548773 0.000442644 0.338932352 0 0
0 1.056766171 1.065285628 0 0
0 0.000985945 1.189149711 0 0
0 0 4.158005665 0 0
0.003663733 0.019989628 0 0 0
0 3.747866195 1.119931098 0 0
0.012470845 0.007909812 0 0 0
0 0.020977608 1.806177223 0 0
0.033228951 0.008617107 0 0 0
0 0.000440739 0 0 0
0 0.000370401 1.334874394 0 0
];
% adding an arbitrarily low number makes it so no value is 0
%(prevents the beta distrib fit from blowing up, and smooths the edges of beta)
data = data+.0001;
% now turn souce dry-weight values into % of gut content for each fish
for i=1:size(data,1)
data(i,:) = data(i,:)./ sum(data(i,:));
end
% now begin the resampling step
total = size(data,1);%number of fish in the data set
stor=zeros(100000,size(data,2));%storage matrix
for samp = 1:100000
draw = ceil(rand(mix,1)*total); %randomly picks 8 rows from data (or however many you specify by changing the mix value_0
stor(samp,:) = sum(data(draw,:));%gets those rows (fish in our case) and sums gut contents across them
stor(samp,:) = stor(samp,:) ./ sum(stor(samp,:));% finds the avg. % of each item in the gut contents of those fish (this sums to one still)
end
%now fit the dirichlet to the avg. % gut content draws made above (this is
%a multinomial beta function-- easily back transformed to beta distribs)
[a] = dirichlet_fit(stor,[1.1 1.5 2.2 1.01 1.01]);
%now separate out each of the beta distribs from the Dirichlet
for i=1:size(data,2)
betaAB(i,:) = [a(i), sum(a)-a(i)]; %stores max like fit alpha and beta parameters for each source (diet item)
end
%now plots the beta distribs to visualize
prop(1) = 0.01;
bdist(1,:) = betapdf(prop(1),betaAB(:,1),betaAB(:,2));
for i =2:99
prop(i) = prop(i-1)+.01;
bdist(i,:) = betapdf(prop(i),betaAB(:,1),betaAB(:,2));
end
plot(prop,bdist');