data correctedsurvivalcounty;
infile "Bioassay - corrected survival data.csv" dlm=',' firstobs=2;
input pop $ year fieldtype $ corrected;
run;
proc print;
run;

data correctedsurvivalcounty2;
	set correctedsurvivalcounty;
	if corrected > 1 then corrected=1;
	asr_corrected=arsin(sqrt(corrected));
	run;
proc print data=correctedsurvivalcounty2;
run;

proc sort;
by year fieldtype pop;
run;
ods output summary=with_stackods;
proc means data=correctedsurvivalcounty stackodsoutput n mean stderr;
class year fieldtype pop;
var corrected;
run;

/* no difference in corrected survival between the two field types */

proc glm data=correctedsurvivalcounty2 plots(only)=all;
class fieldtype year;
model asr_corrected=fieldtype|year;
lsmeans fieldtype;
lsmeans year;
where fieldtype not eq "Control";
run;

proc mixed data=correctedsurvivalcounty2;
class year pop;
model arsincorrected=pop;
lsmeans pop;
contrast "(2015) PC5 vs controls" pop 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0;
contrast "(2016) NPC6 vs controls" pop 0 0 -2 0 0 0 0 0 0 0 0 0 1 0 1 0 0;
contrast "(2016) PC7 vs controls" pop 0 0 0 0 0 -2 0 0 0 0 0 0 1 0 1 0 0;
contrast "(2016) PC11 vs controls" pop 0 -2 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0;
contrast "(2017) NPC11 vs controls" pop 0 0 0 0 0 0 0 0 -4 0 0 1 0 1 0 1 1;
contrast "(2017) NPC12 vs controls" pop 0 0 0 0 0 0 0 0 0 -4 0 1 0 1 0 1 1;
contrast "(2017) NPC13 vs controls" pop 0 0 0 0 0 0 0 0 0 0 -4 1 0 1 0 1 1;
contrast "(2017) NPC14 vs controls" pop -4 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 1;
contrast "(2017) PC19 vs controls" pop 0 0 0 0 0 0 -4 0 0 0 0 1 0 1 0 1 1;
contrast "(2017) PC20 vs controls" pop 0 0 0 0 0 0 0 -4 0 0 0 1 0 1 0 1 1;
run;



proc glm data=correctedsurvivalcounty2;
class pop;
model asr_corrected=pop;
contrast "PC5 vs control" pop -1 1;
where year=2015;
run;

proc glm data=correctedsurvivalcounty2;
class year pop;
model asr_corrected=pop;
lsmeans pop;
contrast "PC11 vs controls" pop -2 0 0 1 1;
contrast "NPC6 vs controls" pop 0 -2 0 1 1;
contrast "PC7 vs controls" pop 0 0 -2 1 1;
where year=2016;
run;

proc glm data=correctedsurvivalcounty2;
class year pop;
model asr_corrected=pop;
lsmeans pop;
contrast "NPC14 vs controls" pop -4 0 0 0 0 0 1 1 1 1;
contrast "PC19 vs controls" pop 0 -4 0 0 0 0 1 1 1 1;
contrast "PC20 vs controls" pop 0 0 -4 0 0 0 1 1 1 1;
contrast "NPC11 vs controls" pop 0 0 0 -4 0 0 1 1 1 1;
contrast "NPC12 vs controls" pop 0 0 0 0 -4 0 1 1 1 1;
contrast "NPC13 vs controls" pop 0 0 0 0 0 -4 1 1 1 1;
where year=2017;
run;