Analyzing twin data with mets

Table of Contents

1 Installation

Set repositories (see also chooseCRANmirror, chooseBioCmirror) and install dependencies (R >=2.14)

setRepositories() ## Choose CRAN and BioC Software (BioConductor)
install.packages(c("mets","cmprsk"))

OBS: At this point you might have to restart R to flush the cache of previously installed versions of the packages. If you have previously installed timereg and lava, make sure that you have the current versions installed (timereg: 1.6-5, lava: 1.0-5).

2 Load simulated data

library(mets)
data(prt)

3 Estimation of cumulative incidence

times <- seq(60,100,by=1)
cifmod <- comp.risk(Surv(time,status>0)~+1+cluster(id),data=prt,prt$status,causeS=2,n.sim=0,
                  times=times,conservative=1,max.clust=NULL,model="fg")

theta.des <- model.matrix(~-1+factor(zyg),data=prt) ## design for MZ/DZ status
or1 <- or.cif(cifmod,data=prt,cause1=2,cause2=2,theta.des=theta.des,same.cens=TRUE,
              score.method="fisher.scoring")
summary(or1)
or1$score

pcif <- predict(cifmod,X=1,resample.iid=0,uniform=0,se=0)
OR for dependence for competing risks

OR of cumulative incidence for cause1= 2  and cause2= 2
              log-ratio Coef.    SE     z    P-val Ratio    SE
factor(zyg)DZ            0.80 0.111  7.23 4.86e-13  2.22 0.246
factor(zyg)MZ            2.09 0.138 15.10 0.00e+00  8.07 1.110
             [,1]
[1,] 1.325225e-06
[2,] 3.458932e-06
plot(pcif,multiple=1,se=0,uniform=0,ylim=c(0,0.15))

pcif.png

4 Correcting for country

table(prt$country)

times <- seq(60,100,by=1)
cifmodl <-comp.risk(Surv(time,status>0)~-1+factor(country)+cluster(id),data=prt,
                    prt$status,causeS=2,n.sim=0,times=times,conservative=1,
                    max.clust=NULL,cens.model="aalen")
pcifl <- predict(cifmodl,X=diag(4),se=0,uniform=0)
plot(pcifl,multiple=1,se=0,uniform=0,col=1:4,ylim=c(0,0.2))
legend("topleft",levels(prt$country),col=1:4,lty=1)

pcifl.png

theta.des <- model.matrix(~-1+factor(zyg),data=prt) ## design for MZ/DZ status
or.country <- or.cif(cifmodl,data=prt,cause1=2,cause2=2,theta.des=theta.des,same.cens=TRUE,
                     theta=c(2.8,6.9),score.method="fisher.scoring")
summary(or.country)
or.country$score
OR for dependence for competing risks

OR of cumulative incidence for cause1= 2  and cause2= 2
              log-ratio Coef.    SE     z    P-val Ratio    SE
factor(zyg)DZ           0.754 0.117  6.43 1.26e-10  2.12 0.249
factor(zyg)MZ           1.850 0.139 13.30 0.00e+00  6.36 0.883
              [,1]
[1,] -1.201999e-06
[2,]  1.558011e-06
cifmodlr <-comp.risk(Surv(time,status>0)~+1+const(factor(country))+cluster(id),data=prt,
                    prt$status,causeS=2,n.sim=0,times=times,conservative=1,max.clust=NULL,model="fg",
                    cens.model="aalen",cens.formula=~~factor(country))
pciflr <- predict(cifmodlr,X=rep(1,4),Z=rbind(c(0,0,0),c(1,0,0),c(0,1,0),c(0,0,1)),se=0,uniform=0)
par(mfrow=c(1,2))
plot(pcifl,multiple=1,se=0,uniform=0,col=1:4,ylim=c(0,0.2))
legend("topleft",levels(prt$country),col=1:4,lty=1)
plot(pciflr,multiple=1,se=0,uniform=0,col=1:4,ylim=c(0,0.2))
legend("topleft",levels(prt$country),col=1:4,lty=1)

pcif2.png

or.countryr <- or.cif(cifmodlr,data=prt,cause1=2,cause2=2,theta.des=theta.des,same.cens=TRUE,
                     theta=c(2.8,6.9),score.method="fisher.scoring")
summary(or.countryr)
OR for dependence for competing risks

