growth.R

################################################################################
#    This file is part of growthrate.                                          #
#                                                                              #
#    growthrate is free software: you can redistribute it and/or modify        #
#    it under the terms of the GNU General Public License as published by      #
#    the Free Software Foundation, either version 3 of the License, or         #
#    (at your option) any later version.                                       #
#                                                                              #
#    growthrate is distributed in the hope that it will be useful,             #
#    but WITHOUT ANY WARRANTY; without even the implied warranty of            #
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             #
#    GNU General Public License for more details.                              #
#                                                                              #
#    You should have received a copy of the GNU General Public License         #
#    along with growthrate.  If not, see <http://www.gnu.org/licenses/>.       #
################################################################################

growth <- function(data,tobs,sigma,d) {
data = as.matrix(data);
obsdata = data;
dim(obsdata);

# obsdata is the data matrix of dimension N (number of subjects) times n (number of time points).
N = dim(obsdata)[1];
n = dim(obsdata)[2];

# calculate the difference of consecutive time points
dobspts = diff(tobs);

# calculate the vector y in the paper
YI = matrix(nrow=N,ncol=(n-1));
dobsdatak = rep(0,(n-1));

for (k in 1:N) {
YI[k,] = diff(obsdata[k,])/dobspts;
}

# Estimate the prior mean based on nonparametric approach and obtain Xtilda (ybar in the paper)
resulmuprior = priormeanobs(YI,tobs);
muprior = resulmuprior$muprior;
Xtilda = resulmuprior$Xtilda;

# Estimate the prior precision based on CLIME
re.clime = clime(t(Xtilda),perturb=TRUE);
re.cv = cv.clime(re.clime,loss=c("likelihood", "tracel2"));
re.clime.opt <- clime(t(Xtilda), re.cv$lambdaopt, perturb =TRUE);
Sigma0inv <- re.clime.opt$Omegalist[[1]];

# Get the posterior mean and variance-covariance matrix at observed time points
resulpost = posteriorobs(Sigma0inv,sigma,muprior,Xtilda,tobs,YI);
muhatMatrix = resulpost$muhatMatrix;
Sigmahat = resulpost$Sigmahat;

# The posterior mean curves (muhatcurve) and variance covariance operator (Khat) evaluated at the fine grid
r = posteriordistribcurve(muhatMatrix,Sigmahat,sigma,tobs,d,YI);
tgrid = r$tgrid;
muhatcurve = r$muhatcurvematrix;
Khat = Khatf(tgrid,tobs,sigma,Sigmahat);

# Return the muhat curves for each subject, variance-covariance operator and the fine grid.
result = list(muhatcurve=muhatcurve,Khat=Khat,tgrid=tgrid);
return(result);
}

Project Homepage: