merged

Author: karchjd

.Rbuildignore

@@ -1,2 +1,7 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
+^data-raw$
+.gitignore
+^\.travis\.yml$
+README.Rmd
+README.md

.gitignore

@@ -2,6 +2,9 @@
 .Rhistory
 .RData
 .Ruserdata
-inst/doc
 README.html
 .DS_Store
+*.pdf
+*.tex
+vignettes/cache
+vignettes/gppm.bbl

.travis.yml

@@ -0,0 +1,5 @@
+# R for travis: see documentation at https://docs.travis-ci.com/user/languages/r
+
+language: R
+sudo: false
+cache: packages

DESCRIPTION

@@ -1,20 +1,28 @@
-Package: gppmr
-Title: Preliminary Implementation of GPPM
-Version: 0.1.0
-Authors@R: c(person(c("Julian", "D."), "Karch", email="[email protected]", role=c("aut", "cre", "cph")))
-Description: Preliminary Implementation of Gaussian Process Panel Modeling.
+Package: gppm
+Version: 0.2.0.9000
+Title: Gaussian Process Panel Modeling
+Description:  Provides an implementation of Gaussian process panel modeling (GPPM), which enables: Gaussian process based modeling of longitudinal panel data (GPPM), regular Gaussian process regression with flexible model specification, and Gaussian process regression with multi-task learning.
+Authors@R: c(person(c("Julian", "D."), "Karch", email="[email protected]", role=c("aut", "cre", "cph")))
 Depends:
-    R (>= 3.3.1),
-    OpenMx (>= 2.6.9),
+    R (>= 3.1.0),
+    Rcpp (>= 0.12.17)
 Imports:
-    MASS
-License: GPL-3
-Encoding: UTF-8
+    rstan (>= 2.17.3),
+    ggplot2 (>= 2.2.1),
+    MASS (>= 7.3-49),
+    ggthemes (>= 3.5.0),
+    mvtnorm (>= 1.0-8),
+    stats,
+    methods
+Suggests: 
+    testthat (>= 2.0.0),
+    knitr (>= 1.20),
+    rmarkdown (>= 1.10),
+    roxygen2 (>= 6.0.1)  
+License: GPL-3 | file LICENSE
 LazyData: true
-URL: https://github.com/karchjd/gppmr
-BugReports: https://github.com/karchjd/gppmr/issues
+URL: https://github.com/karchjd/gppm
+BugReports: https://github.com/karchjd/gppm/issues
 RoxygenNote: 6.0.1
-Suggests: testthat (>= 1.0.2),
-    knitr,
-    rmarkdown,
-VignetteBuilder: knitr
+Roxygen: list(markdown = TRUE)
+Encoding: UTF-8

NAMESPACE

@@ -1,8 +1,42 @@
 # Generated by roxygen2: do not edit by hand
 
-export(generateAR)
-export(generateLGCM)
-export(gppFit)
-export(gppModel)
-export(simulateData)
-import(OpenMx)
+S3method(coef,GPPM)
+S3method(confint,GPPM)
+S3method(fit,GPPM)
+S3method(fitted,GPPM)
+S3method(logLik,GPPM)
+S3method(plot,GPPMPred)
+S3method(plot,LongData)
+S3method(predict,GPPM)
+S3method(print,summary.GPPM)
+S3method(simulate,GPPM)
+S3method(summary,GPPM)
+S3method(vcov,GPPM)
+export(SE)
+export(accuracy)
+export(covf)
+export(createLeavePersonsOutFolds)
+export(crossvalidate)
+export(datas)
+export(fit)
+export(getIntern)
+export(gppm)
+export(gppmControl)
+export(maxnObs)
+export(meanf)
+export(nObs)
+export(nPars)
+export(nPers)
+export(nPreds)
+export(parEsts)
+export(pars)
+export(preds)
+import(Rcpp)
+import(ggplot2)
+import(ggthemes)
+import(rstan)
+import(stats)
+importFrom(MASS,mvrnorm)
+importFrom(methods,is)
+importFrom(mvtnorm,dmvnorm)
+importFrom(utils,capture.output)

R/FitRes.R

