https://github.com/cran/bild
Tip revision: 6661b7703cf2271cb13756a6ff2969ef8e1d88f3 authored by M. Helena Goncalves on 03 November 2023, 14:40:02 UTC
version 1.2-1
version 1.2-1
Tip revision: 6661b77
logL_bin0I_aux.R
######## to compute pij
logL.bin0I.aux<- function(parameters, X, data, trace)
{
loglik1 <- function(param, X, y,trace)
{
npar <-as.integer(length(param))
beta<- as.double(param[1:(npar-1)])
bt<- as.double(param[1:(npar-1)])
log.psi<-as.double(0)
omega<-as.double(param[npar])
y<- as.integer(y)
n <- as.integer(length(y))
x<-matrix(as.double(X),nrow=n,ncol=npar-1)
theta<- work<-pij<- prob<-as.double(rep(0,n))
logL <- as.double(0)
eta<-i.fit<- fit<-as.vector(npar-1)
eta<-x%*%beta
##### Estimate bi
m<-glm(as.numeric(y)~offset(eta), family=binomial)
# bi<-coef(m) # comment 27/03/2020
bi<-coef(m)[1]
beta[1]<-beta[1]+bi
y[is.na(y)]<-(-1)
y<- as.integer(y)
results <- .Fortran("mblik1",logL,prob,beta,log.psi,
npar,x,y,theta,work,n,PACKAGE="bild")
i.fit<-x%*%beta
return(list(loglik=results[[1]],pij=results[[2]], fit=i.fit, bi.est=bi))}
if(trace) cat(paste("\t",(format(parameters[length(parameters)], digit=4)), collapse=" "), "\t")
nparam<-length(parameters)
omega1<-as.double(parameters[nparam])
ti.repl<-data[[1]]
cumti.repl<-cumsum(ti.repl)
n.cases<- length(ti.repl)
y<-data[[2]]
counts<-data[[3]]
mt<-fitted<-as.double(rep(0,length(y)))
bi.estimate<-as.double(rep(0,length(n.cases)))
logL1<-as.double(0)
k1<-1
for (i in 1:n.cases)
{
k2<-cumti.repl[i]
z<-loglik1(param=parameters,X=X[k1:k2,], y=y[k1:k2],trace=trace)
mt[k1:k2]<-z$pij
fitted[k1:k2]<-z$fit
bi.estimate[i]<-z$bi.est
k1<-k2+1
}
for (j in 1:length(y))
{
if (mt[j] > 0.999999999999) mt[j]<- 0.999999999999
# else if (mt[j] < 1e-11) mt[j]<- 1e-11
}
return(list(pij=mt, fit=fitted, bi.est=bi.estimate))}