inst/tests/sim-01.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

lvl <- seq(from=0.25, to=0.75, by=0.25)
mu <- -log(-lvl+1)
theta <- log(mu)
## R> mu
## [1] 0.2876821 0.6931472 1.3862944
##
## Testing
if(FALSE){
ppois(q=0, lambda=mu, lower.tail=FALSE) 
## Interpretation of 0.25 
## 75% of cases are missed (i.e., this system captures 25% of cases)
}

## Use this vector of 'mu' as the null values that gives 0.25%, 0.50%, 0.75%
## of capturing performance
mu0 <- c(0.288, 0.694, 1.386)

##
## Models without covariats
## imprecise gamma prior
## (a,b) = { 0.01 <= a <= 10, 0.01 <= b <= 10}
##
## 'alpha' (shape parameter of gamma-poisson parametrization)
## 'alpha' goes to 'infty' (i.e., 'dispersion' goes to '0');
## thus, 'f(y)' goes to Poisson from NegBinom.
##
## meta-parameters 
## 

rm(list=ls())
library(ipeglim)
options(warn=0) # all warning message as errors 2
## options(warn=-1) >> ignore all warnings
## options(warn=0) >> default
## options(warn=1) >> print warnings immediately (make stack empty) 
## options(warn=2) >> make warning behave as error
## 
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/ztpr.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/summary.ztpr.R")
# source("~/Documents/R-ko/countreg/pkg/R/zerotrunc.R")

mu0 <- c(0.288, 0.694, 1.386)   # performance 
N <- c(5e1,1e2,3e2,1e3)         # pop. size
alpha <- c(2, 5, 10, 100)       # shape (gamma) for over dispsersion
# dispersion = c(0.50, 0.20, 0.10, 0.01)
# nb -> pois
lc.lgamma <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), 
  rhs=c(1e-3, -10, 1e-3, -10))  # constrained M0
niters <- 5e2
id <- 0
op <- list()

for(idx.s in alpha){
  for(idx.m in mu0){
    for(idx.n in N){
      for(i in seq_len(niters)){
      id <- id + 1
      
      # data and its summary stats
      tmp <- simulateYX(N=idx.n, Xreg=FALSE, ztrunc=TRUE, 
              param=idx.m, shape=idx.s)
      y <- tmp$y  #$
      
      # classical ZTP model
      fit.ztp <- try(summary(ztpr(formula=y~1, ztrunc=TRUE, dist="poisson"), 
          LM.test=TRUE, HT.est=TRUE))
      if(class(fit.ztp) == "try-error") fit.ztp <- FALSE
      fit.ztnb <- try(summary(ztpr(formula=y~1, ztrunc=TRUE, dist="nbinom"),
          HT.est=TRUE))
      if(class(fit.ztnb) == "try-error") fit.ztnb <- FALSE 
        
      # imprecise ZTP model
      fit <- model(formula=y~0, ztrunc=TRUE, dist="poisson")
      fit <- impose(obj=fit, eqns=lc.lgamma)
      fit <- update(obj=fit, method="LA", priorType="lgamma")
      fit <- summary(fit, HT.est=TRUE)
      
      rvals <- list(d0=tmp, ztp=fit.ztp, ztnb=fit.ztnb, fit=fit)
      rm(tmp, y, fit.ztp, fit.ztnb, fit)
      op[[id]] <- rvals
      print(paste("id=", id, "; alpha=", idx.s,": mu0=", idx.m, ": N=", idx.n, ": niter=", 
              i, sep=""))
      }
    }
  }
}
warnings()

## Potential errors catched by 'try()'
## 
## Error in optim(par = init, fn = fnllik, gr = grllik, method = method,  : 
##  L-BFGS-B needs finite values of 'fn'
## Error in solve.default(as.matrix(fit$hessian)) : 
##  Lapack routine dgesv: system is exactly singular: U[1,1] = 0

save(op, 
  file="~/Documents/PhD/ipeglim/thesis/chapters/ch04/sim-01-lgamma.RData")
