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