OR of cumulative incidence for cause1= 2  and cause2= 2
              log-ratio Coef.    SE     z    P-val Ratio    SE
factor(zyg)DZ           0.756 0.117  6.48 9.33e-11  2.13 0.249
factor(zyg)MZ           1.850 0.139 13.40 0.00e+00  6.38 0.886

5 Concordance estimation

### ignoring country 
p33 <- bicomprisk(Hist(time,status)~strata(zyg)+id(id),data=prt,cause=c(2,2),return.data=1,robust=1)

p33dz <- p33$model$"DZ"$comp.risk
p33mz <- p33$model$"MZ"$comp.risk
plot(p33dz,se=0,ylim=c(0,0.1))
lines(p33mz$time,p33mz$P1,col=3)
title(main="Concordance Prostate cancer")
lines(pcif$time,pcif$P1^2,col=2)
legend("topleft",c("DZ","MZ","Independence"),lty=rep(1,3),col=c(1,3,2))

p33dz.png

### test for genetic effect 
test.conc(p33dz,p33mz); 

Pepe-Mori type test for H_0: conc_1(t)= conc_2(t)
Assuming independence for estimators
Time.range = 60.9 -- 96.9 

          cum dif.    sd    z     pval
pepe-mori    0.394 0.095 4.15 3.39e-05
data33mz <- p33$model$"MZ"$data
data33mz$zyg <- 1
data33dz <- p33$model$"DZ"$data
data33dz$zyg <- 0
data33 <- rbind(data33mz,data33dz)

library(cmprsk)
ftime <- data33$time
fstatus <- data33$status
table(fstatus)
fstatus
   0    1    2 
9597  106 4519
group <- data33$zyg
graytest <- cuminc(ftime,fstatus,group)
graytest
Tests:
      stat           pv df
1 28.82416 7.925617e-08  1
2 33.79236 6.131919e-09  1
Estimates and Variances:
$est
              20         40           60          80        100
0 1 0.0000000000 0.00000000 0.0001741916 0.006741025 0.01880244
1 1 0.0000000000 0.00000000 0.0006710172 0.017420360 0.05031415
0 2 0.0006970762 0.01974882 0.1141800067 0.504364854 0.93797293
1 2 0.0009363302 0.01655314 0.0948098327 0.443996722 0.90692430

$var
              20           40           60           80          100
0 1 0.000000e+00 0.000000e+00 3.034323e-08 2.115863e-06 9.493584e-06
1 1 0.000000e+00 0.000000e+00 2.250627e-07 9.173278e-06 5.102841e-05
0 2 8.094463e-08 2.487399e-06 1.556735e-05 6.990685e-05 4.769058e-05
1 2 1.752378e-07 3.424511e-06 2.388136e-05 1.271394e-04 1.171775e-04
zygeffect <- comp.risk(Surv(time,status==0)~const(zyg),
                  data=data33,data33$status,causeS=1,
                  cens.model="aalen",model="logistic",conservative=1)
summary(zygeffect)
Competing risks Model 

Test for nonparametric terms 

Test for non-significant effects 
            Supremum-test of significance p-value H_0: B(t)=0
(Intercept)                          25.5                   0

Test for time invariant effects 
                  Kolmogorov-Smirnov test p-value H_0:constant effect
(Intercept)                          2.23                           0
                    Cramer von Mises test p-value H_0:constant effect
(Intercept)                          36.2                           0

Parametric terms : 
           Coef.   SE Robust SE    z    P-val
const(zyg) 0.977 0.22      0.22 4.44 9.06e-06
   
  Call: 
comp.risk(Surv(time, status == 0) ~ const(zyg), data = data33, 
    data33$status, causeS = 1, cens.model = "aalen", model = "logistic", 
    conservative = 1)
case33mz <- conc2case(p33mz,pcif)
case33dz <- conc2case(p33dz,pcif)

plot(case33mz$casewise,se=0,col=3)
lines(case33dz$casewise$time,case33dz$casewise$P1)
title(main="Probandwise concordance")
legend("topleft",c("MZ","DZ","Independence"),lty=rep(1,3),col=c(3,1,2))
lines(pcif$time,pcif$P1,col=2)

casewise.png

6 Effect of zygosity correcting for country

p33l <- bicomprisk(Hist(time,status)~country+strata(zyg)+id(id),
                data=prt,cause=c(2,2),return.data=1,robust=1)

