%\VignetteIndexEntry{Analysis of bivariate binomial data: Twin analysis}
%\VignetteEngine{R.rsp::tex}
%\VignetteKeyword{R}
%\VignetteKeyword{package}
%\VignetteKeyword{vignette}
%\VignetteKeyword{LaTeX}
\PassOptionsToPackage{usenames}{xcolor}
\PassOptionsToPackage{hidelinks,colorlinks,linkcolor={blue!50!black},urlcolor={blue!50!black},citecolor={blue!50!black}}{hyperref}
\documentclass[a4paper]{tufte-handout}
\usepackage{etex}
\usepackage{etoolbox}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{mathtools}
\usepackage[strict=true]{csquotes}
\usepackage{xcolor}
\usepackage{natbib}
\usepackage{amsmath,amssymb,textcomp}
\usepackage{bm}
\usepackage{graphicx}
\usepackage{wrapfig}
\usepackage{rotating}
\usepackage{capt-of}
\usepackage{listings}
\usepackage{booktabs}
\usepackage{multicol}
\usepackage{longtable}
\usepackage{zlmtt}
\usepackage{float}
\usepackage{fancyvrb}
\usepackage{verbatim}
\usepackage[english]{babel}
\usepackage{hyperref}
\usepackage{environ}
\NewEnviron{mnote}{\marginnote{\BODY}}
\NewEnviron{snote}{\sidenote{\BODY}}
\usepackage{xspace}
\usepackage{url}
\usepackage{amsthm}
\newtheorem{thm}{Theorem.}
\newtheorem{prop}{Proposition.}
\theoremstyle{definition}
\newtheorem{defn}[thm]{Definition.}
\newtheorem{exa}[thm]{Example}
\newcommand{\R}{\ensuremath{\mathbb{R}}} 
\newcommand{\Real}{\ensuremath{\mathbb{R}}} 
\newcommand{\Complex}{\ensuremath{\mathbb{C}}} 
\newcommand{\Z}{\ensuremath{\mathbb{Z}}} 
\newcommand{\N}{\ensuremath{\mathbb{N}}} 
\newcommand{\norm}[1]{\ensuremath{\left\Vert#1\right\Vert}} 
\newcommand{\var}{\ensuremath{\mathbb{V}\text{ar}}} 
\newcommand{\cov}{\ensuremath{\mathbb{C}\text{ov}}} 
\newcommand{\cor}{\ensuremath{\mathbb{C}\text{or}}} 
\newcommand{\E}{\ensuremath{\mathbb{E}}} 
\newcommand{\pr}{\ensuremath{\mathbb{P}}} 
\newcommand{\from}{\leftarrow} 
\newcommand{\To}[2]{\ensuremath{#1\!\to\!#2}}
\newcommand{\From}[2]{\ensuremath{#1\!\from\!#2}}
\newcommand{\chain}[3]{\ensuremath{#1\!\to\!#2\to\!#3}}
\newcommand{\ichain}[3]{\ensuremath{#1\!\from\!#2\from\!#3}}
\newcommand{\fork}[3]{\ensuremath{#1\!\from\!#2\to\!#3}}
\newcommand{\ifork}[3]{\ensuremath{#1\!\to\!#2\from\!#3}}
\newcommand{\pa}{\text{pa}}
\newcommand{\abs}[1]{\ensuremath{\left\vert#1\right\vert}} 
\newcommand{\ipr}[1]{\langle#1\rangle} 
\newcommand{\set}[1]{\left{#1\right}} 
\newcommand{\seq}[1]{\left<#1\right>} 
\renewcommand{\subset}{\subseteq}  
\renewcommand{\supset}{\supseteq} 
\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
\DeclareMathOperator{\sgn}{sgn}
\DeclareMathOperator{\mvec}{vec}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\im}{im}
\DeclareMathOperator{\diag}{diag}
\DeclareMathOperator{\bias}{bias}
\DeclareMathOperator{\rank}{rank}
\DeclareMathOperator{\logit}{logit}
\DeclareMathOperator{\expit}{expit}
\makeatletter
 \def\mathcolor#1#{\@mathcolor{#1}}
 \def\@mathcolor#1#2#3{%
   \protect\leavevmode
   \begingroup
     \color#1{#2}#3%
  \endgroup
 }
\newcommand{\Dto}{\overset{\mathcal{D}}{\longrightarrow}}
\newcommand{\Pto}{\overset{P}{\longrightarrow}}
\newcommand{\Wto}{\overset{W}{\longrightarrow}}
\newcommand{\VV}{\bm{\Omega}_\theta}
\newcommand{\independenT}[2]{\mathrel{\setbox0\hbox{$#1#2$}\copy0\kern-\wd0\mkern4mu\box0}} 
\newcommand{\indep}{\protect\mathpalette{\protect\independenT}{\perp}}
\DefineVerbatimEnvironment{verbatim}{Verbatim}{xleftmargin=0.5em,xrightmargin=0em,numbers=none,frame=none,fontsize=\footnotesize,formatcom={\color[rgb]{0.2,0.2,0.2}}}
\setlength{\parindent}{0pt} % Kills annoying indents. 
\let\iint\relax 
\let\iiint\relax

\lstset{basicstyle=\ttfamily\small,keywordstyle=\color{black},commentstyle=\color{gray}\ttfamily\itshape,stringstyle=\color[rgb]{0,0,0.5},columns=fullflexible,alsoletter=.,texcl=true,escapeinside={*@}{@*)},escapebegin=\lst@commentstyle\,,breaklines=true,breakatwhitespace=false,numbers=left,numberstyle=\ttfamily\tiny\color{gray},stepnumber=1,numbersep=10pt,backgroundcolor=\color{white},tabsize=4,showspaces=false,showstringspaces=false,xleftmargin=.23in,frame=lines
,rulesepcolor=\color[rgb]{0.85,0.85,0.85},basewidth={0.5em,0.42em},language=r,label=
,caption= ,captionpos=b}
\lstset{basicstyle=\ttfamily\small,keywordstyle=\color{black},commentstyle=\color{gray}\ttfamily\itshape,stringstyle=\color[rgb]{0,0,0.5},columns=fullflexible,alsoletter=.,texcl=true,escapeinside={*@}{@*)},escapebegin=\lst@commentstyle\,,breaklines=true,breakatwhitespace=false,numbers=left,numberstyle=\ttfamily\tiny\color{gray},stepnumber=1,numbersep=10pt,backgroundcolor=\color{white},tabsize=4,showspaces=false,showstringspaces=false,xleftmargin=.23in,frame=lines
,rulesepcolor=\color[rgb]{0.85,0.85,0.85},basewidth={0.5em,0.42em},language=python,label=
,caption= ,captionpos=b}

\setlength{\parindent}{0pt} % Kills annoying indents.
\newcommand{\n}{}
\author{Klaus Holst \& Thomas Scheike}
\date{\today}
\title{Analysis of bivariate binomial data: Twin analysis}
\hypersetup{
 pdfauthor={Klaus Holst \& Thomas Scheike},
 pdftitle={Analysis of bivariate binomial data: Twin analysis},
 pdfkeywords={},
 pdfsubject={},
 pdfcreator={Emacs 26.1 (Org mode 9.1.14)}, 
 pdflang={English}}
\begin{document}

\maketitle
\noindent\rule{\textwidth}{0.5pt}


\section*{Overview}
\label{sec:orgc699b56}

When looking at bivariate binomial data with the aim of learning about the 
dependence that is present, possibly after correcting for some covariates many
models are available. 

\begin{itemize}
\item Random-effects models logistic regression covered elsewhere (glmer in lme4).
\end{itemize}

in the mets package you can fit the 

\begin{itemize}
\item Pairwise odds ratio model
\item Bivariate Probit model 
\begin{itemize}
\item With random effects
\item Special functionality for polygenic random effects modelling 
such as ACE, ADE ,AE and so forth.
\end{itemize}

\item Additive gamma random effects model 
\begin{itemize}
\item Special functionality for polygenic random effects modelling 
such as ACE, ADE ,AE and so forth.
\end{itemize}
\end{itemize}


Typically it can be hard or impossible 
to specify random effects models with special 
structure among the parameters of the random effects. This is possible in
our models. 

To be concrete about the model structure assume that we have paired binomial 
data \(Y_1, Y_2, X_1, X_2\) where the responses are \(Y_1, Y_2\) and we
have covariates \(X_1, X_2\).

We start by giving a brief description of these different models.  First we
for bivariate data one can specify the marginal probability using logistic 
regression models 
\[
logit(P(Y_i=1|X_i)) = \alpha_i + X_i^T \beta  i=1,2.
\]
These model can be estimated under working independence 
\cite{zeger-liang-86}.  


A typical twin analysis will typically consist of  looking at both 

\begin{itemize}
\item Pairwise odds ratio model
\item Bivariate Probit model
\end{itemize}

The additive gamma can be used for the same as the bivariate probit model but 
is more restrictive in terms of dependence structure, but is nevertheless 
still valuable to have also as a check of results of the bivariate probit
model. 


\subsection*{Biprobit with random effects}
\label{sec:org86402a4}

For these model we assume that given random effects \(Z\) and a covariate vector 
\(V_{12}\)
we have independent logistic regression models 
\[
probit(P(Y_i=1|X_i, Z)) = \alpha_i + X_i^T \beta + V_{12}^T Z  i=1,2.
\]
where \(Z\) is a bivariate normal distribution with some covariance 
\(\Sigma\). The general covariance structure 
\(\Sigma\) makes the model very flexible. 

We note that 

\begin{itemize}
\item Paramters \(\beta\) are subject specific
\item The \(\Sigma\) will reflect dependence
\end{itemize}


The more standard link function \(logit\) rather than the \(probit\) link
is often used and implemented in for example 
\cite{mm}.  The advantage is that one now gets an odds-ratio interpretation 
of the subject specific effects, but one then needs  numerical integration to
fit the model. 

\#We note that 


\subsection*{Pairwise odds ratio model}
\label{sec:orgf29d361}

Now the pairwise odds ratio model the specifies that given \(X_1, X_2\)
the marginal models are 
\[
logit(P(Y_i=1|X_i)) = \alpha_i + X_i^T \beta  i=1,2
\]

The primary object of interest are the odds ratio between \(Y_{1}\) and \(Y_{2}\)
\[
\gamma_{12} = \frac{ P(  Y_{ki} =1 , Y_{kj} =1) P(  Y_{ki} =0 , Y_{kj} =0) }{ 
  P(  Y_{ki} =1 , Y_{kj} =0) P(  Y_{ki} =0 , Y_{kj} =1) }
\]
given \(X_{ki}\), \(X_{kj}\), and \(Z_{kji}\). 

We model the odds ratio with the regression 
\[
\gamma_{12} = \exp( Z_{12}^T \lambda)
\]
Where \(Z_{12}\) are some covarites that may influence the odds-ratio 
between between \(Y_{1}\) and \(Y_{2}\) and contains the marginal covariates,
\cite{carey-1993,dale1986global,palmgren1989,molenberghs1994marginal}. 
This odds-ratio is given covariates as well as marginal covariates. 
The odds-ratio and marginals specify the joint bivariate distribution via
the so-called Placckett-distribution. 

One way of fitting this model is the ALR algoritm, the alternating 
logistic regression ahd this has been described in several papers
\cite{kuk2004permutation,kuk2007hybrid,qaqish2012orthogonalized}.
We here simply estimate the parameters in a two stage-procedure

\begin{itemize}
\item Estimating the marginal parameters via GEE
\item Using marginal estimates, estimate dependence parameters
\end{itemize}

This gives efficient estimates of the dependence parameters because of
orthogonality, but some efficiency may be gained for the marginal parameters 
by using the full likelihood or iterative fitting such as for the ALR. 


The pairwise odds-ratio model is very useful, but one do not have a random 
effects model. 


\subsection*{Additive gamma model}
\label{sec:orgdb54e35}

Again we operate under  marginal logistic regression models are 
\[
logit(P(Y_i=1|X_i)) = \alpha_i + X_i^T \beta  i=1,2
\]

First with just one random effect \(Z\) we assume that  conditional
on \(Z\) the responses are independent  and follow the model 
\[
logit(P(Y_i=1|X_i,Z)) = exp( -Z \cdot \Psi^{-1}(\lambda_{\bullet},\lambda_{\bullet},P(Y_i=1|X_i)) )  
\]
where \(\Psi\) is the laplace transform of \(Z\) where we assume that
\(Z\) is gamma distributed with variance \(\lambda_{\bullet}^{-1}\) and mean 1. 
In general \(\Psi(\lambda_1,\lambda_2)\) is the laplace transform of  a Gamma distributed random 
effect with \(Z\) with mean \(\lambda_1/\lambda_2\) and variance \(\lambda_1/\lambda_2^2\).

We fit this model by 

\begin{itemize}
\item Estimating the marginal parameters via GEE
\item Using marginal estimates, estimate dependence parameters
\end{itemize}

To deal with multiple random effects we consider random effects 
\(Z_i  i=1,...,d\)   such that  \(Z_i\) is gamma distributed with 
mean \(\lambda_j/\lambda_{\bullet}\) and variance \(\lambda_j/\lambda_{\bullet}^2\), where we define the scalar \(\lambda_{\bullet}\) below. 

Now given a cluster-specific design vector \(V_{12}\) we assume that 
\[
V_{12}^T Z
\]
is gamma distributed with mean 1 and variance \(\lambda_{\bullet}^{-1}\) 
such that critically the random effect variance is the same for all clusters.
That is 
\[
 \lambda_{\bullet} = V_{12}^T (\lambda_1,...,\lambda_d)^T 
\]
We return to some specific models below, and show how to fit the ACE and AE 
model using this set-up. 

One last option in the model-specification is to specify how the 
parameters \(\lambda_1,...,\lambda_d\) are related. We thus can specify a 
matrix \(M\) of dimension \(p \times d\) such that 
\[
 (\lambda_1,...,\lambda_d)^T  = M \theta
\]
where \(\theta\) is d-dimensional.  If \(M\) is diagonal we have no 
restrictions on parameters. 

This parametrization is obtained with the var.par=0 option that thus estimates
\(\theta\).

The DEFAULT parametrization instead estimates the variances of the random effecs (var.par=1)
via the parameters \(\nu\) 
\[
 M \nu = ( \lambda_1/\lambda_{\bullet}^2, ...,\lambda_d/\lambda_{\bullet}^2)^T
\]


The basic modelling assumption is now that given random effects 
\(Z=(Z_1,...,Z_d)\) we have independent probabilites 
\[
logit(P(Y_i=1|X_i,Z)) = exp( -V_{12,i}^T Z \cdot \Psi^{-1}(\lambda_{\bullet},\lambda_{\bullet},P(Y_i=1|X_i)) )   i=1,2
\]

We fit this model by 

\begin{itemize}
\item Estimating the marginal parameters via GEE
\item Using marginal estimates, estimate dependence parameters
\end{itemize}

Even though the model not formaly in this formulation allows negative 
correlation in practice the paramters can be negative and this reflects
negative correlation. An advanatage is that no numerical integration is 
needed. 


\section*{The twin-stutter data}
\label{sec:orgf3a7afd}

We consider the twin-stutter where for pairs of twins that are 
either dizygotic or monozygotic we have recorded whether the twins
are stuttering \cite{twinstut-ref}

We here consider MZ and same sex DZ twins. 

Looking at the data 

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
library(mets)
data(twinstut)
twinstut$binstut <- 1*(twinstut$stutter=="yes")
twinsall <- twinstut
twinstut <- subset(twinstut,zyg%in%c("mz","dz"))
head(twinstut)
\end{lstlisting}

\begin{verbatim}
Loading required package: timereg
Loading required package: survival
Loading required package: lava
lava version 1.5.1
mets version 1.2.1.2

Attaching package: ‘mets’

The following object is masked _by_ ‘.GlobalEnv’:

    object.defined

Warning message:
failed to assign RegisteredNativeSymbol for cor to cor since cor is already defined in the ‘mets’ namespace
   tvparnr zyg stutter    sex age nr binstut
1  2001005  mz      no female  71  1       0
2  2001005  mz      no female  71  2       0
3  2001006  dz      no female  71  1       0
8  2001012  mz      no female  71  1       0
9  2001012  mz      no female  71  2       0
11 2001015  dz      no   male  71  1       0
\end{verbatim}


\section*{Pairwise odds ratio model}
\label{sec:orgbf74e81}

We start by fitting an overall dependence OR for both MZ and DZ even though 
the dependence is expected to be different across zygosity.

The first step is to fit the marginal model adjusting for marginal covariates. 
We here note that there is a rather strong gender effect in the risk of
stuttering. 

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
margbin <- glm(binstut~factor(sex)+age,data=twinstut,family=binomial())
summary(margbin)
\end{lstlisting}

\begin{verbatim}

Call:
glm(formula = binstut ~ factor(sex) + age, family = binomial(), 
    data = twinstut)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.4419  -0.4078  -0.2842  -0.2672   2.6395  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -3.027625   0.104012 -29.108  < 2e-16 ***
factor(sex)male  0.869826   0.062197  13.985  < 2e-16 ***
age             -0.005983   0.002172  -2.754  0.00588 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9328.6  on 21287  degrees of freedom
Residual deviance: 9117.0  on 21285  degrees of freedom
AIC: 9123

Number of Fisher Scoring iterations: 6
\end{verbatim}

Now estimating the OR parameter. We see a strong dependence with an OR
at around 8 that is clearly significant. 

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
bina <- binomial.twostage(margbin,data=twinstut,var.link=1,
		       clusters=twinstut$tvparnr,detail=0)
summary(bina)
\end{lstlisting}

\begin{verbatim}
Dependence parameter for Odds-Ratio (Plackett) model 
With log-link 
$estimates
               theta        se
dependence1 2.085347 0.1274536

$or
            Estimate Std.Err 2.5% 97.5% P-value
dependence1     8.05    1.03 6.04  10.1 4.3e-15

$type
[1] "plackett"

attr(,"class")
[1] "summary.mets.twostage"
\end{verbatim}

Now, and more interestingly, we consider an OR that depends on zygosity and
note that MZ have a much larger OR than DZ twins. This type of trait is 
somewhat complicated to interpret, but clearly, one option is that 
that there is a genetic effect, alternatively there might be
a stronger environmental effect for MZ twins. 


\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
# design for OR dependence 
theta.des <- model.matrix( ~-1+factor(zyg),data=twinstut)
bin <- binomial.twostage(margbin,data=twinstut,var.link=1,
			  clusters=twinstut$tvparnr,theta.des=theta.des)
summary(bin)
\end{lstlisting}

\begin{verbatim}
Dependence parameter for Odds-Ratio (Plackett) model 
With log-link 
$estimates
                  theta        se
factor(zyg)dz 0.5221651 0.2401355
factor(zyg)mz 3.4853933 0.1866076

$or
              Estimate Std.Err   2.5% 97.5%  P-value
factor(zyg)dz     1.69   0.405  0.892  2.48 3.12e-05
factor(zyg)mz    32.64   6.090 20.699 44.57 8.38e-08

$type
[1] "plackett"

attr(,"class")
[1] "summary.mets.twostage"
\end{verbatim}


We now consider further regression modelling of the OR structure by
considering possible interactions between sex and zygozsity.
We see that MZ has a much higher dependence and that males have
a much lower dependence. We tested for interaction in this model and 
these were not significant. 

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
twinstut$cage <- scale(twinstut$age)
theta.des <- model.matrix( ~-1+factor(zyg)+factor(sex),data=twinstut)
bina <- binomial.twostage(margbin,data=twinstut,var.link=1,
			  clusters=twinstut$tvparnr,theta.des=theta.des)
summary(bina)
\end{lstlisting}

\begin{verbatim}
Dependence parameter for Odds-Ratio (Plackett) model 
With log-link 
$estimates
                     theta        se
factor(zyg)dz    0.8098841 0.3138423
factor(zyg)mz    3.7318076 0.2632250
factor(sex)male -0.4075409 0.3055349

$or
                Estimate Std.Err   2.5% 97.5%  P-value
factor(zyg)dz      2.248   0.705  0.865  3.63 0.001441
factor(zyg)mz     41.755  10.991 20.213 63.30 0.000145
factor(sex)male    0.665   0.203  0.267  1.06 0.001064

$type
[1] "plackett"

attr(,"class")
[1] "summary.mets.twostage"
\end{verbatim}


\subsection*{Alternative syntax}
\label{sec:org3f01206}

We now demonstrate how the models can fitted jointly and with anohter
syntax, that ofcourse just fits the marginal model and subsequently fits
the pairwise OR model. 

First noticing as before that MZ twins have a much higher dependence. 

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
 # refers to zygosity of first subject in eash pair : zyg1
 # could also use zyg2 (since zyg2=zyg1 within twinpair's)
 out <- easy.binomial.twostage(stutter~factor(sex)+age,data=twinstut,
		response="binstut",id="tvparnr",var.link=1,
		theta.formula=~-1+factor(zyg1))
summary(out)
\end{lstlisting}

\begin{verbatim}
Dependence parameter for Odds-Ratio (Plackett) model 
With log-link 
$estimates
                   theta        se
factor(zyg1)dz 0.5221651 0.2401355
factor(zyg1)mz 3.4853933 0.1866076

$or
               Estimate Std.Err   2.5% 97.5%  P-value
factor(zyg1)dz     1.69   0.405  0.892  2.48 3.12e-05
factor(zyg1)mz    32.64   6.090 20.699 44.57 8.38e-08

$type
[1] "plackett"

attr(,"class")
[1] "summary.mets.twostage"
\end{verbatim}

Now considering all data and estimating separate effects for the OR
for opposite sex DZ twins and same sex twins. 
We here find that os twins are not markedly different from the same sex DZ 
twins. 

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
# refers to zygosity of first subject in eash pair : zyg1
# could also use zyg2 (since zyg2=zyg1 within twinpair's))

desfs<-function(x,num1="zyg1",num2="zyg2")
	c(x[num1]=="dz",x[num1]=="mz",x[num1]=="os")*1

margbinall <- glm(binstut~factor(sex)+age,data=twinsall,family=binomial())
out3 <- easy.binomial.twostage(binstut~factor(sex)+age,
      data=twinsall,response="binstut",id="tvparnr",var.link=1,
      theta.formula=desfs,desnames=c("dz","mz","os"))
summary(out3)
\end{lstlisting}


\begin{verbatim}
Dependence parameter for Odds-Ratio (Plackett) model 
With log-link 
$estimates
       theta        se
dz 0.5278527 0.2396796
mz 3.4850037 0.1864190
os 0.7802940 0.2894394

$or
   Estimate Std.Err   2.5% 97.5%  P-value
dz     1.70   0.406  0.899  2.49 3.02e-05
mz    32.62   6.081 20.703 44.54 8.13e-08
os     2.18   0.632  0.944  3.42 5.50e-04

$type
[1] "plackett"

attr(,"class")
[1] "summary.mets.twostage"
\end{verbatim}



\section*{Bivariate Probit model}
\label{sec:org1b646d8}

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
library(mets)
data(twinstut)
twinstut <- subset(twinstut,zyg%in%c("mz","dz"))
twinstut$binstut <- 1*(twinstut$stutter=="yes")
head(twinstut)
\end{lstlisting}

\begin{verbatim}
   tvparnr zyg stutter    sex age nr binstut
1  2001005  mz      no female  71  1       0
2  2001005  mz      no female  71  2       0
3  2001006  dz      no female  71  1       0
8  2001012  mz      no female  71  1       0
9  2001012  mz      no female  71  2       0
11 2001015  dz      no   male  71  1       0
\end{verbatim}

First testing for same dependence in MZ and DZ  that we recommend doing by 
comparing the correlations of MZ and DZ twins. Apart from regression 
correction in the mean this is an un-structured model, and the useful
concordance and casewise concordance estimates can be reported from this
analysis. 

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
b1 <- bptwin(binstut~sex,data=twinstut,id="tvparnr",zyg="zyg",DZ="dz",type="un")
summary(b1)
\end{lstlisting}

\begin{verbatim}

                Estimate    Std.Err          Z p-value
(Intercept)    -1.794823   0.023289 -77.066728  0.0000
sexmale         0.401432   0.030179  13.301813  0.0000
atanh(rho) MZ   1.096916   0.073574  14.909087  0.0000
atanh(rho) DZ   0.132458   0.062516   2.118800  0.0341

 Total MZ/DZ Complete pairs MZ/DZ
 8777/12511  3255/4058           

                           Estimate 2.5%    97.5%  
Tetrachoric correlation MZ 0.79939  0.74101 0.84577
Tetrachoric correlation DZ 0.13169  0.00993 0.24960

MZ:
                     Estimate 2.5%     97.5%   
Concordance           0.01698  0.01411  0.02042
Casewise Concordance  0.46730  0.40383  0.53185
Marginal              0.03634  0.03287  0.04016
Rel.Recur.Risk       12.85882 10.87510 14.84253
log(OR)               3.75632  3.37975  4.13289
DZ:
                     Estimate 2.5%    97.5%  
Concordance          0.00235  0.00140 0.00393
Casewise Concordance 0.06456  0.03937 0.10413
Marginal             0.03634  0.03287 0.04016
Rel.Recur.Risk       1.77662  0.92746 2.62577
log(OR)              0.63527  0.09013 1.18040

                         Estimate 2.5% 97.5%
Broad-sense heritability   1      NaN  NaN
\end{verbatim}


\subsection*{Polygenic modelling}
\label{sec:org8788a7d}

We now turn attention to specific polygenic modelling where special random 
effects are used to specify ACE, AE, ADE models and so forth. This is very
easy with the bptwin function. The key parts of the output are the sizes of 
the genetic component A and the environmental component, and we can compare 
with the results of the unstructed model above. Also formally we can test 
if this submodel is acceptable by a likelihood ratio test. 


\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
b1 <- bptwin(binstut~sex,data=twinstut,id="tvparnr",zyg="zyg",DZ="dz",type="ace")
summary(b1)
\end{lstlisting}

\begin{verbatim}

             Estimate   Std.Err         Z p-value
(Intercept)  -3.70371   0.24449 -15.14855       0
sexmale       0.83310   0.08255  10.09201       0
log(var(A))   1.18278   0.17179   6.88512       0
log(var(C)) -29.99519        NA        NA      NA

 Total MZ/DZ Complete pairs MZ/DZ
 8777/12511  3255/4058           

                   Estimate 2.5%    97.5%  
A                  0.76545  0.70500 0.82590
C                  0.00000  0.00000 0.00000
E                  0.23455  0.17410 0.29500
MZ Tetrachoric Cor 0.76545  0.69793 0.81948
DZ Tetrachoric Cor 0.38272  0.35210 0.41253

MZ:
                     Estimate 2.5%     97.5%   
Concordance           0.01560  0.01273  0.01912
Casewise Concordance  0.42830  0.36248  0.49677
Marginal              0.03643  0.03294  0.04027
Rel.Recur.Risk       11.75741  9.77237 13.74246
log(OR)               3.52382  3.13466  3.91298
DZ:
                     Estimate 2.5%    97.5%  
Concordance          0.00558  0.00465 0.00670
Casewise Concordance 0.15327  0.13749 0.17050
Marginal             0.03643  0.03294 0.04027
Rel.Recur.Risk       4.20744  3.78588 4.62900
log(OR)              1.69996  1.57262 1.82730

                         Estimate 2.5%    97.5%  
Broad-sense heritability 0.76545  0.70500 0.82590
\end{verbatim}


\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}

b0 <- bptwin(binstut~sex,data=twinstut,id="tvparnr",zyg="zyg",DZ="dz",type="ae")
summary(b0)
\end{lstlisting}

\begin{verbatim}

             Estimate   Std.Err         Z p-value
(Intercept)  -3.70371   0.24449 -15.14855       0
sexmale       0.83310   0.08255  10.09201       0
log(var(A))   1.18278   0.17179   6.88512       0

 Total MZ/DZ Complete pairs MZ/DZ
 8777/12511  3255/4058           

                   Estimate 2.5%    97.5%  
A                  0.76545  0.70500 0.82590
E                  0.23455  0.17410 0.29500
MZ Tetrachoric Cor 0.76545  0.69793 0.81948
DZ Tetrachoric Cor 0.38272  0.35210 0.41253

MZ:
                     Estimate 2.5%     97.5%   
Concordance           0.01560  0.01273  0.01912
Casewise Concordance  0.42830  0.36248  0.49677
Marginal              0.03643  0.03294  0.04027
Rel.Recur.Risk       11.75741  9.77237 13.74246
log(OR)               3.52382  3.13466  3.91298
DZ:
                     Estimate 2.5%    97.5%  
Concordance          0.00558  0.00465 0.00670
Casewise Concordance 0.15327  0.13749 0.17050
Marginal             0.03643  0.03294 0.04027
Rel.Recur.Risk       4.20744  3.78588 4.62900
log(OR)              1.69996  1.57262 1.82730

                         Estimate 2.5%    97.5%  
Broad-sense heritability 0.76545  0.70500 0.82590
\end{verbatim}



\section*{Additive gamma random effects}
\label{sec:orgff761d3}

Fitting first a model with different size random effects for MZ and DZ. We 
note that as before in the OR and biprobit model the dependence is much
stronger for MZ twins. We also test if these are the same by parametrizing the
OR model with an intercept. This clearly shows a significant difference. 


\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
theta.des <- model.matrix( ~-1+factor(zyg),data=twinstut)
margbin <- glm(binstut~sex,data=twinstut,family=binomial())
bintwin <- binomial.twostage(margbin,data=twinstut,model="gamma",
     clusters=twinstut$tvparnr,detail=0,theta=c(0.1)/1,var.link=1,
     theta.des=theta.des)
summary(bintwin)

# test for same dependence in MZ and DZ 
theta.des <- model.matrix( ~factor(zyg),data=twinstut)
margbin <- glm(binstut~sex,data=twinstut,family=binomial())
bintwin <- binomial.twostage(margbin,data=twinstut,model="gamma",
     clusters=twinstut$tvparnr,detail=0,theta=c(0.1)/1,var.link=1,
     theta.des=theta.des)
summary(bintwin)
\end{lstlisting}

\begin{verbatim}
Dependence parameter for Clayton-Oakes model
Variance of Gamma distributed random effects 
With log-link 
$estimates
                    theta        se
factor(zyg)dz -2.61194495 0.4854454
factor(zyg)mz -0.01817181 0.1030735

$vargam
              Estimate Std.Err    2.5% 97.5%  P-value
factor(zyg)dz   0.0734  0.0356 0.00356 0.143 3.94e-02
factor(zyg)mz   0.9820  0.1012 0.78361 1.180 2.96e-22

$type
[1] "gamma"

attr(,"class")
[1] "summary.mets.twostage"
Dependence parameter for Clayton-Oakes model
Variance of Gamma distributed random effects 
With log-link 
$estimates
                  theta        se
(Intercept)   -2.611945 0.4854454
factor(zyg)mz  2.593773 0.4962675

$vargam
              Estimate Std.Err    2.5%  97.5% P-value
(Intercept)     0.0734  0.0356 0.00356  0.143  0.0394
factor(zyg)mz  13.3802  6.6401 0.36573 26.395  0.0439

$type
[1] "gamma"

attr(,"class")
[1] "summary.mets.twostage"
\end{verbatim}


\subsection*{Polygenic modelling}
\label{sec:orgc97abc0}

   First setting up the random effects design for the random effects and 
   the the relationship between variance parameters.
   We see that the genetic random effect has size one for MZ and 0.5 for DZ subjects, 
   that have shared and non-shared genetic components with variance 0.5 such that the total 
   genetic variance is the same for all subjects. The shared environmental effect is the samme for 
all. Thus two parameters with these bands. 

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
out <- twin.polygen.design(twinstut,id="tvparnr",zygname="zyg",zyg="dz",type="ace")
head(cbind(out$des.rv,twinstut$tvparnr),10)
out$pardes
\end{lstlisting}

\begin{verbatim}
   MZ DZ DZns1 DZns2 env        
1   1  0     0     0   1 2001005
2   1  0     0     0   1 2001005
3   0  1     1     0   1 2001006
8   1  0     0     0   1 2001012
9   1  0     0     0   1 2001012
11  0  1     1     0   1 2001015
12  0  1     1     0   1 2001016
13  0  1     0     1   1 2001016
15  0  1     1     0   1 2001020
18  0  1     1     0   1 2001022
     [,1] [,2]
[1,]  1.0    0
[2,]  0.5    0
[3,]  0.5    0
[4,]  0.5    0
[5,]  0.0    1
\end{verbatim}


Now, fitting the ACE model, we see that the variance of the genetic, 
component, is 1.5 and the environmental variance is -0.5. Thus suggesting that 
the ACE model does not fit the data.  When the random design is given we 
automatically use the gamma fralty model. 

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}

margbin <- glm(binstut~sex,data=twinstut,family=binomial())
bintwin1 <- binomial.twostage(margbin,data=twinstut,
     clusters=twinstut$tvparnr,detail=0,theta=c(0.1)/1,var.link=0,
     random.design=out$des.rv,theta.des=out$pardes)
summary(bintwin1)
\end{lstlisting}

\begin{verbatim}
Dependence parameter for Clayton-Oakes model
Variance of Gamma distributed random effects 
$estimates
                 theta        se
dependence1  1.5261839 0.2475041
dependence2 -0.5447955 0.1942159

$type
[1] "clayton.oakes"

$h
            Estimate Std.Err   2.5%  97.5%  P-value
dependence1    1.555   0.187  1.189  1.922 9.11e-17
dependence2   -0.555   0.187 -0.922 -0.189 2.99e-03

$vare
NULL

$vartot
   Estimate Std.Err  2.5% 97.5%  P-value
p1    0.981   0.102 0.781  1.18 8.29e-22

attr(,"class")
[1] "summary.mets.twostage"
\end{verbatim}


For this model we estimate the concordance and casewise concordance as well 
as the marginal rates of stuttering for females. 

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
concordanceTwinACE(bintwin1,type="ace")
\end{lstlisting}

\begin{verbatim}
$MZ
                     Estimate Std.Err   2.5%  97.5%  P-value
concordance            0.0182 0.00147 0.0153 0.0211 2.61e-35
casewise concordance   0.5033 0.03256 0.4395 0.5672 6.49e-54
marginal               0.0362 0.00188 0.0325 0.0399 7.15e-83

$DZ
                     Estimate  Std.Err   2.5%   97.5%  P-value
concordance           0.00235 0.000589 0.0012 0.00351 6.45e-05
casewise concordance  0.06501 0.015836 0.0340 0.09604 4.04e-05
marginal              0.03620 0.001877 0.0325 0.03988 7.15e-83
\end{verbatim}


The E component was not consistent with the fit of the data and we
now consider instead the AE model. 


\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
out <- twin.polygen.design(twinstut,id="tvparnr",zygname="zyg",zyg="dz",type="ae")

bintwin <- binomial.twostage(margbin,data=twinstut,
     clusters=twinstut$tvparnr,detail=0,theta=c(0.1)/1,var.link=0,
     random.design=out$des.rv,theta.des=out$pardes)
summary(bintwin)
\end{lstlisting}

\begin{verbatim}
Dependence parameter for Clayton-Oakes model
Variance of Gamma distributed random effects 
$estimates
                theta         se
dependence1 0.9094847 0.09536268

$type
[1] "clayton.oakes"

$h
            Estimate Std.Err 2.5% 97.5% P-value
dependence1        1       0    1     1       0

$vare
NULL

$vartot
   Estimate Std.Err  2.5% 97.5%  P-value
p1    0.909  0.0954 0.723   1.1 1.47e-21

attr(,"class")
[1] "summary.mets.twostage"
\end{verbatim}

Again, the concordance can be computed: 

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
concordanceTwinACE(bintwin,type="ae")
\end{lstlisting}

\begin{verbatim}
$MZ
                     Estimate Std.Err   2.5%  97.5%  P-value
concordance            0.0174 0.00143 0.0146 0.0202 5.00e-34
casewise concordance   0.4795 0.03272 0.4154 0.5437 1.20e-48
marginal               0.0362 0.00188 0.0325 0.0399 7.15e-83

$DZ
                     Estimate  Std.Err   2.5%   97.5%   P-value
concordance           0.00477 0.000393 0.0040 0.00554  5.94e-34
casewise concordance  0.13175 0.005417 0.1211 0.14237 1.14e-130
marginal              0.03620 0.001877 0.0325 0.03988  7.15e-83
\end{verbatim}
\end{document}