# R/growth.R In growthrate: Bayesian reconstruction of growth velocity

#### Documented in growth

```################################################################################
#    This file is part of growthrate.                                          #
#                                                                              #
#    growthrate is free software: you can redistribute it and/or modify        #
#    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);
n = dim(obsdata);

## 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[];

## 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);
}
```

## Try the growthrate package in your browser

Any scripts or data that you put into this service are public.

growthrate documentation built on May 2, 2019, 3:43 p.m.