data33mz <- p33l$model$"MZ"$data
data33mz$zyg <- 1
data33dz <- p33l$model$"DZ"$data
data33dz$zyg <- 0
data33 <- rbind(data33mz,data33dz)
zygeffectl <- comp.risk(Surv(time,status==0)~const(country)+const(zyg),
                  data=data33,data33$status,causeS=1,
                  cens.model="aalen",model="logistic",conservative=1)
summary(zygeffectl)
Competing risks Model 

Test for nonparametric terms 

Test for non-significant effects 
            Supremum-test of significance p-value H_0: B(t)=0
(Intercept)                          16.1                   0

Test for time invariant effects 
                  Kolmogorov-Smirnov test p-value H_0:constant effect
(Intercept)                          2.01                           0
                    Cramer von Mises test p-value H_0:constant effect
(Intercept)                          35.9                           0

Parametric terms : 
                      Coef.    SE Robust SE    z    P-val
const(country)Finland 1.160 0.419     0.419 2.77 5.54e-03
const(country)Norway  0.655 0.458     0.458 1.43 1.53e-01
const(country)Sweden  0.796 0.372     0.372 2.14 3.23e-02
const(zyg)            0.932 0.230     0.230 4.05 5.15e-05
   
  Call: 
comp.risk(Surv(time, status == 0) ~ const(country) + const(zyg), 
    data = data33, data33$status, causeS = 1, cens.model = "aalen", 
    model = "logistic", conservative = 1)
zygeffectpl <- comp.risk(Surv(time,status==0)~const(country)+const(zyg),
                  data=data33,data33$status,causeS=1,
                  cens.model="aalen",model="fg",conservative=1)
print(summary(zygeffectpl))
Competing risks Model 

Test for nonparametric terms 

Test for non-significant effects 
            Supremum-test of significance p-value H_0: B(t)=0
(Intercept)                          2.83               0.012

Test for time invariant effects 
                  Kolmogorov-Smirnov test p-value H_0:constant effect
(Intercept)                        0.0101                           0
                    Cramer von Mises test p-value H_0:constant effect
(Intercept)                       0.00115                       0.004

Parametric terms : 
                      Coef.    SE Robust SE    z    P-val
const(country)Finland 1.140 0.412     0.412 2.77 5.63e-03
const(country)Norway  0.646 0.452     0.452 1.43 1.53e-01
const(country)Sweden  0.785 0.368     0.368 2.14 3.27e-02
const(zyg)            0.916 0.226     0.226 4.05 5.22e-05
   
  Call: 
comp.risk(Surv(time, status == 0) ~ const(country) + const(zyg), 
    data = data33, data33$status, causeS = 1, cens.model = "aalen", 
    model = "fg", conservative = 1)

NULL
zygeffectll <- comp.risk(Surv(time,status==0)~country+const(zyg),
                         data=data33,data33$status,causeS=1,
                         cens.model="aalen",model="logistic",conservative=1)
print(summary(zygeffectll))
Competing risks Model 

Test for nonparametric terms 

Test for non-significant effects 
               Supremum-test of significance p-value H_0: B(t)=0
(Intercept)                            75.70                   0
countryFinland                        441.00                   0
countryNorway                           6.09                   0
countrySweden                         703.00                   0

Test for time invariant effects 
                     Kolmogorov-Smirnov test p-value H_0:constant effect
(Intercept)                             6.59                       0.000
countryFinland                          6.24                       0.000
countryNorway                           1.31                       0.574
countrySweden                           6.39                       0.000
                       Cramer von Mises test p-value H_0:constant effect
(Intercept)                            200.0                         0.0
countryFinland                        1180.0                         0.0
countryNorway                           17.6                         0.4
countrySweden                         1300.0                         0.0

Parametric terms : 
           Coef.   SE Robust SE    z    P-val
const(zyg) 0.939 0.23      0.23 4.08 4.58e-05
   
WARNING problem with convergence for time points:
64.88587 66.74123
Readjust analyses by removing points

  Call: 
comp.risk(Surv(time, status == 0) ~ country + const(zyg), data = data33, 
    data33$status, causeS = 1, cens.model = "aalen", model = "logistic", 
    conservative = 1)

NULL

7 Liability model, ignoring censoring

(M <- with(prt, table(cancer,zyg)))
      zyg
