[R] Maximum likelihood for censored geometric distribution

A.D.Higginson A.D.Higginson at bristol.ac.uk
Wed Nov 16 13:44:45 CET 2011

Hi all,

I need to check for a difference between treatment groups in the 
parameter of the geometric distribution, but with a cut-off (i.e. right 
censored). In my experiment I stimulated animals to see whether I got a 
response, and stopped stimulating if the animal responded OR if I had 
stimulated 10 times. Since the response could only be to a stimulation, 
the distribution of response times is discrete (so I can't use survival 
analysis). I am interested in assessing at what stimulation number the 
animals responded, and whether this differed between groups. Assuming 
that the animals have a fixed probability (p) of responding to each 
stimulation the distribution of responses will follow a geometric with 
parameter p, but censored at 10 stimulations. I hope to estimate p for 
my two treatment groups (p1 and p2), and then see whether the treatment 
affects p (i.e. whether p1 and p2 differ).

# A subsample of the data, where 11 represent no response to 10 stimulations
# create the geometric likelihood function
geometric.loglikelihood <- function(pi, y, n) {
     loglikelihood <- n*log(pi)-n*log(1-pi)+log(1-pi)*sum(y)
# optimise the likelihood function
# print parameter value at maximum likelihood

The above does give me a reasonable estime of p in par (except I do get 
report of NaNs produced), but counts 11 as stimulation numbers rather 
than censored, so is incorrect. I can't find any information on how to 
express the cut-off distribution (piecewise?).

Also, I plan to calculate the SE for p1 and p2, and then do a t.test to 
see if they differ. Will this be correct?

Thanks in advance,


More information about the R-help mailing list