@@ -0,0 +1,57 @@
+new_FitRes <- function(paraEsts,vcov,LL,nPar,mu,Sigma){
+  stopifnot(is.double(paraEsts)) #maximum likelihood estimates
+  stopifnot(is.matrix(vcov)) #covariance of sampling distribution of maximum likelihood estimator
+  stopifnot(is.double(LL) && length(LL)==1)
+  stopifnot(is.integer(nPar) && length(nPar)==1)
+  stopifnot(is.list(mu))
+  stopifnot(is.list(Sigma))
+
+  structure(list(
+    paraEsts=paraEsts,
+    vcov=vcov,
+    LL=LL,
+    nPar=nPar,
+    mu=mu,
+    Sigma=Sigma
+  ),
+  class='StanData'
+  )
+}
+
+extractMoments <- function(stanOutParas,dataStats){
+  mu <- rep(list(numeric(dataStats$maxTime)),dataStats$nPer)
+  Sigma <- rep(list(matrix(nrow=dataStats$maxTime,ncol = dataStats$maxTime)),dataStats$nPer)
+  for (iPer in seq_len(dataStats$nPer)){
+    mu[[iPer]] <- stanOutParas$mu[iPer,1:dataStats$nTime[iPer]]
+    Sigma[[iPer]] <- stanOutParas$Sigma[iPer,1:dataStats$nTime[iPer],1:dataStats$nTime[iPer]]
+  }
+  res <- list(mu=mu,Sigma=Sigma,IDs=attr(dataStats,'IDs'))
+}
+
+
+extractLL <- function(meanCov,dataStats){
+  Y <- dataStats$Y
+  ll <- 0
+  for (i in 1:length(Y)){
+    ll <- ll + mvtnorm::dmvnorm(Y[[i]][1:dataStats$nTime[i]],mean=meanCov$mu[[i]],sigma=as.matrix(meanCov$Sigma[[i]]),log=TRUE)
+  }
+  ll
+}
+
+extractFitRes <- function(stanOut,parsedModel,dataStats){
+  paraEsts <- stanOut$par[parsedModel$params]
+  theNames <- names(paraEsts)
+  paraEsts <- as.numeric(paraEsts)
+  names(paraEsts) <- theNames
+  vcov <- as.matrix(NA)
+  if (!is.null(stanOut$hessian)){
+    tryCatch(vcov <- solve(-stanOut$hessian),error=function(e){})
+  }
+
+  nPar <- length(parsedModel$params)
+  meanCov <- extractMoments(stanOut$par,dataStats)
+  LL <- extractLL(meanCov,dataStats)
+  mu=meanCov$mu
+  Sigma=meanCov$Sigma
+  new_FitRes(paraEsts,vcov,LL,nPar,mu,Sigma)
+}

R/GPPM.R

@@ -0,0 +1,155 @@
+new_GPPM <- function(mFormula,cFormula,myData,control){
+  stopifnot(is.character(mFormula))
+  stopifnot(is.character(cFormula))
+  stopifnot(is.data.frame(myData))
+  stopifnot(class(control)=='GPPMControl')
+
+  structure(list(
+    mFormula=mFormula, #formula for the mean
+    cFormula=cFormula, #formula for the covariance
+    data=myData,       #data must be a data frame
+    control=control,    #list of controls
+    parsedModel=NA,    #model in a parsed format
+    dataForStan=NA,    #data as used for stan
+    stanModel=NA,      #generated stan Model
+    stanOut =NA,       #stan output
+    fitRes=NA          #all the fitting results
+  ),class='GPPM')
+}
+
+
+#' Define a Gaussian process panel model
+#'
+#' This function is used to specify a Gaussian process panel model (GPPM),
+#' which can then be fit using \code{\link{fit.GPPM}}.
+#'
+#' @param mFormula character string. Contains the specification of the mean function. See details for more information.
+#'
+#' @param cFormula character string. Contains the specification of the covariance function. See details for more information.
+#'
+#' @param myData data frame. Contains the data to which the model is fitted. Must be in the long-format.
+#'
+#' @param ID character string. Contains the column label in myData which describes the subject ID.
+#'
+#' @param DV character string. Contains the column label in myData which contains the to be modeled variable.
+#'
+#' @param control object of class GPPMControl. Used for storing technical settings. Default should only be changed by advanced users. Generated via \code{\link{gppmControl}}.
+#'
+#' @return A (unfitted) Gaussian process panel model, which is an object of class 'GPPM'
+#' @details
+#' mFormula and cFormula contain the specification of the mean and the covariance function respectively.
+#' These formulas are defined using character strings. Within these strings there are four basic elements:
+#'  \itemize{
+#'   \item Parameters
+#'   \item Functions and operators
+#'   \item References to observed variables in the data frame myData
+#'    \item Mathematical constants
+#' }
+#' The gppm function automatically recognizes which part of the string refers to which elements. To be able to do this certain relatively common rules need to be followed:
+#'
+#' Parameters: Parameters may not have the same name as any of the columns in myData to avoid confusing them with a reference to an observed variable.
+#'  Furthermore, to avoid confusing them with functions, operators, or constants, parameter labels must always begin with a lower case letter and only contain letters and digits.
+#'
+#' Functions and operators: All functions and operators that are supported by stan can be used; see \url{http://mc-stan.org/users/documentation/} for a full list. In general, all basic operators and functions are supported.
+#'
+#' References: A reference must be the same as one of the elements of the output of \code{names(myData)}. For references, the same rules apply as for parameters. That is, the column names of myData may only contain letters and digits and must start with a letter.
+#'
+#' Constants: Again, all constants that are supported by stan can be used and in general the constants are available by their usual name.
+#' @seealso \code{\link{fit.GPPM}} for how to fit a GPPM
+#' @examples
+#' # Defintion of a latent growth curve model
+#' \dontrun{
+#' data("demoLGCM")
+#' lgcm <- gppm('muI+muS*t','varI+covIS*(t+t#)+varS*t*t#+(t==t#)*sigma',
+#'         demoLGCM,'ID','y')
+#' }
+#' @import rstan
+#' @import ggplot2
+#' @import ggthemes
+#' @import Rcpp
+#' @importFrom mvtnorm dmvnorm
+#' @import stats
+#' @importFrom MASS mvrnorm
+#' @importFrom methods is
+#' @importFrom utils capture.output
+#' @export
+gppm <- function(mFormula,cFormula,myData,ID,DV,control=gppmControl()){
+  myData <- as_LongData(myData,ID,DV)
+  validate_gppm(mFormula,cFormula,myData,control)
+
+  theModel <- new_GPPM(mFormula,cFormula,myData,control)
+  theModel$dataForStan <- as_StanData(myData)
+  theModel$parsedModel <- parseModel(theModel$mFormula,theModel$cFormula,theModel$dataForStan)
+  theModel$stanModel <- toStan(theModel$parsedModel,control)
+  return(theModel)
+}
+
+subsetData <- function(gpModel,rowIdxs){
+  newModel <- gpModel
+  newModel$data <- gpModel$data[rowIdxs,]
+  newModel$dataForStan <- as_StanData(newModel$data)
+  return(newModel)
+}
+
+updateData <- function(gpModel,newData){
+  newModel <- gpModel
+  oldData <- datas(newModel)
+  stopifnot(identical(names(oldData),names(newData)))
+  newModel$data <- as_LongData(newData,getID(oldData),getDV(oldData))
+  newModel$dataForStan <- as_StanData(newModel$data)
+  return(newModel)
+}
+
+validate_gppm <- function(mFormula,cFormula,myData,control){
+  ID <- attr(myData,'ID')
+  DV <- attr(myData,'DV')
+  stopifnot(!is.null(ID))
+  stopifnot(!is.null(DV))
+  #type checks
+  if (!is.character(mFormula)){
+    stop('mFormula must contain a string')
+  } else if(length(mFormula)!=1){
+
+  }
+
+  if (!(is.character(cFormula)&& length(cFormula)==1)){
+    stop('cFormula must contain a string')
+  }
+
+  if (!is.data.frame(myData)){
+    stop('myData must be a data frame')
+  }
+
+  if (!(is.character(ID))){
+    browser()
+    stop('ID must contain a string')
+  }
+
+  if (!is.character(DV)){
+    stop('DV must contain a string')
+  }
+
+  if (!class(control)=='GPPMControl'){
+    stop('control must be of class GPPMControl')
+  }
+
+  varNames <- names(myData)
+  allValid <- grepl("^[A-Za-z]+[0-9A-Za-z]*$",varNames)
+  if (any(!allValid)){
+    stop(sprintf('Invalid variable name %s in your data frame. See ?gppm for naming conventions\n',varNames[!allValid]))
+  }
+
+
+  if (!"longData" %in% class(myData)){
+    if (!ID %in% names(myData)){
+      stop(sprintf('ID variable %s not in data frame',ID))
+    }
+
+    if (!DV %in% names(myData)){
+      stop(sprintf('DV variable %s not in data frame',DV))
+    }
+  }
+}
+
+
+