cancer    DZ    MZ
     0 17408 10872
     1   583   359
coef(lm(cancer~-1+zyg,prt))
     zygDZ      zygMZ 
0.03240509 0.03196510
## Saturated model
bpmz <- 
    biprobit(cancer~1 + cluster(id), 
             data=subset(prt,zyg=="MZ"), eqmarg=TRUE)

logLik(bpmz) # Log-likelihood
AIC(bpmz) # AIC
coef(bpmz) # Parameter estimates
vcov(bpmz) # Asymptotic covariance
summary(bpmz) # concordance, case-wise, tetrachoric correlations, ...
'log Lik.' -1472.972 (df=2)
[1] 2949.943
(Intercept)  atanh(rho) 
 -1.8539454   0.8756506
             (Intercept)   atanh(rho)
(Intercept) 0.0007089726 0.0003033296
atanh(rho)  0.0003033296 0.0044023587

              Estimate    Std.Err          Z p-value
(Intercept)  -1.853945   0.026627 -69.627727       0
atanh(rho)    0.875651   0.066350  13.197393       0

    n pairs 
11231  5473 
Score: -3.453e-05 5.123e-06
logLik: -1472.972 
Variance of latent residual term = 1 (standard probit link) 

                        Estimate 2.5%    97.5%  
Tetrachoric correlation 0.70423  0.63252 0.76398
Concordance             0.01131  0.00886 0.01443
Case-wise/Conditional   0.35487  0.29391 0.42094
Marginal                0.03187  0.02834 0.03583
bp0 <- biprobit(cancer~1 + cluster(id)+strata(zyg), data=prt)
summary(bp0)
------------------------------------------------------------
Strata 'DZ'

              Estimate    Std.Err          Z p-value
(Intercept)  -1.846841   0.019247 -95.955243       0
atanh(rho)    0.418065   0.050421   8.291446       0

    n pairs 
17991  8749 
Score: -0.001841 -0.0006879
logLik: -2536.242 
Variance of latent residual term = 1 (standard probit link) 

                        Estimate 2.5%    97.5%  
Tetrachoric correlation 0.39530  0.30882 0.47529
Concordance             0.00486  0.00361 0.00655
Case-wise/Conditional   0.15019  0.11459 0.19443
Marginal                0.03239  0.02976 0.03523

------------------------------------------------------------
Strata 'MZ'

              Estimate    Std.Err          Z p-value
(Intercept)  -1.853945   0.026627 -69.627727       0
atanh(rho)    0.875651   0.066350  13.197393       0

    n pairs 
11231  5473 
Score: -3.453e-05 5.123e-06
logLik: -1472.972 
Variance of latent residual term = 1 (standard probit link) 

                        Estimate 2.5%    97.5%  
Tetrachoric correlation 0.70423  0.63252 0.76398
Concordance             0.01131  0.00886 0.01443
Case-wise/Conditional   0.35487  0.29391 0.42094
Marginal                0.03187  0.02834 0.03583
## Eq. marginals MZ/DZ
bp1 <- bptwin(cancer~1,zyg="zyg",DZ="DZ",id="id",type="u",data=prt)
summary(bp1) # Components (concordance,cor,...) can be extracted from returned list
                 Estimate     Std.Err           Z p-value
(Intercept)     -1.849284    0.015601 -118.539777       0
atanh(rho) MZ    0.877667    0.065815   13.335456       0
atanh(rho) DZ    0.417475    0.050276    8.303615       0

 Total MZ/DZ Complete pairs MZ/DZ
 11231/17991 5473/8749           

                           Estimate 2.5%    97.5%  
Tetrachoric correlation MZ 0.70525  0.63436 0.76438
Tetrachoric correlation DZ 0.39480  0.30854 0.47462

MZ:
                        Estimate 2.5%    97.5%  
Concordance             0.01149  0.00942 0.01400
Probandwise Concordance 0.35672  0.29764 0.42049
Marginal                0.03221  0.03007 0.03449
DZ:
                        Estimate 2.5%    97.5%  
Concordance             0.00482  0.00363 0.00640
Probandwise Concordance 0.14956  0.11441 0.19315
Marginal                0.03221  0.03007 0.03449

                         Estimate 2.5%    97.5%  
Broad-sense Heritability 0.62090  0.40145 0.79997
compare(bp0,bp1) # LRT

      Likelihood ratio test

