When looking at multivariate binomial data with the aim of learning about the dependence that is present, possibly after correcting for some covariates many models are available.
in the mets package you can fit the
Pairwise odds ratio model
Bivariate Probit model
Additive gamma random effects model
These last three models are all fitted in the mets package using composite likelihoods for pairs of data. The models can be fitted specifically based on specifying which pairs one wants to use for the composite score.
The models are described in futher details in the binomial-twin vignette.
We start by simulating family data with and additive gamma structure on ACE form. Here 10000 families consisting of two parents and two children. The response is ybin and there is one covariate x.
library(mets)
set.seed(100)
data <- simbinClaytonOakes.family.ace(1000,2,1,beta=NULL,alpha=NULL)
data$number <- c(1,2,3,4)
data$child <- 1*(data$number==3)
head(data)
#> ybin x type cluster number child
#> 1 1 1 mother 1 1 0
#> 2 0 1 father 1 2 0
#> 3 1 1 child 1 3 1
#> 4 0 1 child 1 4 0
#> 5 1 0 mother 2 1 0
#> 6 0 0 father 2 2 0We fit the marginal models, and here find a covariate effect at 0.3 for x. The marginals can be specified excatly as one wants.
aa <- margbin <- glm(ybin~x,data=data,family=binomial())
summary(aa)
#>
#> Call:
#> glm(formula = ybin ~ x, family = binomial(), data = data)
#>
#> Deviance Residuals:
#> Min 1Q Median 3Q Max
#> -1.5202 -1.4036 0.8697 0.9669 0.9669
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 0.51765 0.04676 11.07 <2e-16 ***
#> x 0.25959 0.06673 3.89 1e-04 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 5146.6 on 3999 degrees of freedom
#> Residual deviance: 5131.5 on 3998 degrees of freedom
#> AIC: 5135.5
#>
#> Number of Fisher Scoring iterations: 4For the additive gamma of this type we set-up the random effects included in such a family to make the ACE valid using some special functions for this.
The model is constructe with one enviromental effect shared by all in the family and 8 genetic random effects with size (1/4) genetic variance. Looking at the first family we see that the mother and father both share half the genes with the children and that the two children also share half their genes with this specification. Below we also show an alternative specification of this model using all pairs.
# make ace random effects design
out <- ace.family.design(data,member="type",id="cluster")
out$pardes
#> [,1] [,2]
#> [1,] 0.25 0
#> [2,] 0.25 0
#> [3,] 0.25 0
#> [4,] 0.25 0
#> [5,] 0.25 0
#> [6,] 0.25 0
#> [7,] 0.25 0
#> [8,] 0.25 0
#> [9,] 0.00 1
head(out$des.rv,4)
#> m1 m2 m3 m4 f1 f2 f3 f4 env
#> [1,] 1 1 1 1 0 0 0 0 1
#> [2,] 0 0 0 0 1 1 1 1 1
#> [3,] 1 1 0 0 1 1 0 0 1
#> [4,] 1 0 1 0 1 0 1 0 1We can now fit the model calling the two-stage function
# fitting ace model for family structure
ts <- binomial.twostage(margbin,data=data,clusters=data$cluster,
theta=c(2,1),
random.design=out$des.rv,theta.des=out$pardes)
summary(ts)
#> Dependence parameter for Clayton-Oakes model
#> Variance of Gamma distributed random effects
#> $estimates
#> theta se
#> dependence1 2.5726358 0.6212573
#> dependence2 0.8903891 0.2487322
#>
#> $type
#> [1] "clayton.oakes"
#>
#> $h
#> Estimate Std.Err 2.5% 97.5% P-value
#> dependence1 0.7429 0.07228 0.6012 0.8845 8.833e-25
#> dependence2 0.2571 0.07228 0.1155 0.3988 3.747e-04
#>
#> $vare
#> NULL
#>
#> $vartot
#> Estimate Std.Err 2.5% 97.5% P-value
#> p1 3.463 0.6574 2.175 4.751 1.381e-07
#>
#> attr(,"class")
#> [1] "summary.mets.twostage"
# true variance parameters
c(2,1)
#> [1] 2 1
# total variance
3
#> [1] 3We now specify the same model via extracting all pairs.
The random effecs structure is typically simpler when just looking at
pairs, but we start by specifying the random effects jointly for whole
family. A special function writes up all combinations of pairs. There
are 6 pairs within each family, and we keep track of who belongs to the
different families. We first simply give the pairs and we then should
get the same result as before.
mm <- familycluster.index(data$cluster)
head(mm$familypairindex,n=20)
#> [1] 1 2 1 3 1 4 2 3 2 4 3 4 5 6 5 7 5 8 6 7
pairs <- mm$pairs
dim(pairs)
#> [1] 6000 2
head(pairs,12)
#> [,1] [,2]
#> [1,] 1 2
#> [2,] 1 3
#> [3,] 1 4
#> [4,] 2 3
#> [5,] 2 4
#> [6,] 3 4
#> [7,] 5 6
#> [8,] 5 7
#> [9,] 5 8
#> [10,] 6 7
#> [11,] 6 8
#> [12,] 7 8Now with the pairs we fit the model
tsp <- binomial.twostage(margbin,data=data,clusters=data$cluster,theta=c(2,1),detail=0,
random.design=out$des.rv,theta.des=out$pardes,pairs=pairs)
summary(tsp)
#> Dependence parameter for Clayton-Oakes model
#> Variance of Gamma distributed random effects
#> $estimates
#> theta se
#> dependence1 2.5726358 0.6212573
#> dependence2 0.8903891 0.2487322
#>
#> $type
#> [1] "clayton.oakes"
#>
#> $h
#> Estimate Std.Err 2.5% 97.5% P-value
#> dependence1 0.7429 0.07228 0.6012 0.8845 8.833e-25
#> dependence2 0.2571 0.07228 0.1155 0.3988 3.747e-04
#>
#> $vare
#> NULL
#>
#> $vartot
#> Estimate Std.Err 2.5% 97.5% P-value
#> p1 3.463 0.6574 2.175 4.751 1.381e-07
#>
#> attr(,"class")
#> [1] "summary.mets.twostage"Here a random sample of pairs are given instead and we get other estimates.
set.seed(100)
ssid <- sort(sample(1:nrow(pairs),nrow(pairs)/2))
tsd <- binomial.twostage(aa,data=data,clusters=data$cluster,
theta=c(2,1),step=1.0,
random.design=out$des.rv,iid=1,Nit=10,
theta.des=out$pardes,pairs=pairs[ssid,])
summary(tsd)
#> Dependence parameter for Clayton-Oakes model
#> Variance of Gamma distributed random effects
#> $estimates
#> theta se
#> dependence1 1.9018993 0.7610494
#> dependence2 0.9473264 0.2636491
#>
#> $type
#> [1] "clayton.oakes"
#>
#> $h
#> Estimate Std.Err 2.5% 97.5% P-value
#> dependence1 0.6675 0.1217 0.42901 0.906 4.122e-08
#> dependence2 0.3325 0.1217 0.09399 0.571 6.289e-03
#>
#> $vare
#> NULL
#>
#> $vartot
#> Estimate Std.Err 2.5% 97.5% P-value
#> p1 2.849 0.7315 1.415 4.283 9.823e-05
#>
#> attr(,"class")
#> [1] "summary.mets.twostage"To specify such a model when only the pairs are availble we show how to specify the model. We here use the same marginal “aa” to make the results comparable. The marginal can also be fitted based on available data.
We start by selecting the data related to the pairs, and sets up new id’s and to start we specify the model using the full design with 9 random effects. Below we show how one can use with only the random effects needed for each pair, which is typically simpler.
head(pairs[ssid,])
#> [,1] [,2]
#> [1,] 1 2
#> [2,] 2 3
#> [3,] 3 4
#> [4,] 5 6
#> [5,] 9 11
#> [6,] 10 11
ids <- sort(unique(c(pairs[ssid,])))
pairsids <- c(pairs[ssid,])
pair.new <- matrix(fast.approx(ids,c(pairs[ssid,])),ncol=2)
head(pair.new)
#> [,1] [,2]
#> [1,] 1 2
#> [2,] 2 3
#> [3,] 3 4
#> [4,] 5 6
#> [5,] 7 9
#> [6,] 8 9
dataid <- dsort(data[ids,],"cluster")
outid <- ace.family.design(dataid,member="type",id="cluster")
outid$pardes
#> [,1] [,2]
#> [1,] 0.25 0
#> [2,] 0.25 0
#> [3,] 0.25 0
#> [4,] 0.25 0
#> [5,] 0.25 0
#> [6,] 0.25 0
#> [7,] 0.25 0
#> [8,] 0.25 0
#> [9,] 0.00 1
head(outid$des.rv)
#> m1 m2 m3 m4 f1 f2 f3 f4 env
#> [1,] 1 1 1 1 0 0 0 0 1
#> [2,] 0 0 0 0 1 1 1 1 1
#> [3,] 1 1 0 0 1 1 0 0 1
#> [4,] 1 0 1 0 1 0 1 0 1
#> [5,] 1 1 1 1 0 0 0 0 1
#> [6,] 0 0 0 0 1 1 1 1 1Now fitting the model with the data set up
tsdid <- binomial.twostage(aa,data=dataid,clusters=dataid$cluster,
theta=c(2,1),random.design=outid$des.rv,theta.des=outid$pardes,pairs=pair.new)
summary(tsdid)
#> Dependence parameter for Clayton-Oakes model
#> Variance of Gamma distributed random effects
#> $estimates
#> theta se
#> dependence1 1.9018993 0.7573240
#> dependence2 0.9473264 0.2464157
#>
#> $type
#> [1] "clayton.oakes"
#>
#> $h
#> Estimate Std.Err 2.5% 97.5% P-value
#> dependence1 0.6675 0.1211 0.4302 0.9048 3.514e-08
#> dependence2 0.3325 0.1211 0.0952 0.5698 6.027e-03
#>
#> $vare
#> NULL
#>
#> $vartot
#> Estimate Std.Err 2.5% 97.5% P-value
#> p1 2.849 0.7111 1.455 4.243 6.161e-05
#>
#> attr(,"class")
#> [1] "summary.mets.twostage"We now specify the design specifically using the pairs. The random.design and design on the parameters are now given for each pair, as a 3 dimensional matrix. with a direct specification of random.design and the design on the parameters theta.design. In addition we need also to give the number of random effects for each pair. These basic things are constructed by certain functions for the ACE design.
pair.types <- matrix(dataid[c(t(pair.new)),"type"],byrow=T,ncol=2)
head(pair.new,7)
#> [,1] [,2]
#> [1,] 1 2
#> [2,] 2 3
#> [3,] 3 4
#> [4,] 5 6
#> [5,] 7 9
#> [6,] 8 9
#> [7,] 8 10
head(pair.types,7)
#> [,1] [,2]
#> [1,] "mother" "father"
#> [2,] "father" "child"
#> [3,] "child" "child"
#> [4,] "mother" "father"
#> [5,] "mother" "child"
#> [6,] "father" "child"
#> [7,] "father" "child"
###
theta.des <- rbind( c(rbind(c(1,0), c(1,0), c(0,1), c(0,0))),
c(rbind(c(0.5,0),c(0.5,0),c(0.5,0),c(0,1))))
random.des <- rbind(
c(1,0,1,0),c(0,1,1,0),
c(1,1,0,1),c(1,0,1,1))
mf <- 1*(pair.types[,1]=="mother" & pair.types[,2]=="father")
## pair, rv related to pairs, theta.des related to pair
pairs.new <- cbind(pair.new,(mf==1)*1+(mf==0)*3,(mf==1)*2+(mf==0)*4,(mf==1)*1+(mf==0)*2,(mf==1)*3+(mf==0)*4)pairs.new is matrix with
columns 1:2 giving the indeces of the data points
columns 3:4 giving the indeces of the random.design for the different pairs
columns 5 giving the indeces of the theta.des written as rows
columns 6 giving the number of random variables for this pair
The length of all rows of theta.des are the maximum number of random effects \(\times\) the number of parameters. These two numbers are given in the call. In this case 4 \(\times\) 2. So theta.des has rows of length \(8\), possibly including some 0’s for rows not relevant due to fewer random effects, as is the case here for pairs that do not share genetic effects.
For pair 1 that is a mother/farther pair, we see that they share 1 environmental random effect of size 1. There are also two genetic effects that are unshared between the two. So a total of 3 random effects are needed here. The theta.des relates the 3 random effects to possible relationships in the parameters. Here the genetic effects are full and so is the environmental effect. In contrast we also consider a mother/child pair that share half the genes, now with random effects with (1/2) gene variance. We there need 4 random effects, 2 non-shared half-gene, 1 shared half-gene, and one shared full environmental effect.
# 3 rvs here
random.des
#> [,1] [,2] [,3] [,4]
#> [1,] 1 0 1 0
#> [2,] 0 1 1 0
#> [3,] 1 1 0 1
#> [4,] 1 0 1 1
theta.des
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] 1.0 1.0 0.0 0 0 0 1 0
#> [2,] 0.5 0.5 0.5 0 0 0 0 1
head(pairs.new)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 1 2 1 2 1 3
#> [2,] 2 3 3 4 2 4
#> [3,] 3 4 3 4 2 4
#> [4,] 5 6 1 2 1 3
#> [5,] 7 9 3 4 2 4
#> [6,] 8 9 3 4 2 4Now fitting the model, and we see that it is a lot quicker due to the fewer random effects needed for pairs. We need to also specify the number of parameters in this case.
tsdid2 <- binomial.twostage(aa,data=dataid,clusters=dataid$cluster, theta=c(2,1),
random.design=random.des,theta.des=theta.des,pairs=pairs.new,dim.theta=2)
summary(tsdid2)
#> Dependence parameter for Clayton-Oakes model
#> Variance of Gamma distributed random effects
#> $estimates
#> theta se
#> dependence1 1.9018993 0.7573240
#> dependence2 0.9473264 0.2464157
#>
#> $type
#> [1] "clayton.oakes"
#>
#> $h
#> Estimate Std.Err 2.5% 97.5% P-value
#> dependence1 0.6675 0.1211 0.4302 0.9048 3.514e-08
#> dependence2 0.3325 0.1211 0.0952 0.5698 6.027e-03
#>
#> $vare
#> NULL
#>
#> $vartot
#> Estimate Std.Err 2.5% 97.5% P-value
#> p1 2.849 0.7111 1.455 4.243 6.161e-05
#>
#> attr(,"class")
#> [1] "summary.mets.twostage"The same model can be specifed even simpler via the kinship coefficient. For this speicification there are 4 random effects for each pair, but some have variance 0. The mother-father pair, here shares a random effect with variance 0, and have two non-shared genetic effects with full variance, in addition to a fully shared environmental effect.
kinship <- rep(0.5,nrow(pair.types))
kinship[pair.types[,1]=="mother" & pair.types[,2]=="father"] <- 0
head(kinship,n=10)
#> [1] 0.0 0.5 0.5 0.0 0.5 0.5 0.5 0.5 0.0 0.5
out <- make.pairwise.design(pair.new,kinship,type="ace") tsdid3 <- binomial.twostage(aa,data=dataid,clusters=dataid$cluster,
theta=c(2,1)/9,random.design=out$random.design,
theta.des=out$theta.des,pairs=out$new.pairs,dim.theta=2)
summary(tsdid3)
#> Dependence parameter for Clayton-Oakes model
#> Variance of Gamma distributed random effects
#> $estimates
#> theta se
#> dependence1 1.9018993 0.7573240
#> dependence2 0.9473264 0.2464157
#>
#> $type
#> [1] "clayton.oakes"
#>
#> $h
#> Estimate Std.Err 2.5% 97.5% P-value
#> dependence1 0.6675 0.1211 0.4302 0.9048 3.514e-08
#> dependence2 0.3325 0.1211 0.0952 0.5698 6.027e-03
#>
#> $vare
#> NULL
#>
#> $vartot
#> Estimate Std.Err 2.5% 97.5% P-value
#> p1 2.849 0.7111 1.455 4.243 6.161e-05
#>
#> attr(,"class")
#> [1] "summary.mets.twostage"library(mets)
set.seed(1000)
data <- simbinClaytonOakes.family.ace(10000,2,1,beta=NULL,alpha=NULL)
head(data)
#> ybin x type cluster
#> 1 1 1 mother 1
#> 2 0 1 father 1
#> 3 1 1 child 1
#> 4 1 0 child 1
#> 5 1 1 mother 2
#> 6 1 1 father 2
data$number <- c(1,2,3,4)
data$child <- 1*(data$number==3)
mm <- familycluster.index(data$cluster)
head(mm$familypairindex,n=20)
#> [1] 1 2 1 3 1 4 2 3 2 4 3 4 5 6 5 7 5 8 6 7
pairs <- mm$pairs
dim(pairs)
#> [1] 60000 2
head(pairs,12)
#> [,1] [,2]
#> [1,] 1 2
#> [2,] 1 3
#> [3,] 1 4
#> [4,] 2 3
#> [5,] 2 4
#> [6,] 3 4
#> [7,] 5 6
#> [8,] 5 7
#> [9,] 5 8
#> [10,] 6 7
#> [11,] 6 8
#> [12,] 7 8 dtypes <- interaction( data[pairs[,1],"type"], data[pairs[,2],"type"])
dtypes <- droplevels(dtypes)
table(dtypes)
#> dtypes
#> child.child father.child mother.child mother.father
#> 10000 20000 20000 10000
dm <- model.matrix(~-1+factor(dtypes))Now with the pairs we fit the model
aa <- margbin <- glm(ybin~x,data=data,family=binomial())
tsp <- binomial.twostage(aa,data=data, clusters=data$cluster,
theta.des=dm,pairs=cbind(pairs,1:nrow(dm)))
summary(tsp)
#> Dependence parameter for Odds-Ratio (Plackett) model
#> $estimates
#> theta se
#> factor(dtypes)child.child 4.426619 0.20125594
#> factor(dtypes)father.child 4.439122 0.15953961
#> factor(dtypes)mother.child 4.348930 0.15659723
#> factor(dtypes)mother.father 2.045588 0.08999444
#>
#> $log.or
#> Estimate Std.Err 2.5% 97.5% P-value
#> factor(dtypes)child.child 1.4876 0.04546 1.3985 1.5767 7.978e-235
#> factor(dtypes)father..... 1.4905 0.03594 1.4200 1.5609 0.000e+00
#> factor(dtypes)mother..... 1.4699 0.03601 1.3994 1.5405 0.000e+00
#> factor(dtypes)mother......1 0.7157 0.04399 0.6295 0.8019 1.675e-59
#>
#> $type
#> [1] "plackett"
#>
#> attr(,"class")
#> [1] "summary.mets.twostage"To fit the pairwise odds-ratio model in the case of a pair-specification there are two options for fitting the model.
Starting by the second option. We need to start by specify the design of the odds-ratio of each pair. We set up the data and find all combinations within the pairs. Subsequently, we remove all the empty groups, by grouping together the factor levels 4:9, and then we construct the design.
tdp <-cbind( dataid[pair.new[,1],],dataid[pair.new[,2],])
names(tdp) <- c(paste(names(dataid),"1",sep=""),
paste(names(dataid),"2",sep=""))
tdp <-transform(tdp,tt=interaction(type1,type2))
dlevel(tdp)
#> tt #levels=:6
#> [1] "child.child" "father.child" "mother.child" "child.father"
#> [5] "father.father" "mother.father"
#> -----------------------------------------
drelevel(tdp,newlevels=list(mother.father=4:9)) <- obs.types~tt
dtable(tdp,~tt+obs.types)
#>
#> obs.types mother.father child.child father.child mother.child
#> tt
#> child.child 0 509 0 0
#> father.child 0 0 986 0
#> mother.child 0 0 0 1014
#> child.father 0 0 0 0
#> father.father 0 0 0 0
#> mother.father 491 0 0 0
tdp <- model.matrix(~-1+factor(obs.types),tdp)We then can fit the pairwise model using the pairs and the pair-design for descrbing the OR. The results are consistent with the the ACE model as the mother-father have a lower dependence as is due only the environmental effects. All other combinations should have the same dependence as also seem to be the case.
To fit the OR model it is generally recommended to use the var.link to use the parmetrization with log-odd-ratio regression.
###porpair <- binomial.twostage(aa,data=dataid,clusters=dataid$cluster,
### theta.des=tdp,pairs=pair.new,model="or",var.link=1)
###summary(porpair)sessionInfo()
#> R version 4.2.1 (2022-06-23)
#> Platform: x86_64-apple-darwin21.5.0 (64-bit)
#> Running under: macOS Monterey 12.5.1
#>
#> Matrix products: default
#> BLAS: /usr/local/Cellar/openblas/0.3.21/lib/libopenblasp-r0.3.21.dylib
#> LAPACK: /usr/local/Cellar/r/4.2.1_2/lib/R/lib/libRlapack.dylib
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] mets_1.3.0 lava_1.7.0 timereg_2.0.2 survival_3.3-1
#>
#> loaded via a namespace (and not attached):
#> [1] Rcpp_1.0.9 bslib_0.3.1 compiler_4.2.1
#> [4] jquerylib_0.1.4 tools_4.2.1 digest_0.6.29
#> [7] jsonlite_1.8.0 evaluate_0.15 lattice_0.20-45
#> [10] rlang_1.0.2 Matrix_1.4-1 cli_3.3.0
#> [13] yaml_2.3.5 parallel_4.2.1 mvtnorm_1.1-3
#> [16] xfun_0.30 fastmap_1.1.0 stringr_1.4.0
#> [19] knitr_1.37 sass_0.4.0 globals_0.16.1
#> [22] grid_4.2.1 listenv_0.8.0 R6_2.5.1
#> [25] future.apply_1.9.0 parallelly_1.32.1 rmarkdown_2.13
#> [28] magrittr_2.0.2 codetools_0.2-18 htmltools_0.5.2
#> [31] splines_4.2.1 future_1.27.0 numDeriv_2016.8-1.1
#> [34] stringi_1.7.6