BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.
hajsalmohs0
Fluorite | Level 6

Dear all,

 

I appreciate that if you let me know if there is any module implemeted in SAS for Estimating Mean of a Poisson Distribution Based on Interval censoring? and if yes, could you please name it:) 

 

Kind regards,

Mohsen

1 ACCEPTED SOLUTION

Accepted Solutions
JacobSimonsen
Barite | Level 11

Maybe you can use this little example, where I made the censoring very extreme to verify that it actually works. I generate outcomes from a Poisson distribution with mean 10, and lower censoring at 8 and upper censoring at 12. Using the functions logsdf and logcdf (log of probability to respectively upper and lower tail) inside proc nlmixed I get almost exact the parameter used generate data (beta = log(10)).

 

data abc;
  do i=1 to 10000;
    y=rand('poisson',10);
	output;
  end;
ruN;

/*lower censoring at 3 and upper censoring at 15*/
proc nlmixed data=abc;
  parameters beta 2.3; 
  if y<=8 then l=logcdf('poisson',8,exp(beta));
  else if y>12 then l=logsdf('poisson',12,exp(beta));
  else l=logpdf('poisson',y,exp(beta));
  model y ~ general(l);
run;

View solution in original post

14 REPLIES 14
Norman21
Lapis Lazuli | Level 10

Something like this?

 

http://support.sas.com/documentation/cdl/en/statug/63962/HTML/default/viewer.htm#statug_mcmc_sect062...

 

Norman.

Norman.
SAS 9.4 (TS1M6) X64_10PRO WIN 10.0.17763 Workstation

JacobSimonsen
Barite | Level 11

It sounds like you have time to event data which are interval censored. In that case you can use the new PROC ICPHREG, which assume some some parametric model on the baseline hazard and you can will get estimates of both the baseline hazard rate and the parameters of interes (relative rates).

 

If you have time to event data it is not correct use the term "Poisson distribution". Poisson regression comes into play because the likelihood function will be the same as if the number of events were Poisson distributed, and that is only true if you have exact knowledge of time to event. For interval censored data you can not use Poisson regression because the likelihood doesnt get the form similar to that from Poisson distributed data..

hajsalmohs0
Fluorite | Level 6

Hi Jacob

 

Thanks for your explanation and also the link and sorry for my late reply. I was reading some books and also studing the proc that you suggested. It looks helpful, but still I do not know how to use that. It seems that I have to define everything for possion in this proc. I am working on that.

 

Regards,

Mohsen

JacobSimonsen
Barite | Level 11

I still havent understood what kind of problem you are dealing with. Is it observations from a Poisson distribution (count data) where you only know a lower bound and an upper bound?

Or, is it, as I think, time to event data, where you can not use the Poisson regression technique because of interval censoring?

hajsalmohs0
Fluorite | Level 6

Hey Jakob, 

 

Thanks for your reply. I used R and solved my problem. It has a package with the name of fitdist that gave me the estimation. However, I would also like to get the result from SAS and compare it with R and I really appreciate your support.

 

Please let me explain that. We collected the data from some devices and they are able to count the data within a minutes. For example, we had 10 observation in 12:01, 15 in 12:02, and so forth. I randomly distributed the data (10 obs) in the first 60 sec, and the 15 obs in (60 to 120 sec) and count the interarrival of obs. In this case, I got a result like this:

1, 200

2, 400

3, 150

The first row means I had 200 objects that arrived with the difference of 1 second, 400 objects with two seconds arrival rate, and so on. to do that, I wrote a function to repeat 1, two hundred times, repeat 2, four hundered times, and ...

 

Now, I would like to investigate which distribition they follow. Since some of the devices were faulty, I had some wrong data that needed to be consored, and those are obs<3 sec and obs>500 sec.  I would like to do interval censoring for ❤️ and probably right cens for obs>500. 

 

I hope I could have explained clearly. If you need more information, please let me know.

 

Cheers,

Mohsen

Norman21
Lapis Lazuli | Level 10

You might like to look at this:

 

http://blogs.sas.com/content/iml/2012/04/04/fitting-a-poisson-distribution-to-data-in-sas.html

 

Norman.

Norman.
SAS 9.4 (TS1M6) X64_10PRO WIN 10.0.17763 Workstation

hajsalmohs0
Fluorite | Level 6

Hey Norman, 

 

Thanks for the link, I already studied that. It did not cover the censoring which is my main concern.

 

Many thanks

Mohsen

JacobSimonsen
Barite | Level 11

Maybe you can use this little example, where I made the censoring very extreme to verify that it actually works. I generate outcomes from a Poisson distribution with mean 10, and lower censoring at 8 and upper censoring at 12. Using the functions logsdf and logcdf (log of probability to respectively upper and lower tail) inside proc nlmixed I get almost exact the parameter used generate data (beta = log(10)).

 

data abc;
  do i=1 to 10000;
    y=rand('poisson',10);
	output;
  end;
ruN;

/*lower censoring at 3 and upper censoring at 15*/
proc nlmixed data=abc;
  parameters beta 2.3; 
  if y<=8 then l=logcdf('poisson',8,exp(beta));
  else if y>12 then l=logsdf('poisson',12,exp(beta));
  else l=logpdf('poisson',y,exp(beta));
  model y ~ general(l);
run;
hajsalmohs0
Fluorite | Level 6

Awesome Jacob.

 

I will try that and inform you about the result. I really appreciate that.

 

Cheers,

Mohsen

 

 

hajsalmohs0
Fluorite | Level 6

Hi Jacob,

Thank you for your code. Please correct me if I am making any possible mistake,

You generated the data with mean 10 and then gave the mean to nlmixed because you knew that the mean is 10. For me, since my data contains some corrupted data, I can not get my mean(lambda) and give straight away to nlmixed. I need to do censoring then get the lambda. Does this make sense to you?

Thanks again for your thoughtfulness and time.

hajsalmohs0
Fluorite | Level 6

Jacob, I changed the beta to a random number and it returned the estimated value for me, It just needs a starting value. Awesome, Sincerely appreciate that.

JacobSimonsen
Barite | Level 11
Maybe you can also find a good starting point for your beta value by just ignoring that some observations was censored. I mean, if the lower bound is 3, then just use the number 3 as a normal observation. The starting point is then easily found in with proc genmod.

proc genmod data=mydata;
model count=/dist=poisson;
run;

Then put the estimated beta-value into proc nlmixed.
hajsalmohs0
Fluorite | Level 6

Yes good idea,. Please accept my deepest gratitude. I really appreciate that.

 

Thanks a million, 

Mohsen

sas-innovate-2024.png

Join us for SAS Innovate April 16-19 at the Aria in Las Vegas. Bring the team and save big with our group pricing for a limited time only.

Pre-conference courses and tutorials are filling up fast and are always a sellout. Register today to reserve your seat.

 

Register now!

What is ANOVA?

ANOVA, or Analysis Of Variance, is used to compare the averages or means of two or more populations to better understand how they differ. Watch this tutorial for more.

Find more tutorials on the SAS Users YouTube channel.

Discussion stats
  • 14 replies
  • 2867 views
  • 4 likes
  • 3 in conversation