### Likelihoods ### # rm(list = ls()) # Clean the workspace ### Define function 'expit' for back transformation from logit scale expit <- function(x) 1/(1 + exp(-x)) ## SITE CONFIRMATION DESIGN # Multiple-Methods # MultipleMethod <- function(beta,data,occ){ psi = expit(beta[1]) p11 = expit(beta[2]) p10 = expit(beta[3]) r11 = expit(beta[4]) y = data[,1] w = data[,2] Lik = psi*dbinom(y,occ[1],p11)*dbinom(w,occ[2],r11) + (1-psi)*dbinom(y,occ[1],p10)*dbinom(w,occ[2],0) nLL = sum( -log(Lik) ) return(nLL) } # Multiple-States # MultipleState <- function(beta,data,occ){ psi = expit(beta[1]) p11 = expit(beta[2]) p10 = expit(beta[3]) b = expit(beta[4]) pi.z0 <- c( 1-p10, p10, 0 ) ; sum(pi.z0) pi.z1 <- c( 1-p11 , (1-b)*p11 , b*p11 ) ; sum(pi.z1) y = data Lik = ( psi*apply(t(pi.z1^t(y)), 1, prod) ) + ( (1-psi)*apply(t(pi.z0^t(y)), 1, prod) ) nLL = sum( -log( Lik ) ) return(nLL) } ## CALIBRATION DESIGN Calibration <- function(beta,dataField,occ,dataLab){ psi = expit(beta[1]) p11 = expit(beta[2]) p10 = expit(beta[3]) y = dataField x = dataLab[1,1] n1 = dataLab[1,2] w = dataLab[2,1] n0 = dataLab[2,2] nLL = ( -log(dbinom(x,n1,p11)) ) + ( -log(dbinom(w,n0,p10)) ) + sum( -log( psi*dbinom(y,occ,p11) + (1-psi)*dbinom(y,occ,p10) ) ) return(nLL) } ## OBSERVATION CONFIRMATION DESIGN ObsConf <- function(beta,data,occ){ psi = expit(beta[1]) s1 = expit(beta[2]) s0 = expit(beta[3]) p11 <- s1 + s0 - (s1*s0) p10 <- s0 pi.z1 <- c( (1-s1-s0+s1*s0) , (s0-s1*s0) , (s1-s1*s0), (s1*s0) ) pi.z0 <- c( (1-s0) , s0 , 0 , 0 ) y = data[[1]] w = data[[2]] nLL = sum( -log( psi*(p11^y)*((1-p11)^(occ-y)) + (1-psi)*(p10^y)*((1-p10)^(occ-y)) ) ) + sum( -log( psi*apply(t(pi.z1^t(w)), 1, prod) + (1-psi)*apply(t(pi.z0^t(w)), 1, prod) ) ) return(nLL) } ### Load Data ### load("DataSiteConf-MM.RData") load("DataSiteConf-MS.RData") load("DataCal.RData") load("DataObsConf.RData") ### Optimization of the Likelihood functions ### beta.1 <- c(1,0,-1,-1) # starting values -- set 1 (4 parameters) beta.2 <- c(1,0,-1) # starting values -- set 2 (3 parameters) ### SITE CONFIRMATION DESIGN ## Multiple-Methods out.SiteConf.MM = optim(beta.1, MultipleMethod, data = data.SiteConf.MM, occ = occ.SiteConf.MM, hessian = T) # MLE, SE and true parameter value round(expit(out.SiteConf.MM$par), 3) ## estimates: parameters are psi, p11, p10, r11 round(sqrt(diag(solve(out.SiteConf.MM$hessian))), 3) ## SE: parameters are psi, p11, p10, r11 TrueValues.SiteConf.MM ## Multiple-States out.SiteConf.MS = optim(beta.1, MultipleState, data = data.SiteConf.MS, hessian = T) # MLE and SE round(expit(out.SiteConf.MS$par), 3) ## estimates: parameters are psi, p11, p10, b round(sqrt(diag(solve(out.SiteConf.MS$hessian))), 3) ## SE: parameters are psi, p11, p10, b TrueValues.SiteConf.MS ### CALIBRATION DESIGN out.Cal = optim(beta.2, Calibration, dataField = dataField.Cal, occ = occ.Cal, dataLab = dataLab.Cal, hessian = T) # MLE and SE round(expit(out.Cal$par), 3) ## estimates: parameters are psi, p11, p10 round(sqrt(diag(solve(out.Cal$hessian))), 3) ## SE: parameters are psi, p11, p10 TrueValues.Cal ### OBSERVATION CONFIRMATION DESIGN out.ObsConf = optim(beta.2, ObsConf, data = data.ObsConf, occ = occ.ObsConf, hessian = T) # MLE and SE round(expit(out.ObsConf$par), 3) ## estimates: parameters are psi, s1, s0 round(sqrt(diag(solve(out.ObsConf$hessian))), 3) ## SE: parameters are psi, s1, s0 TrueValues.ObsConf #### END ####