* Activate the next 2 lines to avoid the problem of having the SAS Log Window becoming full; * PROC PRINTTO LOG='C:\Windows\Temp\LOGFILE.Txt'; * RUN; * beginning of the macro program 'RPT'; %MACRO RPT (reps=,numsites=,numnests=,trueB0=,trueB1=,s2site=); /* The macro parameters are:; # of reps, # of sites, # of nests, true Beta values, and true random effects (s**2); This macro is for a single year but can be modified to allow for multiple years. */ %Do L=1 %to &reps; * how many replications are to be done; DATA nests; /*** set parameters ***/ retain seed_site (2); * This can be changed to get different simulations; ssite=sqrt(&s2site); * Calculate the std dev for the random effect of site; site=0; * Initializes the counter for site number; /*** end of parameters ***/ NestID=0; DO s=1 TO &numsites; * number of sites per dataset; site=site+1; call rannor(seed_site,x2); * x2 is used to store the random #, which is ~N(0,1); siterndm=0+ssite*x2; PpnGrass=ranuni(0); * site-level covariate condition is uniformly distributed; * Other site-level covariates can be added here; /* generate basic properties of each nest on the site */ DO n=1 TO &numnests; * number of nests per site; init=round((ranuni(0)*59)+1); * initiation dates are uniformly distributed; agestart=round((ranuni(0)*22)+1); * age of nest when found is uniformly distributed; NestID=NestID+1; * Assigns a unique NestID to each nest in sample; * Other nest-site level covariates can be added here; output; end; end; run; /* Generate interval lengths and fates for each nest */ data intervals; set nests; do j=1 to &numnests; * if file is modified to allow multiple years, also modify this line; interval=0; * initial condition; datestart=init+agestart; * day in nesting season when the nest is found; fate=1; * initial condition; do until (check=0); * keep checking nest until nest fails or nest attempt is successful; * check is altered below with x and z; interval=interval+1; a=ranuni(0); if a <0.25 then intleng=5; if 0.25 <= a < 0.5 then intleng=6; if 0.5 <= a < 0.75 then intleng=7; if a >= 0.75 then intleng=8; drop a; sum=agestart+intleng; if sum>35 then intleng=(35-agestart); * all nests complete at age 35 days in this example; do k=0 TO intleng-1; * if i=0 then Ob=1; * this line is used in observer-effects modeling; * else Ob=0; * this line is used in observer-effects modeling; logit=(&trueB0+&trueB1*PpnGrass+siterndm); dsr=exp(logit)/(1+exp(logit)); fate=fate*ranbin(0,1,dsr); end; ageend=agestart+intleng; dateend=datestart+intleng; drop j k sum ; output ; agestart=ageend; datestart=dateend; if agestart>=35 then x=0; else x=1; if fate=0 then z=0; else z=1; check=x*z; if check=0 then return; end; end; run; * Run the Constant DSR model; title 'Constant DSR model'; Proc Nlmixed data=intervals tech=quanew method=gauss maxiter=1000; parms B0=0; p=1; do i=0 TO intleng-1; if i=0 then Ob=1; else Ob=0; logit=B0; p=p*(exp(logit)/(1+exp(logit))); end; model fate~binomial(1,p); ods output ParameterEstimates=B0Hat0 (where=(Parameter='B0') rename=(Estimate=Est_B0 StandardError=SE_B0 Lower=LCI_B0 Upper=UCI_B0)); RUN; data ests0; set B0Hat0; B0True=&TrueB0; if LCI_B0 <= &TrueB0 <= UCI_B0 then CovB0=1; else covB0=0; keep B0True Est_B0 SE_B0 LCI_B0 UCI_B0 CovB0; run; data _Null_; set ests0; proc append base=work.Model_0_Constant_DSR; run; * Run the B0+B1*PpnGrass model; title 'B0+B1*PpnGrass model'; Proc Nlmixed data=intervals tech=quanew method=gauss maxiter=1000; parms B0=0, B1=0; p=1; do i=0 TO intleng-1; if i=0 then Ob=1; else Ob=0; logit=B0+B1*PpnGrass; p=p*(exp(logit)/(1+exp(logit))); end; model fate~binomial(1,p); ods output ParameterEstimates=B0Hat1 (where=(Parameter='B0') rename=(Estimate=Est_B0 StandardError=SE_B0 Lower=LCI_B0 Upper=UCI_B0)); ods output ParameterEstimates=B1Hat1 (where=(Parameter='B1') rename=(Estimate=Est_B1 StandardError=SE_B1 Lower=LCI_B1 Upper=UCI_B1)); RUN; data ests1; set B0Hat1; set B1Hat1; B0True=&TrueB0; B1True=&TrueB1; if LCI_B0 <= &TrueB0 <= UCI_B0 then CovB0=1; else covB0=0; if LCI_B1 <= &TrueB1 <= UCI_B1 then CovB1=1; else covB1=0; keep B0True Est_B0 SE_B0 LCI_B0 UCI_B0 CovB0 B1True Est_B1 SE_B1 LCI_B1 UCI_B1 CovB1; run; data _Null_; set ests1; proc append base=work.Model_1_Grass; run; * Run the B0+random-effect-of-site model; title 'B0+random-effect-of-site model'; Proc sort data=intervals; * data must be sorted by the random effect variable; by site; run; Proc Nlmixed data=intervals tech=quanew method=gauss maxiter=1000; parms B0=4, vsite=.4; p=1; do i=0 TO intleng-1; if i=0 then Ob=1; else Ob=0; logit=B0+u; p=p*(exp(logit)/(1+exp(logit))); end; model fate~binomial(1,p); random u~normal(0,vsite) subject=site; ods output ParameterEstimates=B0Hat2 (where=(Parameter='B0') rename=(Estimate=Est_B0 StandardError=SE_B0 Lower=LCI_B0 Upper=UCI_B0)); ods output ParameterEstimates=vsiteHat2 (where=(Parameter='vsite') rename=(Estimate=Est_vsite StandardError=SE_vsite Lower=LCI_vsite Upper=UCI_vsite)); RUN; data ests2; merge B0Hat2 vsiteHat2; B0True=&TrueB0; s2siteTrue=&s2site; if LCI_B0 <= &TrueB0 <= UCI_B0 then CovB0=1; else CovB0=0; if LCI_vsite <= s2siteTrue <= UCI_vsite then Covvsite=1; else Covvsite=0; keep B0True Est_B0 SE_B0 LCI_B0 UCI_B0 CovB0 s2siteTrue Est_vsite SE_vsite LCI_vsite UCI_vsite Covvsite ; run; data _Null_; set ests2; proc append base=work.Model_2_v_Site; run; * Run the (B0+random-effect-of-site)+B1*PpnGrass model; title '(B0+random-effect-of-site)+B1*PpnGrass model'; Proc sort data=intervals; * data must be sorted by the random effect variable; by site; run; Proc Nlmixed data=intervals tech=quanew method=gauss maxiter=1000; parms B0=3.5, B1=0.09, vsite=.3; p=1; do i=0 TO intleng-1; if i=0 then Ob=1; else Ob=0; logit=(B0+u)+B1*PpnGrass; p=p*(exp(logit)/(1+exp(logit))); end; model fate~binomial(1,p); random u~normal(0,vsite) subject=site; ods output ParameterEstimates=B0Hat3 (where=(Parameter='B0') rename=(Estimate=Est_B0 StandardError=SE_B0 Lower=LCI_B0 Upper=UCI_B0)); ods output ParameterEstimates=B1Hat3 (where=(Parameter='B1') rename=(Estimate=Est_B1 StandardError=SE_B1 Lower=LCI_B1 Upper=UCI_B1)); ods output ParameterEstimates=vsiteHat3 (where=(Parameter='vsite') rename=(Estimate=Est_vsite StandardError=SE_vsite Lower=LCI_vsite Upper=UCI_vsite)); RUN; data ests3; merge B0Hat3 B1Hat3 vsiteHat3; B0True=&TrueB0; B1True=&TrueB1; s2siteTrue=&s2site; if LCI_B0 <= &TrueB0 <= UCI_B0 then CovB0=1; else CovB0=0; if LCI_B1 <= &TrueB1 <= UCI_B1 then CovB1=1; else covB1=0; if LCI_vsite <= s2siteTrue <= UCI_vsite then Covvsite=1; else Covvsite=0; keep B0True Est_B0 SE_B0 LCI_B0 UCI_B0 CovB0 B1True Est_B1 SE_B1 LCI_B1 UCI_B1 CovB1 s2siteTrue Est_vsite SE_vsite LCI_vsite UCI_vsite Covvsite; run; data _Null_; set ests3; proc append base=work.Model_3_Grass_and_Site; run; %End; %MEND RPT ; %RPT (reps=5,numsites=20,numnests=15,trueB0=2.75,trueB1=0.40,s2site=0.25); * The number of reps is set low so that all output can be accomodated in the log window. However, when large numbers of reps are run it may be necessary to activate the 1st 2 lines of the file. The output files are stored in the work directory; Run;