rapi.saemix: Rutgers Alcohol Problem Index

In saemix: Stochastic Approximation Expectation Maximization (SAEM) Algorithm

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)
}
countsimulate.poisson<-function(psi, id, xidep) {
   time<-xidep[,1]
   y<-xidep[,2]
   ymax<-max(y)
   intercept<-psi[id,1]
   slope<-psi[id,2]
   lambda<- exp(intercept + slope*time)
   y<-rpois(length(time), lambda=lambda)
   y[y>ymax]<-ymax+1 # truncate to maximum observed value to avoid simulating aberrant values
   return(y)
}
# Gender effect on intercept and slope
rapimod.poisson<-saemixModel(model=count.poisson, simulate.function=countsimulate.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)
}
countsimulate.poissonzip<-function(psi, id, xidep) {
  time<-xidep[,1]
  y<-xidep[,2]
  ymax<-max(y)
  intercept<-psi[id,1]
  slope<-psi[id,2]
  p0<-psi[id,3] # Probability of zero's
  lambda<- exp(intercept + slope*time)
  prob0<-rbinom(length(time), size=1, prob=p0)
  y<-rpois(length(time), lambda=lambda)
  y[prob0==1]<-0
  y[y>ymax]<-ymax+1 # truncate to maximum observed value to avoid simulating aberrant values
  return(y)
}
rapimod.zip<-saemixModel(model=count.poissonzip, simulate.function=countsimulate.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)

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

saemix documentation built on May 29, 2024, 10:43 a.m.