data turtle2011;
infile 'C:\Users\abronikolab\Downloads\Age_Immune_Data.csv' dlm="," firstobs=2;
informat date $10. DBid $17. id $11.;
input date $ jdays season $ DBid $ id $ Sex $ CCL CCW CL CW PL PW Nesting $ iCapYear Age KnownAge $ 
	ClutchSize Agglut Lysis rawBK Bknoneg;


	/*Remove 2 outlier values for Agglut*/
	If id = "8.9-12" then Agglut = ".";
	If id = "2.3.10.12-1" then Agglut = ".";

If (sex = "Male") then zAGE = ((Age - 8.36)/3.26);
If (sex = "Female") then zAGE = ((Age-14.01)/6.39);
If (sex = "Juv") then zAGE = ((Age - 2.22)/.76);
run;

proc print data = turtle2011;
run;


/*Look at distribution of dependent variables*/
Proc univariate;
var Agglut Lysis Bknoneg;
Histogram;
qqplot /normal(mu=est color=red);
run; 
/*Bknoneg not normal, transformation does not make it normal, so use Bknoneg*/

proc corr;
var Agglut Lysis Bknoneg;
run;
/*No immune variables are correlated*/


/*All Individuals*/

/*Agglutination*/
/*Full Model*/
proc mixed data = turtle2011;
class Sex;
model Agglut = zAGE|Sex|jdays @2;
lsmeans Sex / pdiff;
run;
/*Final Model*/
proc mixed data = turtle2011;
class Sex;
model Agglut = zAGE|Sex|jdays @1;
lsmeans Sex / pdiff;
run;

/*Lysis*/
/*Full Model*/
proc mixed data = turtle2011;
class Sex;
model Lysis =  zAGE|Sex|jdays @2;
run;
/*Final Model*/
proc mixed data = turtle2011;
class Sex;
model Lysis =  zAGE|Sex|jdays @1;
lsmeans Sex / pdiff;
run;

/*Bactericidal Competence*/
/*Full Model*/
proc mixed data = turtle2011;
class Sex;
model Bknoneg =  zAGE|Sex|jdays @2;
run;
/*Final Model*/
proc mixed data = turtle2011;
class Sex;
model Bknoneg =  zAGE|Sex|jdays @1;
run;




/*Reproductive Females Only*/

/*Agglutination*/
/*Full Model*/
proc mixed data = turtle2011;
class Sex;
where Sex eq "Female";
model Agglut = ClutchSize|Age|jdays @2;
run;
/*Final Model*/
proc mixed data = turtle2011;
class Sex;
where Sex eq "Female";
model Agglut = ClutchSize|Age|jdays @1;
run;

/*Lysis*/
/*Full Model*/
proc mixed data = turtle2011;
class Sex;
where Sex eq "Female";
model Lysis = ClutchSize|Age|jdays @2;
run;
/*Final Model*/
proc mixed data = turtle2011;
class Sex;
where Sex eq "Female";
model Lysis = ClutchSize Age jdays ClutchSize*Age;
run;

/*Bactericidal Competence*/
/*Full Model*/
proc mixed data = turtle2011;
class Sex;
where Sex eq "Female";
model Bknoneg = ClutchSize|Age|jdays @2;
run;
/*Final Model*/
proc mixed data = turtle2011;
class Sex;
where Sex eq "Female";
model Bknoneg = ClutchSize Age jdays ClutchSize*jdays Age*jdays;
run;








/*Figure 2*/
/*Create residuals for agglutination model for female data only without Age to assess relationship between female age and agglut*/
proc mixed data = turtle2011;
class Sex;
model Agglut = Sex jdays /  outp=zAGEresid residual;
run; 
proc reg data=zAGEresid;
model resid = zAGE;
run;




/*Figure 3*/
proc mixed data = turtle2011;
class Sex ClutchSize;
where Sex eq "Female";
model Lysis = ClutchSize Age jdays ClutchSize*Age;
lsmeans ClutchSize / at Age = 6 pdiff;
lsmeans ClutchSize / at Age = 8 pdiff;
lsmeans ClutchSize / at Age = 13 pdiff;
lsmeans ClutchSize / at Age = 18 pdiff;
lsmeans ClutchSize / at Age = 33 pdiff;
run;


/*Figure 4*/
proc mixed data = turtle2011;
class Sex Season ;
where Sex eq "Female";
model Bknoneg = Season Age ClutchSize Season*ClutchSize Age*Season;
lsmeans Season / at Age = 6 pdiff;
lsmeans Season / at Age = 8 pdiff;
lsmeans Season / at Age = 13 pdiff;
lsmeans Season / at Age = 18 pdiff;
lsmeans Season / at Age = 33 pdiff;
run;