The lava and mets packages

Table of Contents

Multivariate Event Times (mets)

Implementation of various statistical models for multivariate event history data. Including multivariate cumulative incidence models, and bivariate random effects probit models (Liability models)

Installation

install.packages("mets",dependencies=TRUE)

The development version may be installed directly from R-forge (requires Rtools on windows and development tools (+Xcode) for Mac OS X):

install.packages("mets", repos="http://R-Forge.R-project.org")

Citation

To cite the mets package please use the following reference

Thomas H. Scheike and Klaus K. Holst and Jacob B. Hjelmborg (2013). 
Estimating heritability for cause specific mortality based on twin studies.
Lifetime Data Analysis. http://dx.doi.org/10.1007/s10985-013-9244-x

BibTeX
@Article{
  title={Estimating heritability for cause specific mortality based on twin studies},
  author={Scheike, Thomas H. and Holst, Klaus K. and Hjelmborg, Jacob B.},
  year={2013},
  issn={1380-7870},
  journal={Lifetime Data Analysis},
  doi={10.1007/s10985-013-9244-x},
  url={http://dx.doi.org/10.1007/s10985-013-9244-x},
  publisher={Springer US},
  keywords={Cause specific hazards; Competing risks; Delayed entry; 
            Left truncation; Heritability; Survival analysis},
  pages={1-24},
  language={English}
}

Examples

library(mets)

data(prt) ## Prostate data example (sim)

## Bivariate competing risk, concordance estimates
p33 <- bicomprisk(Hist(time,status)~strata(zyg)+id(id),
                  data=prt,cause=c(2,2),return.data=1,prodlim=TRUE)

p33dz <- p33$model$"DZ"$comp.risk
p33mz <- p33$model$"MZ"$comp.risk

## Probability weights based on Aalen's additive model 
prtw <- ipw(Surv(time,status==0)~country, data=prt,
            cluster="id",weightname="w")


## Marginal model (wrongly ignoring censorings)
bpmz <- biprobit(cancer~1 + cluster(id), 
                 data=subset(prt,zyg=="MZ"), eqmarg=TRUE)

## Extended liability model
bpmzIPW <- biprobit(cancer~1 + cluster(id), 
                    data=subset(prtw,zyg=="MZ"), 
                    weight="w")
smz <- summary(bpmzIPW)

## Concordance
plot(p33mz,ylim=c(0,0.1),axes=FALSE,automar=FALSE,atrisk=FALSE)
axis(2); axis(1)

abline(h=smz$prob["Concordance",],lwd=c(2,1,1),col="darkblue")
## Wrong estimates:
abline(h=summary(bpmz)$prob["Concordance",],lwd=c(2,1,1),col="lightgray")

mets1.png

Latent Variable Models (lava)

Estimation and simulation of latent variable models

Installation

install.packages("lava",dependencies=TRUE)
library("lava")
demo("lava")

For graphical capabilities the Rgraphviz package is needed

source("http://bioconductor.org/biocLite.R");
biocLite("BiocUpgrade") ## Upgrade previous bioconductor installation
biocLite("Rgraphviz")

or the igraph package

install.packages("igraph",dependencies=TRUE)

The development version may be installed directly from R-forge:

install.packages("lava", repos="http://R-Forge.R-project.org")

Citation

To cite the lava package please use the following reference

Klaus K. Holst and Esben Budtz-Joergensen (2013). 
Linear Latent Variable Models: The lava-package. 
Computational Statistics 28 (4), pp 1385-1453. 
http://dx.doi.org/10.1007/s00180-012-0344-y

BibTeX:
@Article{lava,
  title = {Linear Latent Variable Models: The lava-package},
  author = {Klaus K. Holst and Esben Budtz-Joergensen},
  year = {2013},
  volume = {28},
  number = {4},
  pages = {1385-1452},
  journal = {Computational Statistics},
  note = {http://dx.doi.org/10.1007/s00180-012-0344-y},
}

Examples

## Specify structurual equation models with two factors
m <- lvm()
regression(m) <- c(y1,y2,y3) ~ eta1
regression(m) <- c(z1,z2,z3) ~ eta2
latent(m) <- ~eta1+eta2
regression(m) <- eta2~eta1+x
regression(m) <- eta1~x

labels(m) <- c(eta1=expression(eta[1]),eta2=expression(eta[2]))
plot(m)

lava1.png

## Simulation
set.seed(1)
d <- sim(m,100)
## Estimation
e <- estimate(m,d)
e
                    Estimate Std. Error  Z-value   P-value
Measurements:                                             
   y2~eta1           0.95462    0.08083 11.80993    <1e-12
   y3~eta1           0.98476    0.08922 11.03722    <1e-12
    z2~eta2          0.97038    0.05368 18.07714    <1e-12
    z3~eta2          0.95608    0.05643 16.94182    <1e-12
Regressions:                                              
   eta1~x            1.24587    0.11486 10.84694    <1e-12
    eta2~eta1        0.95608    0.18008  5.30910 1.102e-07
    eta2~x           1.11495    0.25228  4.41951 9.893e-06
Intercepts:                                               
   y2               -0.13896    0.12458 -1.11537    0.2647
   y3               -0.07661    0.13869 -0.55241    0.5807
   eta1              0.15801    0.12780  1.23644    0.2163
   z2               -0.00441    0.14858 -0.02969    0.9763
   z3               -0.15900    0.15731 -1.01076    0.3121
   eta2             -0.14143    0.18380 -0.76949    0.4416
Residual Variances:                                       
   y1                0.69684    0.14858  4.69004          
   y2                0.89804    0.16630  5.40026          
   y3                1.22456    0.21182  5.78109          
   eta1              0.93620    0.19623  4.77084          
   z1                1.41422    0.26259  5.38570          
   z2                0.87569    0.19463  4.49934          
   z3                1.18155    0.22640  5.21883          
   eta2              1.24430    0.28992  4.29195
## Assessing goodness-of-fit (linearity between eta2 and eta1)
library(gof)
g <- cumres(e,eta2~eta1)
plot(g)

gof1.png

Created: 2014-07-30 Wed 11:39