if (FALSE) {
##
library(JMbayes2)
library(lme4)
library(insight)
library(numDeriv)
library(stats)
lmer(log(alkaline)~drug+age+year+(1|id),data=na.omit(pbc2))
data1<-pbc2
data1$alkaline<-log(data1$alkaline)
names(pbc2)
apply(pbc2,2,function(x){sum(is.na(x))})
r.ij<-ifelse(is.na(data1$alkaline)==T,0,1)
data1<-cbind.data.frame(data1,r.ij)
data1$drug<-factor(data1$drug,levels=c("placebo","D-penicil"),labels = c(0,1))
data1$sex<-factor(data1$sex,levels=c('male','female'),labels=c(1,0))
data1$drug<-as.numeric(as.character(data1$drug))
data1$sex<-as.numeric(as.character(data1$sex))
r.ij~year+age+sex+drug+serBilir+(1|id)
model.r<-glmer(r.ij~year+age+sex+drug+serBilir+(1|id),family=binomial(link='logit'),data=data1)
model.y<-lmer(alkaline~year+age+sex+drug+serBilir+(1|id),data=na.omit(data1))
nu<-model.r@beta
psi<-model.y@beta
sigma<-get_variance_residual(model.y)
sigmaR<-get_variance(model.r)$var.random
sigmaMiss<-get_variance(model.y)$var.random
m11<-lmeaipw(data=data1,id='id',
analysis.model = alkaline~year,
wgt.model=~year+age+sex+drug+serBilir+(1|id),
imp.model = ~year+age+sex+drug+serBilir+(1|id),
psiCov = vcov(model.y),nu=nu,psi=psi,
sigma=sigma,sigmaMiss=sigmaMiss,sigmaR=sigmaR,dist='gaussian',link='identity',
maxiter = 200)
m11
##
}
Run the code above in your browser using DataLab