Can you specify more detail about what the problem is. I see often problems with log-binomial models, due to the fact that there is no mechanism that ensure that all estimated probabilites is in the interval (0,1) at the maximum-likelihood value. I thiink  an error will occur at the first iteration in the mcmc-algortitm where one of the probabilities is estimated outside the interval. So, maybe a change to logit-scale can solve the problem. ... View more

I am still confused about available fileformats with ODS graphic. And I would therefore like to repeat the original question: Does anyone know how to produce graphics in EPS format with ODS? From the documentation it seems that EPS should work with the ODS LISTING destination, but in practice it doesn't. The error message is "ERROR: Unsupported image file format for output file." Also, EMF format shold be possible with the RTF-destination. Then I get this error message: "ERROR: Trying to open a key that does not exist: CORE\PRINTING\PRINTERS\\DEFAULT SETTINGS" In practice, I therefore often use the old method where I produce the graphics with annotate and gplot. And I can put it to almost whatever format I want with the device option. Jacob ... View more

Thanks Pierre, your solution is very elegant. The very good thing is that with this solution is that the data is never written on disk, which give good performance (for large data of course). Your program can be extended even further if one will have the statistical test, but then it starts to get more difficult. /*access to a function library with matrixalgebra function is required, see attached file. Here I just have the matrixalgebra functions in the work library*/ options cmplib=work.func; Data _null_;   infile cards missover eof=done;   array x{1,2} _temporary_;   array xt{2,1} _temporary_;   array xx_temp{2,2} _temporary_;   array xx{2,2} _temporary_;   array varians{2,2} _temporary_;   array xy{2,1} _temporary_;   array beta{2,1} _temporary_;   Input YEAR ACTUAL;   x[1,1]=1;   x[1,2]=year;   N+1;   yearsum+year;   yearuss+year**2;   actualsum+actual;   actualuss+actual**2;   productsum+actual*year;   call trans(x,xt);   call multiplicer(xt,x,xx_temp);   if _N_=1 then do;     call zero(xx);     call zero(xy);     error=0;   end;   call add(xx,xx_temp,xx);   pred=0;   do i=1 to dim(x,2);     xy[i,1]+x[1,i]*actual;     pred+x[1,i]*actual;   end;   return;   done:   error=(actualuss-(actualsum**2)/N -( (productsum -yearsum*actualsum/N)**2)/(yearuss-yearsum**2/N))/(N-dim(x,2));   call invers(xx,varians);   call multiplicer(varians,xy,beta);   call show(beta);   put @1 'Estimate' @20 'Std Err' @ 40 'P-value';   do i=1 to dim(beta,1);     stderr=sqrt(varians[i,i]*error);     *pvalue=sdf('chisq',(beta[i,1]/stderr)**2,1);     pvalue=2*sdf('t',abs(beta[i,1]/stderr),N-dim(beta,1));     put @1 beta[i,1] @20 stderr @40 pvalue pvalue6.4;   end;   call symputx("Trend", beta[2,1]); cards; 1995 4356.26 1996 8520.93 1997 8161.92 1998 8715.83 1999 4060.38 2000 5446.24 2001 2212.86 2002 4348.00 2003 3803.34 2004 3095.96 2005 8455.64 2006 3553.20 ; run; ... View more

The question is already solved above. But I also like this solution because it runs within one single datastep: Data Test;   infile cards missover ;    Input YEAR   ACTUAL best8.;   retain productsum yearsum actualsum yearuss; if _N_=1 then do;   yearsum=0;   productsum=0;   actualsum=0;   yearuss=0; end;   productsum=productsum+actual*year;   yearsum=yearsum+year;   actualsum=actualsum+actual;   yearuss=yearuss+year**2; if _N_=12 then do;   yearcss=yearuss-(yearsum)**2/_N_;   put productsum= yearsum= actualsum= yearuss= yearcss=;   beta=(productsum-yearsum*actualsum/_N_)/yearcss;   call symputx("Trend", Estimate); end; cards; 1995 4356.26 1996 8520.93 1997 8161.92 1998 8715.83 1999 4060.38 2000 5446.24 2001 2212.86 2002 4348.00 2003 3803.34 2004 3095.96 2005 8455.64 2006 3553.20 ; run; ... View more

Hi, There are several procedures which can estimate a linear trend. One of them are GLM proc glm data=test;   model actual=year; run; in macro form: %macro trend(data,dependent, independent);   ods output parameterestimates=parameterestimates;   proc glm data=&data.;   model &dependent=&independent; quit; %mend; %trend(test,actual,year); Then you can find the trend estimate in the output window and in the data "parameterestimates". You also calculate the number manually, which of course is more cumbersome, but also more fun:-). Then you also will need a variable which contain the product of the two variables: Data Test; infile cards missover; Input YEAR   ACTUAL best8.; product=actual*year; cards; . . . run; proc means data=test uss sum n;   var year product actual;   output out=out sum=yearsum productsum actualsum css=yearcss productcss actualcss n=yearn procuctn actualn; run; *write the number in the log window:; data _NULL_;   set out;   beta=(productsum-yearsum*actualsum/yearn)/yearcss;   put beta; run; ... View more

