[R] Maximum likelihood for censored geometric distribution

[A.D.Higginson]

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
Nstim=c(1,1,1,1,1,1,2,2,3,4,5,6,8,11,11,11,11)
# create the geometric likelihood function
geometric.loglikelihood <- function(pi, y, n) {
     loglikelihood <- n*log(pi)-n*log(1-pi)+log(1-pi)*sum(y)
     return(loglikelihood)
}
# optimise the likelihood function
test<-optim(c(0.5),geometric.loglikelihood,method="BFGS",hessian=TRUE,control=list(fnscale=-1),y=Nstim,n=length(Nstim))
# print parameter value at maximum likelihood
print(test)

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,

Andy