op1 <- do.call(rbind, op)
# lapply(op, with, ztp$N)

if(FALSE){

    # z <- range(sfit$est)
    # zlim <- c(z[1]-0.05, z[2]+0.05)
    # plot(sfit, zlim=zlim)
    # cmf1 <- setGrid(cmf, len=10)
    # fit1 <- update(obj=cmf1, method="LA", priorType="lgamma")
    # sfit1 <- summary(fit1, HT.est=TRUE)
    # z1 <- range(sfit1$est)
    # zlim1 <- c(z1[1]-0.05, z1[2]+0.05)
    # plot(sfit1, zlim=zlim1, xlim=c(-1,11), ylim=c(-1,11), clevels=10)

    ################################################################################
    ## Conventional model
    ################################################################################
    ## Final simulation test
    # log(mu) = 1 + x1 + x2
    # 
    # What if the skewness about X? (in multivariate)
    # the multivariate skew-normal (MSN) 
    rm(list=ls())
    library(ipeglim)

    # source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/simulateYX.R")
    # btb <- expand.grid(N=c(100), b0=seq(-1.0, 1.0, len=3), b1=1, b2=-1, x1sd=seq(0.5, 1.5, len=3), x2sd=seq(0.5, 1.5, len=3), r=seq(0.0, 0.8, len=3))
    # show a robustnesss over correlations
    btb <- expand.grid(N=c(100), b0=seq(-1.0, 1.0, len=3), b1=1, b2=-1, x1sd=c(1), x2sd=c(1), r=seq(0.0, 0.8, len=3))
    op <- list()
    niters <- 1
    i <- 1

    for(j in 1:nrow(btb)){
      while(i <= niters){
      
        Xmean <- c(1, 0, 0)
        Xsd <- c(0, btb$x1sd[j], btb$x2sd[j])
        XR <- matrix(c(1,0,0, 0,1,btb$r[j], 0,btb$r[j],1), ncol=length(Xmean))
        N0 <- btb$N[j]
        bcoef <- as.vector(c(btb$b0[j], btb$b1[j], btb$b2[j]))
        Xvcv <- diag(Xsd) %*% XR %*% t(diag(Xsd))
        mu <- exp((Xvcv+diag(Xmean))%*%bcoef)[1] 
        miss <- dpois(0, lambda=mu)
        
        tmp <- simulateYX(N=N0, Xreg=TRUE, param=c(btb$b0[j], btb$b1[j], btb$b2[j]), Xstr=list(type="mvnorm", mean=Xmean, sd=Xsd, R=XR), ztrunc=TRUE)
        
        dat <- tmp$yX
        yt <- c(tmp$y)    
        
        fit <- ztpr(formula=y~x1+x2, data=dat, dist="poisson", ztrunc=TRUE, method="BFGS")
        fit <- summary(fit)
        
        sds <- sqrt(diag(fit$vcv))
        cover <- all(c(fit$cil < N0, N0 < fit$ciu))
        op[[niters*(j-1)+i]] <- c(N0=N0, b0=btb$b0[j], b1=btb$b1[j], b2=btb$b2[j], x1sd=btb$x1sd[j], x2sd=btb$x2sd[j], r=btb$r[j], n=nrow(dat), nr=nrow(dat)/N0, scm1=mean(yt), scm2=var(yt), sg1=skewness(yt), scv=cv(yt), N=fit$N, Nr=(fit$N)/N0, cil=fit$cil, cilr=fit$cil/N0, ciu=fit$ciu, ciur=fit$ciu/N0, cover=cover, cf0=fit$cfs[1], cf1=fit$cfs[2], cf2=fit$cfs[3], cf0se=sds[1], cf1se=sds[2], cf2se=sds[3])
        i <- i + 1
      }
      i <- 1
      print(j)
    }
    out <- as.data.frame(do.call(rbind, op))
    out

    ss <- subset(out, subset=(b0==0 & x1sd==1 & x2sd==1 & r==0))
}