data:  
chisq = 0.0468, df = 1, p-value = 0.8288
sample estimates:
log likelihood (model 1) log likelihood (model 2) 
               -4009.213                -4009.237

Polygenic Libability model via te bptwin function (type can be a subset of "acde", or "flex" for stratitified, "u" for random effects model with same marginals for MZ and DZ)

## Polygenic model
args(bptwin)
function (formula, data, id, zyg, DZ, OS, weight = NULL, biweight = function(x) 1/min(x), 
    strata = NULL, messages = 1, control = list(trace = 0), type = "ace", 
    eqmean = TRUE, pairsonly = FALSE, samecens = TRUE, allmarg = samecens & 
        !is.null(weight), stderr = TRUE, robustvar = TRUE, p, 
    indiv = FALSE, constrain, bound = FALSE, debug = FALSE, ...) 
NULL
bp2 <- bptwin(cancer~1,zyg="zyg",DZ="DZ",id="id",type="ace",data=prt)
summary(bp2)
             Estimate   Std.Err         Z p-value
(Intercept)  -3.40624   0.19032 -17.89736  0.0000
log(var(A))   0.74503   0.25710   2.89787  0.0038
log(var(C))  -1.25112   1.04238  -1.20024  0.2300

 Total MZ/DZ Complete pairs MZ/DZ
 11231/17991 5473/8749           

                   Estimate 2.5%    97.5%  
A                  0.62090  0.40145 0.79997
C                  0.08435  0.00910 0.48028
E                  0.29475  0.23428 0.36343
MZ Tetrachoric Cor 0.70525  0.63436 0.76438
DZ Tetrachoric Cor 0.39480  0.30854 0.47462

MZ:
                        Estimate 2.5%    97.5%  
Concordance             0.01149  0.00942 0.01400
Probandwise Concordance 0.35672  0.29764 0.42049
Marginal                0.03221  0.03007 0.03449
DZ:
                        Estimate 2.5%    97.5%  
Concordance             0.00482  0.00363 0.00640
Probandwise Concordance 0.14956  0.11441 0.19315
Marginal                0.03221  0.03007 0.03449

                         Estimate 2.5%    97.5%  
Broad-sense Heritability 0.70525  0.63657 0.76572

8 Liability model, Inverse Probability Weighting

## Probability weights based on Aalen's additive model 
prtw <- ipw(Surv(time,status==0)~country, data=prt,
            cluster="id",weightname="w") 
plot(0,type="n",xlim=range(prtw$time),ylim=c(0,1),xlab="Age",ylab="Probability")
count <- 0
for (l in unique(prtw$country)) {
    count <- count+1
    prtw <- prtw[order(prtw$time),]
    with(subset(prtw,country==l), 
         lines(time,w,col=count,lwd=2))
}
legend("topright",legend=unique(prtw$country),col=1:4,pch=1)

ipw.png

bpmzIPW <- 
              biprobit(cancer~1 + cluster(id), 
                       data=subset(prtw,zyg=="MZ"), 
                       weight="w")
(smz <- summary(bpmzIPW))
              Estimate    Std.Err          Z p-value
(Intercept)  -1.226276   0.043074 -28.469378       0
atanh(rho)    0.912670   0.100316   9.097911       0

    n pairs 
 2722   997 
Score: 3.318e-05 -2.252e-05
logLik: -6703.246 
Variance of latent residual term = 1 (standard probit link) 

                        Estimate 2.5%    97.5%  
Tetrachoric correlation 0.72241  0.61446 0.80381
Concordance             0.05490  0.04221 0.07113
Case-wise/Conditional   0.49887  0.41321 0.58460
Marginal                0.11005  0.09514 0.12696
## CIF
plot(pcif,multiple=1,se=0,uniform=0,ylim=c(0,0.15))
abline(h=smz$prob["Marginal",],lwd=c(2,1,1))
## Wrong estimates:
abline(h=summary(bpmz)$prob["Marginal",],lwd=c(2,1,1),col="lightgray")

cif2.png

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

conc2.png

bp3 <- bptwin(cancer~1,zyg="zyg",DZ="DZ",id="id",
              type="ace",data=prtw,weight="w")
summary(bp3)
             Estimate   Std.Err         Z p-value
