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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
