Re: [R] conversion from SAS

From: Frank E Harrell Jr <f.harrell_at_vanderbilt.edu>
Date: Fri 29 Jul 2005 - 00:13:00 EST

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

Converting that program into R will be very feasible and the solution will be far more elegant than SAS. But I think you are expecting other people to do your work.

Frank

>
> 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
>

-- 
Frank E Harrell Jr   Professor and Chair           School of Medicine
                      Department of Biostatistics   Vanderbilt University

______________________________________________
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 Fri Jul 29 00:32:30 2005

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