(Intercept)  -2.31618   0.18673 -12.40359   0e+00
log(var(A))   0.85390   0.22689   3.76347   2e-04
log(var(C)) -29.43218   1.13343 -25.96726   0e+00

 Total MZ/DZ Complete pairs MZ/DZ
 2722/5217   997/1809            

                   Estimate 2.5%    97.5%  
A                  0.70138  0.60090 0.78560
C                  0.00000  0.00000 0.00000
E                  0.29862  0.21440 0.39910
MZ Tetrachoric Cor 0.70138  0.59586 0.78310
DZ Tetrachoric Cor 0.35069  0.30328 0.39637

MZ:
                        Estimate 2.5%    97.5%  
Concordance             0.04857  0.03963 0.05940
Probandwise Concordance 0.47238  0.39356 0.55260
Marginal                0.10281  0.09463 0.11161
DZ:
                        Estimate 2.5%    97.5%  
Concordance             0.02515  0.02131 0.02965
Probandwise Concordance 0.24461  0.21892 0.27226
Marginal                0.10281  0.09463 0.11161

                         Estimate 2.5%    97.5%  
Broad-sense Heritability 0.70138  0.60090 0.78560
bp4 <- bptwin(cancer~1,zyg="zyg",DZ="DZ",id="id",
              type="u",data=prtw,weight="w")
summary(bp4)
                Estimate    Std.Err          Z p-value
(Intercept)    -1.266427   0.024091 -52.568381       0
atanh(rho) MZ   0.898548   0.098841   9.090866       0
atanh(rho) DZ   0.312574   0.073668   4.243006       0

 Total MZ/DZ Complete pairs MZ/DZ
 2722/5217   997/1809            

                           Estimate 2.5%    97.5%  
Tetrachoric correlation MZ 0.71559  0.60742 0.79771
Tetrachoric correlation DZ 0.30278  0.16662 0.42760

MZ:
                        Estimate 2.5%    97.5%  
Concordance             0.04974  0.04044 0.06104
Probandwise Concordance 0.48442  0.40185 0.56785
Marginal                0.10268  0.09453 0.11144
DZ:
                        Estimate 2.5%    97.5%  
Concordance             0.02269  0.01667 0.03081
Probandwise Concordance 0.22097  0.16448 0.29013
Marginal                0.10268  0.09453 0.11144

                         Estimate 2.5%    97.5%  
Broad-sense Heritability 0.82563  0.33329 0.97819
score(bp4) ## Check convergence
[1]  2.729971e-07 -8.463577e-08 -5.014015e-09
bp5 <- bptwin(cancer~1,zyg="zyg",DZ="DZ",id="id",
              type="ade",data=prtw,weight="w")
summary(bp5)
             Estimate   Std.Err         Z p-value
(Intercept)  -2.37470   0.20268 -11.71665  0.0000
log(var(A))   0.55519   0.54480   1.01905  0.3082
log(var(D))  -0.25645   1.36092  -0.18844  0.8505

 Total MZ/DZ Complete pairs MZ/DZ
 2722/5217   997/1809            

                   Estimate 2.5%    97.5%  
A                  0.49552  0.10422 0.89238
D                  0.22007  0.01081 0.87931
E                  0.28441  0.19987 0.38740
MZ Tetrachoric Cor 0.71559  0.60742 0.79771
DZ Tetrachoric Cor 0.30278  0.16662 0.42760

MZ:
                        Estimate 2.5%    97.5%  
Concordance             0.04974  0.04044 0.06104
Probandwise Concordance 0.48442  0.40185 0.56785
Marginal                0.10268  0.09453 0.11144
DZ:
                        Estimate 2.5%    97.5%  
Concordance             0.02269  0.01667 0.03081
Probandwise Concordance 0.22097  0.16448 0.29013
Marginal                0.10268  0.09453 0.11144

                         Estimate 2.5%    97.5%  
Broad-sense Heritability 0.71559  0.61260 0.80013

9 Liability model, adjusting for covariates

Main effect of country

bp6 <- bptwin(cancer~country,zyg="zyg",DZ="DZ",id="id",
              type="ace",data=prtw,weight="w")
summary(bp6)
Warning message:
In sqrt(diag(V)) : NaNs produced

                Estimate   Std.Err         Z p-value
