# [R] conversion from SAS

alessandro carletti alxmilton at yahoo.it
Thu Jul 28 14:11:30 CEST 2005

```Hi, I wonder if anybody could help me in converting
this easy SAS program into R.
(I'm still trying to do that!)

PROC IMPORT OUT= WORK.CHLA_italian
DATAFILE= "C:\Documents and
Settings\carleal\My
Documents\REBECCA\stat\sas\All&nutrients.xls"
DBMS=EXCEL2000 REPLACE;
GETNAMES=YES;
RUN;
data chla_italian;
set chla_italian;
year=year(datepart(date));
month=month(datepart(date));
run;

proc sort data=chla_italian; by station; run;
/* Check bloom for seasonal cycle outliers */
data sort_dataset;
set chla_italian;
chla=chl_a;
dayno=date-mdy(1,1,year)+1;
cos1=cos(2*3.14*dayno/365);
sin1=sin(2*3.14*dayno/365);
cos2=cos(4*3.14*dayno/365);
sin2=sin(4*3.14*dayno/365);
cos3=cos(6*3.14*dayno/365);
sin3=sin(6*3.14*dayno/365);
cos4=cos(8*3.14*dayno/365);
sin4=sin(8*3.14*dayno/365);
bloom=0;
w_chla=1/chla/chla;
run;
ODS listing close;
%macro sort_event(cut_off,last=0);
/*proc glm data=sort_dataset;
class year;
model logchla=year cos1 sin1 cos2 sin2 cos3 sin3
cos4 sin4 /solution;
by station;
where bloom=0;
output out=chla_res predicted=pred student=studres
cookd=cookd rstudent=rstudent u95=u95;
lsmeans year / at (cos1 sin1 cos2 sin2 cos3 sin3
cos4 sin4)=(0 0 0 0 0 0 0 0);
ODS output ParameterEstimates=parmest
LSmeans=lsmeans;
run;*/
proc glm data=sort_dataset;
class year month;
model chla=/solution;
by station;
weight w_chla;
where bloom=0;
output out=chla_res predicted=pred student=studres
cookd=cookd
daynumber<-data\$date-mdy(1,1,y)+1
rstudent=rstudent ucl=ucl lcl=lcl / alpha=0.02;
*   lsmeans year / at (cos1 sin1)=(0 0);
*   ODS output ParameterEstimates=parmest
LSmeans=lsmeans;
run;
quit;
data sort_dataset;
set chla_res;
increase=ucl-pred;
if chla>UCL then bloom=1;
w_chla=1/(50+5*pred*pred);
%if &last=0 %then %do; drop ucl lcl cookd rstudent
studres pred; %end;
run;
%mend sort_event;
ODS listing;
%sort_event(2.0);
%sort_event(2.0);
%sort_event(2.0);
%sort_event(2.0);
%sort_event(2.0);
%sort_event(2.0);
%sort_event(2.0);
%sort_event(2.0);
%sort_event(2.0);
%sort_event(2.0);
%sort_event(2.0);
%sort_event(2.0);
%sort_event(2.0);
%sort_event(2.0);
%sort_event(2.0,last=1);
/* Combine bloom information with all chlorophyll
values */
data chla_sep;
merge sort_dataset(keep=station date bloom pred ucl
lcl) chla_italian (rename=(chl_a=chla));
by station date;
*   lcl=exp(lcl);
*   ucl=exp(ucl);
*   pred=exp(pred);
if bloom=. then bloom=1;
if bloom=0 then chla1=chla; else chla1=.;
if bloom=1 then chla2=chla; else chla2=.;
run;

symbol1 i=none value=plus color=red;
symbol2 i=none value=plus color=green;
symbol3 i=join value=none line=1 color=black;
axis1 logbase=10; axis1;
proc gplot data=chla_sep;
plot chla2*date=1 chla1*date=2 (ucl pred
lcl)*date=3 /overlay vaxis=axis1;
by station;
run;
quit;
proc glm data=chla_sep;
class station year month;
model salinity temperature Transparency__m_
Nitrate__mmol_l_1_ Phosphate__mmol_l_1_
Silicate__mmol_l_1_=bloom month/solution;
by station;
run;
quit;

Thanks

```