[R] conversion from SAS

[alessandro carletti]

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