(Intercept)     -2.81553   0.23889 -11.78590   0e+00
countryFinland   0.87558   0.16123   5.43061   0e+00
countryNorway    0.68483   0.17762   3.85567   1e-04
countrySweden    0.77248   0.12350   6.25468   0e+00
log(var(A))      0.77724   0.23186   3.35220   8e-04
log(var(C))    -28.96268        NA        NA      NA

 Total MZ/DZ Complete pairs MZ/DZ
 2722/5217   997/1809            

                   Estimate 2.5%    97.5%  
A                  0.68509  0.58001 0.77411
C                  0.00000  0.00000 0.00000
E                  0.31491  0.22589 0.41999
MZ Tetrachoric Cor 0.68509  0.57428 0.77124
DZ Tetrachoric Cor 0.34254  0.29262 0.39060

MZ:
                        Estimate 2.5%    97.5%  
Concordance             0.02236  0.01588 0.03141
Probandwise Concordance 0.39194  0.30778 0.48305
Marginal                0.05705  0.04654 0.06977
DZ:
                        Estimate 2.5%    97.5%  
Concordance             0.00989  0.00700 0.01394
Probandwise Concordance 0.17329  0.14505 0.20570
Marginal                0.05705  0.04654 0.06977

                         Estimate 2.5%    97.5%  
Broad-sense Heritability 0.68509  0.58001 0.77411

Stratified analysis

bp7 <- bptwin(cancer~country,zyg="zyg",DZ="DZ",id="id",
              type="u",data=prtw,weight="w")
summary(bp7)
                 Estimate    Std.Err          Z p-value
(Intercept)     -1.581478   0.051318 -30.817030   0e+00
countryFinland   0.491725   0.081517   6.032155   0e+00
countryNorway    0.385830   0.094254   4.093497   0e+00
countrySweden    0.433789   0.060648   7.152599   0e+00
atanh(rho) MZ    0.884166   0.099366   8.898113   0e+00
atanh(rho) DZ    0.271770   0.073240   3.710668   2e-04

 Total MZ/DZ Complete pairs MZ/DZ
 2722/5217   997/1809            

                           Estimate 2.5%    97.5%  
Tetrachoric correlation MZ 0.70850  0.59760 0.79280
Tetrachoric correlation DZ 0.26527  0.12752 0.39298

MZ:
                        Estimate 2.5%    97.5%  
Concordance             0.02347  0.01664 0.03300
Probandwise Concordance 0.41255  0.32395 0.50721
Marginal                0.05688  0.04643 0.06953
DZ:
                        Estimate 2.5%    97.5%  
Concordance             0.00794  0.00489 0.01287
Probandwise Concordance 0.13966  0.09312 0.20421
Marginal                0.05688  0.04643 0.06953

                         Estimate 2.5%    97.5%  
Broad-sense Heritability 0.88646  0.22665 0.99521
bp8 <- bptwin(cancer~strata(country),zyg="zyg",DZ="DZ",id="id",
              type="u",data=prtw,weight="w")
summary(bp8)
------------------------------------------------------------
Strata 'Denmark'

                Estimate    Std.Err          Z p-value
(Intercept)    -1.583608   0.051241 -30.904856  0.0000
atanh(rho) MZ   0.992896   0.217349   4.568215  0.0000
atanh(rho) DZ   0.070588   0.186956   0.377566  0.7058

 Total MZ/DZ Complete pairs MZ/DZ
 760/1611    287/589             

                           Estimate 2.5%     97.5%   
Tetrachoric correlation MZ  0.75859  0.51308  0.88937
Tetrachoric correlation DZ  0.07047 -0.28750  0.41117

MZ:
                        Estimate 2.5%    97.5%  
Concordance             0.02611  0.01584 0.04274
Probandwise Concordance 0.46093  0.28426 0.64799
Marginal                0.05664  0.04623 0.06922
DZ:
                        Estimate 2.5%    97.5%  
Concordance             0.00420  0.00110 0.01596
Probandwise Concordance 0.07422  0.01888 0.25037
Marginal                0.05664  0.04623 0.06922

                         Estimate 2.5% 97.5%
Broad-sense Heritability   1      NaN  NaN  

------------------------------------------------------------
Strata 'Finland'

                Estimate    Std.Err          Z p-value
(Intercept)    -1.087902   0.063221 -17.207912  0.0000
atanh(rho) MZ   0.859335   0.302752   2.838410  0.0045
atanh(rho) DZ   0.393145   0.179942   2.184840  0.0289

 Total MZ/DZ Complete pairs MZ/DZ
 392/1001    134/316             

                           Estimate 2.5%    97.5%  
