[R] conversion from SAS

From: alessandro carletti <alxmilton_at_yahoo.it>
Date: Thu 28 Jul 2005 - 22:11:30 EST


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



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Thu Jul 28 22:18:30 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:34:06 EST