demo/xreg-eq3d.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## Overall framework is over;
## ONLY ADD COMMENTARY REMAINS - 2013.10.22
##
## ipeglim/demo/poisxreg3p.R
## 
## y ~ b1x1 + b2x2 + b3x3 (no intercept model)

library(ipeglim)
rm(list=ls())

# structure of X
xm <- c(0,0,0)
xsd <- c(1,1,1)
xr <- diag(3)
xr[1,2] <- xr[2,1] <- r1 <- 0
xr[1,3] <- xr[3,1] <- r2 <- r1
xr[2,3] <- xr[3,2] <- r3 <- r1
xsigma <- t(diag(xsd))%*%xr%*%diag(xsd)

# structure of B
bm <- c(0,0,0)
bsd <- c(1,1,1)
br <- diag(3)
br[1,2] <- br[2,1] <- r1 <- 0
br[1,3] <- br[3,1] <- r2 <- r1
br[2,3] <- br[3,2] <- r3 <- r1
Bmat <- t(diag(bsd))%*%br%*%diag(bsd)
Bmat <- Bmat*1e-2

# no intercept model
tmp <- rvg4yx(N=1e3, pars=c(1.0, -0.5, 0.5), ztrunc=TRUE, xreg=TRUE, kind="mvnorm", mean=xm, sigma=xsigma, center.X=FALSE)
dat <- tmp$yX
print(table(dat$y))

m2fit <- model(formula=y~x1+x2+x3-1, data=dat, ztrunc=TRUE, dist="poisson")
lc0 <- list(lhs=rbind(diag(3), -diag(3)), rhs=c(rep(-1,6)))
cmfit <- iprior(obj=m2fit, eqns=lc0)
plot(cmfit, # for m0xtms
  main="Cube Type of Imprecise Prior \nspecified by a set of inequalities", 
  sub=list(expression(-1<beta[1]~","~beta[2]~","~beta[3]<1), cex=1.25),
  xlab=expression(beta[1]),
  ylab=expression(beta[2]),
  zlab=expression(beta[3]), col="darkblue")
xtms0 <- cmfit$xtms

cmfit <- setGrid(obj=cmfit, len=7)
plot(cmfit, # for m0xtms
  main="Cube Type of Imprecise Prior Specification\nwith Grid Points", 
  sub=list(expression(-1<beta[1]~","~beta[2]~","~beta[3]<1), cex=1.25),
  xlab=expression(beta[1]),
  ylab=expression(beta[2]),
  zlab=expression(beta[3]), col="darkblue")
xtms1 <- cmfit$xtms

sop <- summary(update(obj=cmfit, method="LAF", B=Bmat, verbose=TRUE))
plot(sop, main="Imprecise Posterior Shrinkage")

# contour plot + wireframe (Is this neccessary?)
# problematic with sphere3d
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/cwplot.R")
cwplot(sop, xlab="", ylab="") 

xid <- which(xtms1 %in% xtms0) # convex hull
BETAs <- pbox(x=sop, which.beta=1, which.xtms=xid, smooth=TRUE)
BB <- extractMCHAIN(BETAs)
pbox(x=BB, which.beta=2, which.xtms=xid, smooth=TRUE)
pbox(x=BB, which.beta=3, which.xtms=xid, smooth=TRUE)


##################################3

## PART 2


## Overall framework is over;
## ONLY ADD COMMENTARY REMAINS - 2013.10.22
##
## ipeglim/demo/poisxreg3p.R
## 
## y ~ b1x1 + b2x2 + b3x3 (no intercept model)

library(ipeglim)
rm(list=ls())

# structure of X
xm <- c(0,0,0)
xsd <- c(1,1,1)
xr <- diag(3)
xr[1,2] <- xr[2,1] <- r1 <- 0
xr[1,3] <- xr[3,1] <- r2 <- r1
xr[2,3] <- xr[3,2] <- r3 <- r1
xsigma <- t(diag(xsd))%*%xr%*%diag(xsd)

# structure of B
bm <- c(0,0,0)
bsd <- c(1,1,1)
br <- diag(3)
br[1,2] <- br[2,1] <- r1 <- 0
br[1,3] <- br[3,1] <- r2 <- r1
br[2,3] <- br[3,2] <- r3 <- r1
Bmat <- t(diag(bsd))%*%br%*%diag(bsd)
Bmat <- Bmat*1e-2

# no intercept model
tmp <- rvg4yx(N=1e3, pars=c(1.0, -0.5, 0.5), ztrunc=TRUE, xreg=TRUE, kind="mvnorm", mean=xm, sigma=xsigma, center.X=FALSE)
dat <- tmp$yX
print(table(dat$y))

m2fit <- model(formula=y~x1+x2+x3-1, data=dat, ztrunc=TRUE, dist="poisson")
lc0 <- list(lhs=rbind(diag(3), -diag(3)), rhs=c(rep(-1,6)))
cmfit <- iprior(obj=m2fit, eqns=lc0)
# cmfit <- setGrid(obj=cmfit, len=7)

sop <- summary(update(obj=cmfit, method="LAF", B=Bmat, verbose=TRUE))
plot(sop, main="Imprecise Posterior Shrinkage")

if(FALSE){
cwplot(sop, xlab="", ylab="") 
BETAs <- pbox(x=sop, which.beta=1, which.xtms=xid, smooth=TRUE)
BB <- extractMCHAIN(BETAs)
pbox(x=BB, which.beta=3, which.xtms=xid, smooth=TRUE)
}