Tetrachoric correlation MZ 0.69592  0.25985 0.89623
Tetrachoric correlation DZ 0.37407  0.04044 0.63265

MZ:
                        Estimate 2.5%    97.5%  
Concordance             0.07008  0.03975 0.12064
Probandwise Concordance 0.50666  0.27641 0.73412
Marginal                0.13832  0.11316 0.16801
DZ:
                        Estimate 2.5%    97.5%  
Concordance             0.04160  0.02237 0.07607
Probandwise Concordance 0.30073  0.16558 0.48242
Marginal                0.13832  0.11316 0.16801

                         Estimate 2.5%    97.5%  
Broad-sense Heritability 0.64369  0.04069 0.98717

------------------------------------------------------------
Strata 'Norway'

                Estimate    Std.Err          Z p-value
(Intercept)    -1.192293   0.079124 -15.068598  0.0000
atanh(rho) MZ   0.916471   0.301133   3.043409  0.0023
atanh(rho) DZ   0.533761   0.252070   2.117509  0.0342

 Total MZ/DZ Complete pairs MZ/DZ
 387/618     115/155             

                           Estimate 2.5%    97.5%  
Tetrachoric correlation MZ 0.72422  0.31516 0.90635
Tetrachoric correlation DZ 0.48825  0.03969 0.77303

MZ:
                        Estimate 2.5%    97.5%  
Concordance             0.05918  0.03218 0.10633
Probandwise Concordance 0.50764  0.27633 0.73572
Marginal                0.11657  0.08945 0.15057
DZ:
                        Estimate 2.5%    97.5%  
Concordance             0.03945  0.01840 0.08257
Probandwise Concordance 0.33842  0.15583 0.58636
Marginal                0.11657  0.08945 0.15057

                         Estimate 2.5%    97.5%  
Broad-sense Heritability 0.47195  0.01989 0.97522

------------------------------------------------------------
Strata 'Sweden'

                Estimate    Std.Err          Z p-value
(Intercept)    -1.149412   0.032155 -35.745836  0.0000
atanh(rho) MZ   0.836864   0.125476   6.669520  0.0000
atanh(rho) DZ   0.199677   0.092907   2.149202  0.0316

 Total MZ/DZ Complete pairs MZ/DZ
 1183/1987   461/749             

                           Estimate 2.5%    97.5%  
Tetrachoric correlation MZ 0.68414  0.53057 0.79423
Tetrachoric correlation DZ 0.19706  0.01758 0.36425

MZ:
                        Estimate 2.5%    97.5%  
Concordance             0.06055  0.04659 0.07835
Probandwise Concordance 0.48365  0.38001 0.58872
Marginal                0.12519  0.11277 0.13877
DZ:
                        Estimate 2.5%    97.5%  
Concordance             0.02515  0.01672 0.03766
Probandwise Concordance 0.20088  0.13541 0.28746
Marginal                0.12519  0.11277 0.13877

                         Estimate 2.5%    97.5%  
Broad-sense Heritability 0.97416  0.00000 1.00000
## Wald test
B <- (lava::contrmat(3,4))[-(1:3),]
compare(bp8,contrast=B)

      Wald test

data:  
chisq = 3.4972, df = 6, p-value = 0.7443

10 Cumulative heritability

args(cumh)
function (formula, data, ..., time, timestrata = quantile(data[, 
    time], c(0.25, 0.5, 0.75, 1)), cumulative = TRUE, silent = FALSE) 
NULL
ch1 <- cumh(cancer~1,time="time",zyg="zyg",DZ="DZ",id="id",
            type="ace",data=prtw,weight="w")
summary(ch1)
                      time Heritability    Std.Err      2.5%     97.5%
65.5691955406266  65.56920    0.7038286 0.10969626 0.4586422 0.8695520
76.4446739437236  76.44467    0.6757445 0.06363443 0.5411756 0.7864218
85.8807708995545  85.88077    0.6204174 0.05652481 0.5052219 0.7234726
117.622394945129 117.62239    0.7013847 0.04752116 0.6008962 0.7855993
plot(ch1)

cumh.png


Date: 2012-05-20

Author: Klaus K. Holst and Thomas Scheike