rapi.saemix: Rutgers Alcohol Problem Index

rapi.saemixR Documentation

Rutgers Alcohol Problem Index

Description

The RAPI data studies gender differences across two years in alcohol-related problems, as measured by the Rutgers Alcohol Problem Index (RAPI; White & Labouvie, 1989). The dataset includes 3,616 repeated measures of counts representing the number of alcohol problems reported over six months period, across five time points from 818 individuals.

Format

This data frame contains the following columns:

id

subject identification number

time

time since the beginning of the study (months)

rapi

the number of reported alcohol problems during a six-months period

gender

gender (0=Men, 1=Women

Source

David Atkins, University of Washington

References

D Atkins, S Baldwin, C Zheng, R Gallop, C Neighbors C (2013). A tutorial on count regression and zero-altered count models for longitudinal substance use data. Psychology of Addictive Behaviors, 27(1):166–177.

C Neighbors, Lewis MA, Atkins D, Jensen MM, Walter T, Fossos N, Lee C, Larimer M (2010). Efficacy of web-based personalized normative feedback: A two-year randomized controlled trial. Journal of Consulting and Clinical Psychology 78(6):898-911.

C Neighbors, N Barnett, D Rohsenow, S Colby, P Monti (2010). Cost-Effectiveness of a Motivational lntervention for Alcohol-Involved Youth in a Hospital Emergency Department. Journal of Studies on Alcohol and Drugs 71(3):384-394.

Examples

 data(rapi.saemix)
 saemix.data<-saemixData(name.data=rapi.saemix, name.group=c("id"),
                     name.predictors=c("time","rapi"),name.response=c("rapi"),
                     name.covariates=c("gender"),units=list(x="months",y="",covariates=c("")))
hist(rapi.saemix$rapi, main="", xlab="RAPI score", breaks=30)


# Fitting a Poisson model
count.poisson<-function(psi,id,xidep) { 
  time<-xidep[,1]
  y<-xidep[,2]
  intercept<-psi[id,1]
  slope<-psi[id,2]
  lambda<- exp(intercept + slope*time)
  logp <- -lambda + y*log(lambda) - log(factorial(y))
  return(logp)
}
# Gender effect on intercept and slope
rapimod.poisson<-saemixModel(model=count.poisson,
   description="Count model Poisson",modeltype="likelihood",   
   psi0=matrix(c(log(5),0.01),ncol=2,byrow=TRUE,dimnames=list(NULL, c("intercept","slope"))), 
   transform.par=c(0,0), omega.init=diag(c(0.5, 0.5)),
    covariance.model =matrix(data=1, ncol=2, nrow=2),
    covariate.model=matrix(c(1,1), ncol=2, byrow=TRUE))
saemix.options<-list(seed=632545,save=FALSE,save.graphs=FALSE, displayProgress=FALSE, fim=FALSE)
poisson.fit<-saemix(rapimod.poisson,saemix.data,saemix.options)

# Fitting a ZIP model
count.poissonzip<-function(psi,id,xidep) {
  time<-xidep[,1]
  y<-xidep[,2]
  intercept<-psi[id,1]
  slope<-psi[id,2]
  p0<-psi[id,3] # Probability of zero's
  lambda<- exp(intercept + slope*time)
  logp <- log(1-p0) -lambda + y*log(lambda) - log(factorial(y)) # Poisson
  logp0 <- log(p0+(1-p0)*exp(-lambda)) # Zeroes
  logp[y==0]<-logp0[y==0]
  return(logp)
}
rapimod.zip<-saemixModel(model=count.poissonzip,
   description="count model ZIP",modeltype="likelihood",   
   psi0=matrix(c(1.5, 0.01, 0.2),ncol=3,byrow=TRUE,
   dimnames=list(NULL, c("intercept", "slope","p0"))), 
   transform.par=c(0,0,3), covariance.model=diag(c(1,1,0)), omega.init=diag(c(0.5,0.3,0)),
   covariate.model = matrix(c(1,1,0),ncol=3, byrow=TRUE))
zippoisson.fit<-saemix(rapimod.zip,saemix.data,saemix.options)

# Simulation functions to simulate from the models
saemix.simulatePoisson<-function(psi, id, xidep) {
  time<-xidep[,1]
  y<-xidep[,2]
  intercept<-psi[id,1]
  slope<-psi[id,2]
  lambda<- exp(intercept + slope*time)
  y<-rpois(length(time), lambda=lambda)
  return(y)
}
saemix.simulatePoissonZIP<-function(psi, id, xidep) {
  time<-xidep[,1]
  y<-xidep[,2]
  intercept<-psi[id,1]
  slope<-psi[id,2]
  p0<-psi[id,3] 
  lambda<- exp(intercept + slope*time)
  prob0<-rbinom(length(time), size=1, prob=p0)
  y<-rpois(length(time), lambda=lambda)
  y[prob0==1]<-0
  return(y)
}

# Using simulations to compare the predicted proportion of 0's in the two models
nsim<-100
yfit1<-simulateDiscreteSaemix(poisson.fit, saemix.simulatePoisson, nsim=nsim)
yfit2<-simulateDiscreteSaemix(zippoisson.fit, saemix.simulatePoissonZIP, 100)
{
nobssim<-length(yfit1@sim.data@datasim$ysim)
cat("Observed proportion of 0's", 
   length(yfit1@data@data$rapi[yfit1@data@data$rapi==0])/yfit1@data@ntot.obs,"\n")
cat("      Poisson model, p=",
   length(yfit1@sim.data@datasim$ysim[yfit1@sim.data@datasim$ysim==0])/nobssim,"\n")
cat("          ZIP model, p=",
   length(yfit2@sim.data@datasim$ysim[yfit2@sim.data@datasim$ysim==0])/nobssim,"\n")
}
  
  

saemix documentation built on Aug. 5, 2022, 5:25 p.m.

Related to rapi.saemix in saemix...