proc import out=Yellowfin
datafile="/home/grego1/STAT 705/Yellowfin.txt"
replace;
run;

proc sgplot data=yellowfin;
scatter x=age y=length/group=ocean;
xaxis min=0;
yaxis min=0;
run;

proc nlin data=yellowfin plots=all; where ocean='Pacific';
parms Linf=200 to 300 by 50 x0=-20 to 20 by 10 K=0.001 to 0.005 by 0.001;
expo=exp(-K*(age-x0));
mu=Linf*(1-expo); 
model length = mu;
der.Linf=1-expo;
der.K=Linf*(age-x0)*expo;
der.t0=-Linf*K*expo;
output out=fit pred=Lpred;
run;

*nlmixed doesnt accept PLOTS=ALL;
proc nlmixed data=yellowfin; where ocean='Pacific';
*I like that sigma is a parameter;
parms Linf=250 x0=-10 K=0.004 sigma=5;
mu=Linf*(1-exp(-k*(age-x0)));
model length ~ normal(mu,sigma*sigma);
predict Linf*(1-exp(-k*(age-x0))) out=fit;
run;

proc sort data=fit; by age;
run;

proc sgplot data=fit;
scatter x=Age y=Length;
series x=Age y=pred;
run;

proc nlin data=yellowfin plots=all; where ocean='Atlantic';
*Let's keep a range of x0 and lower K slightly based on the scatterplot;
parms Linf=150 x0=-20 to 20 by 10  K=0.003;
expo=exp(-K*(age-x0));
mu=Linf*(1-expo); 
model length = mu;
der.Linf=1-expo;
der.K=Linf*(age-x0)*expo;
der.t0=-Linf*K*expo;
output out=fit pred=Lpred;
run; 

*Not behaving; 
proc nlmixed data=yellowfin; where ocean='Atlantic';
parms Linf=200 x0=70 K=0.004 sigma=5;
mu=Linf*(1-exp(-k*(age-x0)));
model length ~ normal(mu,sigma*sigma);
predict Linf*(1-exp(-k*(age-x0))) out=fit;
run;

*Placeholder for CE 11;
data snake;
input conc rate @@;
datalines;
31.25 53.01 62.5 81.42 125 122.11 250 304.57 500 376.87
1000 414.13 2000 553.46
;
proc sgscatter; plot rate*conc;
run;

proc nlmixed data=snake;
parms b1= b2= sigma=; * let's find these in class;
mu=b1/(1+(b2/conc));
model rate ~ normal(mu,sigma*sigma);
predict b1/(1+(b2/conc)) out=fit;
estimate "mean rate at conc=750" b1/(1+(b2/750));
run;

proc nlmixed data=snake;
parms b1= b2= b3= sigma=; * let's find these in class;
mu=b1/(1+(b2/conc)**b3);
model rate ~ normal(mu,sigma*sigma);
predict b1/(1+(b2/conc)**b3) out=fit;
estimate "mean rate at conc=750" b1/(1+(b2/750)**b3);
run;

proc sgplot data=fit;
scatter x=conc y=rate;
series x=conc y=pred;