inst/cointack2.R

In skgallagher/TBornotTB: Analyze Tuberculosis Patients in Maryland

# We generate Yi ~ logit(b0+b1 Xi) for i=1,2 Xi~be(1/2)
#
# we get mles of b0,b1 from this data ,
#
# we then scramble up the Xs and derive the correct likelihood
# e.g. P( Y+ = 1 | X+ = 1)  = p0(1-p10 + p1*(1-p0)
#
# we also do the  monte carlo likelihood which fails, for now
#

llFIB<-function(parm){
    b0<-parm[1];b1<-parm[2]
    py<-1/(1+exp(-b0-b1*X))
    - sum((as.vector(Y)*log(py) + as.vector(1-Y)*log(1-py)))
}


llSIB<-function(parm){
    b0<-parm[1];b1<-parm[2]
    py<-1/(1+exp(-b0-b1*X))
    p0<-(1- 1/(1+exp(-b0-b1*0)))*(1- 1/(1+exp(-b0-b1*1))) +
        (1- 1/(1+exp(-b0-b1*0)))*(1- 1/(1+exp(-b0-b1*1)))
    p1<-(1- 1/(1+exp(-b0-b1*0)))*1/(1+exp(-b0-b1*1))      +
        (1- 1/(1+exp(-b0-b1*1)))*1/(1+exp(-b0-b1*0))
    p2<-(1/(1+exp(-b0-b1*0)))   *1/(1+exp(-b0-b1*1))      +
        (1/(1+exp(-b0-b1*1)))   *1/(1+exp(-b0-b1*0))
    - sum(   Y[(SX!=1)]*log(py[(SX!=1)]) + (1-Y[(SX!=1)])*log(1-py[(SX!=1)])) -
        sum(S0)*log(p0)-sum(S1)*log(p1)-sum(S2)*log(p2)
}


llFIGMC<-function(parm){
    b0<-parm[1];b1<-parm[2]
    XVEC<-NULL
    ll<-0
    for(j in 1:nmc){
      XVEC<-MCXVEC[j,]
      py<-1/(1+exp(-b0-b1*XVEC))
      ll<-  - sum(  YY*log(py) + (1-YY)*log(1-py))   + ll
    }
ll
}

n<-2000; MM<-2
p<-.5; B0<- -1; B1<- 2
Y<-matrix(rep(-9,n*MM),nrow=n)
X<-Y


## Bernoulii sampling
for(i in 1:n){
   X[i,]<-rbinom(MM,1,p)
   Y[i,]<-rbinom(MM,1,1/(1+exp(-B0-B1*X[i,])))
}
parm<-c(0,0)
optim(parm,llFIB)

##Proper Likelihood with scrambled Xs
SX<-apply(X,1,sum)
SY<-apply(Y,1,sum)
S0<-( (SX==1)&(SY==0) )
S1<-( (SX==1)&(SY==1) )
S2<-( (SX==1)&(SY==2) )
parm<-c(0,0)
optim(parm,llSIB)


#
#Monte Carlo up the X vector.   Do it nmc times.   Then maximiaxe log-likelihood summed over
#these nmc realizations.
#
nmc<-100
MCXVEC<-matrix(rep(0,nmc*MM*n),nrow=nmc)
for(j in 1:nmc){
    DUM<-NULL;YY<-NULL
    for(i in 1:n){
#     DUM<-c(DUM,sample( X[i,]))
        DUM<-c(DUM, X[i,] )
        YY<-c(YY,Y[i,])
  }
    MCXVEC[j,]<-DUM
    }


## Now maximize over scrambled X
parm<-c(0,0)
optim(parm,llFIGMC)