load WGNScans.txt /ascii;
d=WGNScans;

% data columns
%1) Bag
%2) Julian Day
%3) Sampling round
%4) Uninfected juveniles
%5) Uninfected adults
%6) Metsch-infected juveniles (there were a few)
%7) Metsch-infected adults
%8) Juvenile Density (per L)
%9) Adult Density (per L)
%10) Total Density (per L)
%11) Treatment (0, 0.25, 0.5, 1, or 2)
%12) Secchi depth

bag=d(:,1);
day=d(:,2);
juv=d(:,8);
adult=d(:,9);
pinfectjuv=d(:,6)./(d(:,6)+d(:,4));
pinfectadult=d(:,7)./(d(:,7)+d(:,5));

%replace NaN with 0
pinfectadult(find(isnan(pinfectadult)))=0;

% log transform after replacing 0 with 1/2 lowest observed value
juvmin=min(juv(find(juv>0)));
juv(find(juv==0))=juvmin/2;
juv=log(juv);

adultmin=min(adult(find(adult>0)));
adult(find(adult==0))=adultmin/2;
adult=log(adult);

density=d(:,11)+d(:,12);
densitymin=min(density(find(adult>0)));
density(find(density==0))=densitymin/2;
density=log(density);

numinfected=d(:,6)+d(:,7);
denom=d(:,4)+d(:,5)+d(:,6)+d(:,7);
pinfect=numinfected./denom;
pinfect(find(isnan(pinfect)))=0;

treat=d(:,14);
Secchi=d(:,15);

pinfect=sqrt(pinfect);
pinfect=asin(pinfect);

pinfect=pinfect/cov(pinfectadult)^.5;
density=density/cov(density)^.5;
Secchi=Secchi/cov(Secchi)^.5;

%%%%%%%%%%%%% plot data %%%%%%%%%%%%%%%%%%%%%%%
if 1,
	pinfecteds=[];
	densities=[];
	Secchis=[];
	for i=1:10
		pinfecteds=[pinfecteds pinfect(find(bag==i))];
		densities=[densities density(find(bag==i))];
		Secchis=[Secchis Secchi(find(bag==i))];
	end
	figure(1)
	subplot(3,1,1)
	plot(pinfecteds);

	subplot(3,1,2)
	plot(densities);

	subplot(3,1,3)
	plot(Secchis);
end
	
%%%%%%%%%%%%%%%%% construct data sets %%%%%%%%%%%%%%%%%%%%%%

% set up data for CLS
% This w matrix has the bag, then all of the data from t=1:T, then all of
% the data from t=0:T-1; that is, it has the bag, then Y, then X (using
% notation of Ives et al. 2003)
w=[bag density Secchi pinfect treat];
w=[w(2:end,:) w(1:end-1,:)];

%Removing points where have mismatch between bags:
w=w(find(w(:,1)==w(:,6)),:);

bagX=w(:,1);
Y=w(:,2:4);
X=w(:,7:10);

[n p]=size(Y)

%%%%%%%%%%%%%%%%%% CLS estimates %%%%%%%%%%%%%%%%%
% designmatrix indicates which comparisons are of interest.
% For example, the first row says we're interested in the effects of juvenile density
% and infection prevalence at t on juvenile density at t+1, but not the
% effects of adult density at t on juvenile density at t+1.
designmatrix=[1 0 1 1;
			  1 1 0 1;
              0 0 1 1];
% This B matrix has four rows; the first row is for the intercept; the
% other rows are for the slope parameters:
B=zeros(5,3);	
E=[];
for i=1:3 
    % This sets up just the Y vector of interest (for example, juvenile
    % density):
	y=Y(:,i);
    % This makes a matrix with the first column of ones, and then
    % subsequent columns with the vectors X that are of interest.
    % Continuing the example where y=juvenile density, x will be made up of
    % a vector of ones, the juvenile density vector and the proportion
    % adults infected vector, where, again, the x vectors are from t=0:T-1
    % and the y vectors are from t=1:T
	x=[ones(length(Y),1) X(:,find(designmatrix(i,:)==1))];

    % Slope parameter:
	b=(x'*x)\(x'*y);
    % This plugs in the above to the appropriate part of the B matrix:
	B([1 1+find(designmatrix(i,:)==1)],i)=b;
    % This yields a matrix of Y-Yhat (that is, residuals):
	E=[E (y-x*b)];
end

% Sum of Squares, residual:
SS=(E'*E);
% ML Covariance matrix:
sigma=SS/n;

% The next two lines divide each row and column by the sqrt of the diagonal
% elements to give the correlation matrix.
c=diag(diag(SS).^-.5);
cormatrix=c*SS*c;
MSE=sum(sum(sigma));
lnlike=-n.*(p./2).*log(2.*pi)-(n./2).*log(det(sigma))-n.*p./2;

% parameter count = # non-zero parameters in A, B,C, and some of sigma
par=3 + sum(sum(designmatrix)) + p.*(p+1)./2;

% AIC and BIC
AIC=-2*lnlike + 2*par;

MSE
sigma
cormatrix
AIC

%%%%%%%%%%%%%%%%%% bootstrap CLS estimates %%%%%%%%%%%%%%%%%
Btrue=B';
Etrue=E';
Xtrue=X';
Ytrue=Y';
treat=w(:,10);

'bootstrapped values'

nreps=2000;
alpha=.05;
bootlist=[];
for rep=1:nreps,

	%construct data set
	
	%%%%%%%% NB: The bootstrapping keeps the sets of 3 residuals for
	%%%%%%%% juveniles, adults and infection together, since there is a
	%%%%%%%% correlation between the residuals.  To keep this correlation
	%%%%%%%% in the bootstrap data sets, the residuals are randomly
	%%%%%%%% selected (with replacement) as sets of 3.

	E=Etrue(:,ceil(n*rand(1,n)));
		
	w=[];
	for t=1:n,
        %%Making it so that start the bootstrapped dataset with the actual
        %%t=0 data for that bag:
		if mod(t,14)==1
			x=Xtrue(1:3,t);
		end
		xold=x;
		x=Btrue(:,1) + Btrue(:,2:4)*x + Btrue(:,5).*treat(t)+E(:,t);

		w=[w;x' xold'];
    end
    
	Y=w(:,1:3);
	X=[w(:,4:6) treat];
	
	B=zeros(5,3);	
	E=[];
	for i=1:3
		y=Y(:,i);
		x=[ones(length(Y),1) X(:,find(designmatrix(i,:)==1))];

		b=(x'*x)\(x'*y);
		B([1 1+find(designmatrix(i,:)==1)],i)=b;
		E=[E (y-x*b)];
	end

	SS=(E'*E)/n;
	MSE=sum(sum((E'*E)/n));

	bootlist=[bootlist;B(:)' MSE];
end

[nr nc]=size(bootlist);

%%%%%%%% THIS FINDS THE CONFIDENCE BOUNDS
conf=[];
for i=1:nc,
	sortedrow=sort(bootlist(:,i));
	lower=sortedrow(floor(nr*alpha/2));	
	upper=sortedrow(ceil(nr*(1-alpha/2)));
	conf=[conf;lower mean(sortedrow') upper];
    BSdata(:,i)=sortedrow;
end

save('WGNCLS','BSdata','-ascii');

Btrue
bootmeanB=reshape(conf(1:end-1,2),5,3)'
bootlowerB=reshape(conf(1:end-1,1),5,3)'
bootupperB=reshape(conf(1:end-1,3),5,3)'

bootMSE=conf(end,:)