MRBC_EMV

Conteúdo do arquivo

EMV.RBC <- function(dados, fun.p, fun.rho, chutes, erro = 1e-08, digt = 3, conf = 0.95){

BETAS = as.matrix(chutes)
y <- dados[,1]
n <- dados[,2]
v <- dados[,3]
v[which(v==0)]=0.0000001
X <- dados[,4:ncol(dados)]
m = nrow(X)

IA <- numeric(m)
IB <- numeric(m)
IA[y==0] = 1
IB[y==n] = 1

loglik <- function(beta){
preditor = X %*% beta[-1]
if(fun.p==1){ p = exp(preditor)/(1+exp(preditor)) }
if(fun.p==2){ p = 1 - exp(-exp(preditor)) }
if(fun.p==3){ p = exp(-exp(-preditor)) }
if(fun.p==4){ p = pnorm(preditor) }
if(fun.rho==1) { rho = exp(-beta[1]*v) }
if(fun.rho==2) { rho = beta[1]^v }
if(fun.rho==3) { rho = exp(-(beta[1]*v)^2) }
ll <- sum(log(IA*((1-p)^n*(1-rho)+(1-p)*rho)+IB*(p^n*(1-rho)+p*rho)+(1-IA-IB)*(choose(n,y)*p^y*(1-p)^(n-y)*(1-rho))))
return(ll)
}

res <-  suppressWarnings(optim(chutes, loglik, method="BFGS", hessian=T, control=list(fnscale=-1)))
EMV <- res$par
varbeta <- diag(solve(-res$hessian))

confianca <- c((1-conf)/2, 1-(1-conf)/2) 

resumo <- matrix(ncol=3, nrow=length(EMV), 0)

resumo[,1] <- EMV
resumo[,2] <- EMV + qnorm((1-conf)/2)*sqrt(varbeta)
resumo[,3] <- EMV + qnorm(1-(1-conf)/2)*sqrt(varbeta)

colnames(resumo) <- paste(c("EMV", "ICi", "ICs"))
rownames(resumo) <- paste(c("gamma", replicate((length(EMV)-1),"betas")))

return(round(resumo,digt))
}


### Example:

# Specify as below (do not run):

dados <- cbind(y,n,v,1,x1,x2)
        # y: response variable
        # n: number of individuals within each cluster
        # v: value of the individuals' function
        # 1, x1, x2: list of clusters' covariates

fun.p <- 1 - logit link function
             2 - complementary log-log link function
             3 - log-log link function 
             4 - probit link function
        
fun.rho <- 1 - exponential correlation structure 
                2 - continuous AR correlation structure
                3 - Gaussian correlation structure

chutes: list with initial values of parameters

# Use (run):

EMV.RBC(dados=cbind(y,n,v,1,x1,x2), fun.p=2, fun.rho=2, chutes=c(0.5,0,0,0))

# Output:

- EMV: vector of maximum likelihood estimates
- Asymptotic confidence intervals