R/GPPMControl.R

@@ -0,0 +1,22 @@
+new_GPPMControl <- function(stanModel){
+  stopifnot(is.logical(stanModel))
+
+
+  structure(list(
+    stanModel=stanModel #generate stanModel Object or not
+  ),class='GPPMControl')
+}
+
+#' Define settings for a Gaussian process panel model
+#'
+#' This function is used to specify the settings of a Gaussian process panel model generated by \code{\link{gppm}}.
+#'
+#' @param stanModel boolean. Should the corresponding stan model be created? Should only be set to FALSE for testing purposes.
+#' Not creating the stan model makes model fitting impossible but saves a lot of time.
+
+#' @return Settings for a Gaussian process panel model in an object of class 'GPPMControl'
+#' @seealso \code{\link{gppm}}
+#' @export
+gppmControl <- function(stanModel=TRUE){
+  new_GPPMControl(stanModel)
+}

R/LongData.R

@@ -0,0 +1,30 @@
+new_LongData <- function(myData,ID,DV){
+  stopifnot(is.character(ID) && length(ID)==1)
+  stopifnot(is.character(DV) && length(DV)==1)
+  stopifnot(is.data.frame(myData))
+
+  structure(myData,class=c('LongData',class(myData)),ID=ID,DV=DV)
+}
+
+as_LongData <- function(myData,ID,DV){
+  if (!"LongData" %in% class(myData)){
+    myData <- new_LongData(myData,ID,DV)
+  }else if (!missing(ID)){
+    attr(myData,'ID') <- ID
+  }else if (!missing(DV)){
+    attr(myData,'DV') <- DV
+  }
+  myData
+}
+
+
+getID <- function(longData){
+  stopifnot("LongData" %in% class(longData))
+  return(attr(longData,'ID'))
+}
+
+getDV <- function(longData){
+  stopifnot("LongData" %in% class(longData))
+  return(attr(longData,'DV'))
+}
+