Analyzing twin data with mets
Table of Contents
- 1 Installation
- 2 Load simulated data
- 3 Estimation of cumulative incidence
- 4 Correcting for country
- 5 Concordance estimation
- 6 Effect of zygosity correcting for country
- 7 Liability model, ignoring censoring
- 8 Liability model, Inverse Probability Weighting
- 9 Liability model, adjusting for covariates
- 10 Cumulative heritability
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))
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)
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)
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))
### 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)
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)
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")
## 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")
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)