If there is any biological reason for choosing one model instead of the other then I will choose the biological meaningfull model. Otherwise, I will choose a model which fit the data, and which can produce give meaningfull estimates. I will therefore recommend to do some assesment of the model fit, which is easy done with the assessment statement in proc genmod. Especially, you are interested in assessing the linkfunction, therefore you can add something like this line after the model-statement in the genmod procedure: ASSESS   LINK / nsample=50  nsim=1; There are excellent examples in the sas-documentation of how to interpretate the assesment graphs. You can turn up the number of simulations (nsim and nsample), but start with some small numbers. Unfortunately, the calcuation-time of the assess-statement is O(n^2), so your number of observation must not be very high. If the procedure is still running after an hour, then interrupt the procedure and forget my suggestion until the sas-programmers optimize the assessment method. Jacob ... View more

Hi, I see this problem sometimes. I don't think the error message has anything to do with your stratification on gender. And also not a coding mistake. The problem often occur in binary regression models when the link-function is different from the canonical link-function. The canonical link function in a logistic regression is the logit-function, -not the log function as you specified. There is no mechanism that guarantee that the maximum-likelihood estimate exists within the allowed range of expected values (here the interval (0-1) ), when the log-function is used instead of the logit function (the canonical link). You can see if the problem dissappear when you use the logit function, but then you should also be aware that you estimate log(OR) and not log(RR). This is maybe not what you want, but my disappoint answer is that the log(RR) can not always be estimated with a Maximum likelihood method. It can also be that you can estimate the log(RR) for most of you covariates, and therefore the problem may dissappear when you remove one of you covariates. best, jacob ... View more

Hi, I am not in expert in this, and powercalculations is not something I do often. Though, I think I know what happens, anyone are welcome to correct me if I am wrong. I think you should specifiy "test=diff" instead, because your test is a two-sided test. The equiv_diff allows you to specify two one-sided tests, and therefore you both need to specify the upper alternative and the lower alternative. This is not what you what if I understood you correct. I suggest this code, which was generated by the SAS power and sample size application: proc power;    TwoSampleMeans test=diff       Alpha = 0.05       Sides = 2       MeanDiff = 4       StdDev = 8       Power = 0.8       NTotal = .       NFractional  ; run Best. ... View more

Phreg can't be used to analyse intervalcensored data in general.  Phreg can handle lefttruncated data as specified by fabie10, but lefttruncated data is not the same as intervalcensoring. Though, If the structure of your data is simple enough, then the interval censored data can be analysed as tied data, which can be analyzed with ties=discrete in phreg. In the very newest version there is this proc iclifetest available, which can  anlyse interval-censored data. I haven't tried it, but I don't think it allow fra a random effect. There is also the possibility of using proc nlmixed. There you can specify a any likelyhoodfunction and use the random effects. There is a good example of how to specify a right-censored frailty mdel in the examples in the sas-help section. This method can also be used when data are interval censored. The draw-back of this method is that you need a very detailed knowledge of how your likelihood function is constructed, since you have to specify it yourself. jacob ... View more

I agree. id should be in random, otherwise the model is much too flexible the paramers can not be identified. ... View more

I would like to make a comment more about the model where you put the time to zero at each event. It is a legal model, and easy to handle compared to the model where you use the original timescale. This is so because you avoid the lefttruncations. But, the underlying time is not the same, which it almost no problem if it is acceptable in your context. Though, I see one problem: namely that the time to first event is (perhaps) measured from some arbitrary time, wheareas time to the next events is measured from an event-time. A partly solution can be to stratify or justify for eventnumber (which is legal since as it belong the the past history). Jacob ... View more

There is also this possibility: You can split you data up in intervals, such that it doesn't look like recurrent event data. Then if you have an event, the next record will be left-truncated at the time of the first event. If one more event happens the next record will be left-truncated at the time of the second event, and so on. Notice that the cox-regression model is not violated by the fact that same indiviual contribute with several records, since the intervals for each person will not overlap. the procedure will be something like this, proc phreg data=longtable;    class  groupe id;    model (lefttrunctime time)*Status(0)=groupe  ;    strata id; run; where longtable is a derived dataset, with as many records for each individual as there are events for this individual. There is also the variant where you start the time from zero at each event, and use the time since last event as underlying time. This is not exactly the same model. (you can read about this difference in ("The statistical analysis of recurrent events" by richard J. cook and Jerald F lawless). I hope this make sense. ... View more

I completely agree with Marc's reply. The event should be at the observation where the event took place. The stop-time at that observation should be time of event, and this timeinterval has to be the last one for that individual. I think it was also what you where advised to do by your statistician (the event should be at the last interval for an individual). Intervals after the event took place should be deleted. Jacob ... View more