Nothing
## 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))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.