demo/xre2i.R

## GOOD TO GO TO RELEASE 2013.10.22
##
## ipeglim/demo/xreg-eq2d.R
##
## Poisson regression model with a multivariate normal imprecise prior
## * X is centered;
## * log(E(y)) ~ b0 + b1*x1
## * log(E(y)) ~ b1*x1 + b2*x2
## * log(E(y)) ~ b0 + b1*x1 + b2*x2
## * log(E(y)) ~ b1*x1 + b2*x2 + b3*x3
## tested by methods = 'LA', 'MH', 'IS'.
## 

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

## case 1. log(E(y)) ~ b0 + b1*x1
## default structure of X
xm <- c(0)
xsd <- c(1)

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

# parameter specification
bcfs1 <- c(0.5, 0.5)      # 1 quadrant
bcfs2 <- c(-0.5, 0.5)     # 2 quadrant
bcfs3 <- c(-0.5, -0.5)    # 3 quadrant
bcfs4 <- c(0.5, -0.5)     # 4 quadrant

# random variates generation
tmp <- rvg4yx(N=5e2, cf0=bcfs1[1], pars=bcfs1[2], ztrunc=TRUE, xreg=TRUE, kind="mvnorm", mean=xm, sd=xsd, center.X=TRUE)
dat <- tmp$yX
table(dat$y)
xi1 <- c(1,unlist(dat[1,-1]))

m2fit <- model(formula=y~x1, data=dat, ztrunc=TRUE, dist="poisson")
lc0 <- list(lhs=rbind(diag(2), -diag(2)), rhs=c(rep(-1,4)))
cmfit <- iprior(obj=m2fit, eqns=lc0)

ola <- summary(update(obj=cmfit, method="LA", B=Bmat, xi=xi1, verbose=TRUE))

X11()
plot(ola, # for m0xtms
  main="Shrinkage - Box-Constrained", 
  sub=list(expression(-1<beta[0]~","~beta[1]<1), cex=1.25),
  xlab=expression(beta[1]),
  ylab=expression(beta[2]))

## case 1. log(E(y)) ~ b1*x1 + b2*x2
# structure of X
xm <- c(3,-2)
xsd <- c(1,1)
xr <- diag(2)
xr[2] <- xr[3] <- r <- 0
xsigma <- t(diag(xsd))%*%xr%*%diag(xsd)

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

# random variates generation
tmp <- rvg4yx(N=5e2, pars=c(0.5, -0.5), ztrunc=TRUE, xreg=TRUE, kind="mvnorm", mean=xm, sigma=xsigma, center.X=TRUE)
dat <- tmp$yX
table(dat$y)

m2fit <- model(formula=y~x1+x2-1, data=dat, ztrunc=TRUE, dist="poisson")
lc0 <- list(lhs=rbind(diag(2), -diag(2)), rhs=c(rep(-1,4)))
cmfit <- iprior(obj=m2fit, eqns=lc0)

# plotting extreme points
plot(cmfit, # for m0xtms
  main="Rectangle Type of Imprecise Priors\nconstrained by a set of inequalities", 
  sub=list(expression(-1<beta[1]~","~beta[2]<1), cex=1.25),
  xlab=expression(beta[1]),
  ylab=expression(beta[2]))
xtms0 <- cmfit$xtms

# plotting extreme points on the xy-grid 
cmfit <- setGrid(obj=cmfit, len=10)
plot(cmfit, # for m0xtms
  main="Rectangle Type of Imprecise Priors\nconstrained by a set of inequalities (with grid points)", 
  sub=list(expression(-1<beta[1]~","~beta[2]<1), cex=1.25),
  xlab=expression(beta[1]),
  ylab=expression(beta[2]))

# imprecise model fitting 
xi1 <- unlist(dat[1,-1])
op <- update(obj=cmfit, method="LA", B=Bmat, xi=xi1, verbose=TRUE)
sop <- summary(op, HT.est=TRUE)
sop
sop$sups
sop$infs
sop$infs.idx
sop$sups.idx

plot(sop, # for m0xtms
  main="Imprecise Posterior Shrinkage", 
  sub=list(expression(-1<beta[1]~","~beta[2]<1), cex=1.25),
  xlab=expression(beta[1]),
  ylab=expression(beta[2]))


BETAs <- pbox(x=sop, which.beta=2, which.xtms=c(sop$infs.idx, sop$sups.idx))
BB <- extractMCHAIN(BETAs)
pbox(x=BB, which.beta=1, which.xtms=c(sop$infs.idx, sop$sups.idx))

Try the ipeglim package in your browser

Any scripts or data that you put into this service are public.

ipeglim documentation built on May 2, 2019, 4:31 p.m.