# Integrand for cells with score s for a standard normal density function normIntFunc=function (u,b,sc,sig) { tb=length(b) odds=exp(sig*u*sc) for (i in 1:tb) { odds=odds/(1+exp(sig*u+b[i])) } sam = odds * dnorm(u) return(sam) } #Skip this in latest version # Integrand for empty cell for a standard normal density function pi0normIntFunc=function (u,b,sig) { tb=length(b) odds=1 for (i in 1:tb) { odds=odds/(1+exp(sig*u+b[i])) } sam = odds * dnorm(u) return(sam) } #Conditional likelihood; changed from computing 1-pi_0 to accumulating all the pi's in I_0 logL=function (parm,scor,indx,y) { S=length(parm) T=S-1 ny=length(y) beta=parm[1:T] sig=parm[S] s1=0 piI0=NULL for(i in 1:ny){ w1=integrate(normIntFunc,-4,4,b=beta,sc=scor[i],sig=sig)$value piI0[i]=w1*as.numeric(exp(indx[i,]%*%beta)) } piI0=piI0/sum(piI0) logL=sum(y*log(piI0)) return(-logL) } #trial computation of logL y=c(3,6,0,5,1,0,0,3,2,3,0,0,1,0,0,4,2,3,1,0,1,0,0,1,0,0,0,0,0,0,0,4,1,1,1,2,0,2,0,4,0,3,0,1,0,2,0,2,0,1,0,1,0,1,0,1,1,1,0,0,0,1,2) ny=length(y) T=log2(ny+1) indx=NULL #Construct response vector for(i in 1:T){ indx=cbind(indx,rep(rep(c(0,1),each=2^(i-1)),2^(T-i))) } #Remove the unobserved cell indx=indx[-1,] #Compute score for each cell scor=apply(indx,1,sum) #Column marginals ymarg=apply(indx*y,2,sum) sigmaHat = 1.0 betaHat=-2*(ymarg-ymarg[T])/max(abs(ymarg-ymarg[T])) parmHat = c(betaHat, sigmaHat) #Score group totals ns=tapply(y,scor,sum) logL(parmHat,scor,indx,y) ## non-linear minimization ## soln=nlminb(start=parmHat,logL,lower=c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,0),scor=scor,indx=indx,y=y) ############# Compute Nhat################################ parmhat=soln$par pi0Hat=integrate(pi0normIntFunc,-4,4,b=parmhat[1:6],sig=parmhat[7])$value Nhat=sum(ns)/(1-pi0Hat)