Bayesian P-spline
Belay Birlie
I am working on Generalised Additive Mixed Model (GAMM). I want to implement the B-spline smoothness penalty of Eilers and Marx (1996) using a random walk penalty prior (Lang and Brezger, 2004) for longitudinal count data. The count response is measured for each subject at 20 time points. The time point is the term I want to smooth. Since the subjects are measured at the same time point, I think the basis vectors should stay the same across subjects. Is there a way to implement this in INLA?
Below is my trial. I am afraid I didn’t implement the model well. I even assume an identity penalty matrix because z model didn’t work when I use the random walk penalty.
nyy<-20
n.sub<-100
ti<-seq(from=0,to=1,length=nyy)
data_HPN <-matrix(0,n.sub,nyy)
sig_cluster<-0.5
tru_fun<-(1/3)*dbeta(ti, 20, 5)+(1/3)*dbeta(ti, 12, 12)+(1/3)*dbeta(ti, 7, 30)
a.i1<-rnorm(n.sub,0,sig_cluster)
for(i in 1:n.sub)
{
data_HPN[i,]<-rpois(nyy,lambda =exp(tru_fun+a.i1[i]))
}
y1<-c(data_HPN[1,]) # HPN
for (i in 2:n.sub){y1=c(y1,data_HPN[i,])}
time<-c(rep(ti,n.sub))
id<-rep(1:n.sub,each=nyy)
dayno.Z <- bspline(time,K=10,bdeg=3,cyclic=1)-mean(bspline(time,K=10,bdeg=3,cyclic=1)[,2])# the #function defined below
id.z<-1:nrow(dayno.Z)
formula_HPN <- y~f(id.z,model="z",Z=ZSpline,constr=T,diagonal=1e-3,initial=3,param=c(1, 0.0005))+
f(id,model="iid", prior="loggamma", param=c(0.001,0.001))
INLA_HPN = inla(formula_HPN,data=list(y=y1, id.z=id.z,
id=id,ZSpline=dayno.Z),
control.predictor=list(compute=TRUE,link=1),family="poisson",
control.compute=list(dic=TRUE, waic=TRUE, cpo=TRUE))
#B-spline basis function
# -------------------------- Function to construct the Z matrix
bspline <- function(x,K,bdeg,cyclic,xl=min(x),xr=max(x)){
# x - a covariate vector
# K - the number of knots
# bdeg - the degree of the polynomials in the spline
# cyclic - whether to make the basis periodic
# xl, xr - minimum and maximum values in the covariate space
x <- as.matrix(x,ncol=1) # reshape to be a column
ndx <- K - bdeg # how many knots will be left at the end?
# as outlined in Eilers and Marx (1996)
dx <- (xr - xl) / ndx # step size
t <- xl + dx * (-bdeg:(ndx+bdeg)) # place knots
# form spline basis of order 0
T <- (0 * x + 1) %*% t
X <- x %*% (0 * t + 1)
P = (X - T) / dx
B = (T <= X) & (X < (T + dx))
r = c(2:length(t), 1) # reordering of Indices
# recursive updating of basis, increasing degree each step
for (k in 1:bdeg){
B = (P * B + (k + 1 - P) * B[ ,r]) / k;
}
# only return the first K columns
B <- B[,1:(ndx+bdeg)]
# convert to periodic basis
if (cyclic == 1){
for (i in 1:bdeg){
B[ ,i] = B[ ,i] + B[ ,K-bdeg+i]
}
# get rid of the bits that are being reused
B <- B[ , 1:(K-bdeg)]
}
return(B)
}
Helpdesk
we prefer using the random walk models to do the same, like 'rw2'. You
can read about it in the GMRF book. There is no native support to do p-
splines, so in that case you'll need to use one of the generic models.
computationally, I think you'll be better off with the rw2 models, but
it should not matter much.
Best
H
> You received this message because you are subscribed to the Google
> Groups "R-inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it,
> send an email to r-inla-discussion...@googlegroups.com
> .
> To post to this group, send email to r-inla-discussion-group@googlegr
> oups.com.
> Visit this group at https://groups.google.com/group/r-inla-discussion
> -group.
> For more options, visit https://groups.google.com/d/optout.
--
Håvard Rue
Helpdesk
he...@r-inla.org