[R] conversion from SAS

[Frank E Harrell Jr]

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