% Created 2018-11-19 Mon 19:21
%\VignetteIndexEntry{Twin analysis of quantitative outcomes}
%\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}{}
\usepackage{units}
\author{Klaus Holst \& Thomas Scheike}
\date{\today}
\title{Twin analysis}
\hypersetup{
 pdfauthor={Klaus Holst \& Thomas Scheike},
 pdftitle={Twin analysis},
 pdfkeywords={},
 pdfsubject={},
 pdfcreator={Emacs 26.1 (Org mode 9.1.14)}, 
 pdflang={English}}
\begin{document}

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

\section*{Mets package}
\label{sec:orgeaaf733}

This document provides a brief tutorial to analyzing twin data using
the \textbf{\texttt{mets}} package:
\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
library("mets")
options(warn=-1)
\end{lstlisting}

The development version may be installed from \emph{github}, i.e., with the
\texttt{devtools} package:
\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
devtools::install_github("kkholst/lava")
devtools::install_github("kkholst/mets")
\end{lstlisting}

\section*{Twin analysis, continuous traits}
\label{sec:org0778a55}

  In the following we examine the heritability of Body Mass Index\n{}\cite{korkeila_bmi_1991} \cite{hjelmborg_bmi_2008}, based on
data on self-reported BMI-values from a random sample of 11,411 same-sex twins. First, we will load data

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
data("twinbmi")
head(twinbmi)
\end{lstlisting}

\begin{verbatim}
  tvparnr      bmi      age gender zyg id num
1       1 26.33289 57.51212   male  DZ  1   1
2       1 25.46939 57.51212   male  DZ  1   2
3       2 28.65014 56.62696   male  MZ  2   1
5       3 28.40909 57.73097   male  DZ  3   1
7       4 27.25089 53.68683   male  DZ  4   1
8       4 28.07504 53.68683   male  DZ  4   2
\end{verbatim}

The data is on \emph{long} format with one subject per row.
\begin{mnote}
\begin{description}
\item[{\textbf{\texttt{tvparnr}}}] twin id
\item[{\textbf{\texttt{bmi}}}] Body Mass Index (\(\unitfrac{kg}{m^2}\))
\item[{\textbf{\texttt{age}}}] Age (years)
\item[{\textbf{\texttt{gender}}}] Gender factor (male,female)
\item[{\textbf{\texttt{zyg}}}] zygosity (MZ,DZ)
\end{description}
\end{mnote}


we transpose the data allowing us to do pairwise analyses
\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
twinwide <- fast.reshape(twinbmi, id="tvparnr",varying=c("bmi"))
head(twinwide)
\end{lstlisting}

\begin{verbatim}
   tvparnr     bmi1      age gender zyg id num     bmi2
1        1 26.33289 57.51212   male  DZ  1   1 25.46939
3        2 28.65014 56.62696   male  MZ  2   1       NA
5        3 28.40909 57.73097   male  DZ  3   1       NA
7        4 27.25089 53.68683   male  DZ  4   1 28.07504
9        5 27.77778 52.55838   male  DZ  5   1       NA
11       6 28.04282 52.52231   male  DZ  6   1 22.30936
\end{verbatim}

Next we plot the association within each zygosity group

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

scatterdens <- function(x) {    
    sp <- ggplot(x,  
		aes_string(colnames(x)[1], colnames(x)[2])) +
	theme_minimal() +
	geom_point(alpha=0.3) + geom_density_2d()
    xdens <- ggplot(x, aes_string(colnames(x)[1],fill=1)) + 
	theme_minimal() +
	geom_density(alpha=.5)+ 
	theme(axis.text.x = element_blank(),
	      legend.position = "none") + labs(x=NULL)
    ydens <- ggplot(x, aes_string(colnames(x)[2],fill=1)) +
	theme_minimal() +
	geom_density(alpha=.5) +
	theme(axis.text.y = element_blank(),
	      axis.text.x = element_text(angle=90, vjust=0),
	      legend.position = "none") +
	labs(x=NULL) +
	coord_flip()
    g <- plot_grid(xdens,NULL,sp,ydens,
		  ncol=2,nrow=2,
		  rel_widths=c(4,1.4),rel_heights=c(1.4,4))    
    return(g)
}
\end{lstlisting}

We here show the log-transformed data which is slightly more symmetric
and more appropiate for the twin analysis (see Figure \ref{fig:scatter1} and \ref{fig:scatter2})

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label=scatter1,caption= ,captionpos=b}
\begin{lstlisting}
mz <- log(subset(twinwide, zyg=="MZ")[,c("bmi1","bmi2")])
scatterdens(mz)
\end{lstlisting}

\begin{marginfigure}
\begin{center}
\includegraphics[width=\textwidth]{scatter1.jpg}
\end{center}

\captionof{figure}{Scatter plot of logarithmic BMI measurements in MZ twins.}
\label{fig:scatter1}
\end{marginfigure}

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label=scatter2,caption= ,captionpos=b}
\begin{lstlisting}
dz <- log(subset(twinwide, zyg=="DZ")[,c("bmi1","bmi2")])
scatterdens(dz)
\end{lstlisting}

\begin{marginfigure}
\begin{center}
\includegraphics[width=\textwidth]{scatter2.jpg}
\end{center}

\captionof{figure}{Scatter plot of logarithmic BMI measurements in DZ twins.}
\label{fig:scatter2}
\end{marginfigure}

The plots and raw association measures shows considerable stronger
dependence in the MZ twins, thus indicating genetic influence of the
trait
\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
cor.test(mz[,1],mz[,2], method="spearman")
\end{lstlisting}

\begin{verbatim}

	Spearman's rank correlation rho

data:  mz[, 1] and mz[, 2]
S = 165460000, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.6956209
\end{verbatim}

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
cor.test(dz[,1],dz[,2], method="spearman")
\end{lstlisting}

\begin{verbatim}

	Spearman's rank correlation rho

data:  dz[, 1] and dz[, 2]
S = 2162500000, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.4012686
\end{verbatim}

Ńext we examine the marginal distribution (GEE model with working
independence)

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
l0 <- lm(bmi ~ gender + I(age-40), data=twinbmi)
estimate(l0, id=twinbmi$tvparnr)
\end{lstlisting}

\begin{verbatim}
            Estimate  Std.Err    2.5%   97.5%    P-value
(Intercept)  23.3687 0.054534 23.2618 23.4756  0.000e+00
gendermale    1.4077 0.073216  1.2642  1.5512  2.230e-82
I(age - 40)   0.1177 0.004787  0.1083  0.1271 1.499e-133
\end{verbatim}

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
library("splines")
l1 <- lm(bmi ~ gender*ns(age,3), data=twinbmi)
marg1 <- estimate(l1, id=twinbmi$tvparnr)
\end{lstlisting}

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label=marg1,caption= ,captionpos=b}
\begin{lstlisting}
dm <- Expand(twinbmi,
	    bmi=0,
	    gender=c("male"),
	    age=seq(33,61,length.out=50))
df <- Expand(twinbmi,
	    bmi=0,
	    gender=c("female"),
	    age=seq(33,61,length.out=50))

plot(marg1, function(p) model.matrix(l1,data=dm)%*%p, 
     data=dm["age"], ylab="BMI", xlab="Age",
     ylim=c(22,26.5))
plot(marg1, function(p) model.matrix(l1,data=df)%*%p,  
     data=df["age"], col="red", add=TRUE)
legend("bottomright", c("Male","Female"), 
       col=c("black","red"), lty=1, bty="n")
\end{lstlisting}


\begin{marginfigure}
\begin{center}
\includegraphics[width=\textwidth]{marg1.jpg}
\end{center}


\captionof{figure}{...}
\label{fig:marg1}
\end{marginfigure}

\subsection*{Polygenic model}
\label{sec:org336bb16}

Decompose outcome into
\begin{align*}
Y_i = A_i + D_i + C + E_i, \quad i=1,2
 \end{align*}

\begin{description}
\item[{\(A\)}] Additive genetic effects of alleles
\item[{\(D\)}] Dominante genetic effects of alleles
\item[{\(C\)}] Shared environmental effects
\item[{\(E\)}] Unique environmental genetic effects
\end{description}

Dissimilarity of MZ twins arises from unshared environmental effects
only! \(\cor(E_1,E_2)=0\) and
\begin{align*}
\cor(A_1^{MZ},A_2^{MZ}) = 1, \quad 
\cor(D_1^{MZ},D_2^{MZ}) = 1,
\end{align*}
\begin{align*}
\cor(A_1^{DZ},A_2^{DZ}) = 0.5, \quad
\cor(D_1^{DZ},D_2^{DZ}) = 0.25,
\end{align*}

\begin{align*}
Y_i = A_i + C_i + D_i + E_i
\end{align*}
\begin{align*}
A_i \sim\mathcal{N}(0,\sigma_A^2), C_i
\sim\mathcal{N}(0,\sigma_C^2), D_i
\sim\mathcal{N}(0,\sigma_D^2),
E_i \sim\mathcal{N}(0,\sigma_E^2)
\end{align*}

  \begin{gather*}
    \cov(Y_{1},Y_{2}) = \\
    \begin{pmatrix}
      \sigma_A^2 & 2\Phi\sigma_A^2 \\
      2\Phi\sigma_A^2 & \sigma_A^2
    \end{pmatrix} + 
    \begin{pmatrix}
      \sigma_C^2 & \sigma_C^2 \\
      \sigma_C^2 & \sigma_C^2      
  \end{pmatrix} +
    \begin{pmatrix}
      \sigma_D^2 & \Delta_{7}\sigma_D^2 \\
      \Delta_{7}\sigma_D^2 & \sigma_D^2      
  \end{pmatrix} +
  \begin{pmatrix}
    \sigma_E^2 & 0 \\
    0 & \sigma_E^2      
  \end{pmatrix}
\end{gather*}



\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
dd <- na.omit(twinbmi)
l0 <- twinlm(bmi ~ age+gender, data=dd, 
	    DZ="DZ", zyg="zyg", id="tvparnr", type="sat")

\end{lstlisting}

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
l <- twinlm(bmi ~ ns(age,1)+gender, data=twinbmi, 
	   DZ="DZ", zyg="zyg", id="tvparnr", type="cor", missing=TRUE)
summary(l)
\end{lstlisting}

\begin{verbatim}
____________________________________________________
Group 1
                        Estimate Std. Error   Z value Pr(>|z|)
Regressions:                                                  
   bmi.1~ns(age, 1).1    4.16937    0.16669  25.01334   <1e-12
   bmi.1~gendermale.1    1.41160    0.07284  19.37839   <1e-12
Intercepts:                                                   
   bmi.1                22.53618    0.07296 308.87100   <1e-12
Additional Parameters:                                        
   log(var)              2.44580    0.01425 171.68256   <1e-12
   atanh(rhoMZ)          0.78217    0.02290  34.16186   <1e-12
____________________________________________________
Group 2
                        Estimate Std. Error   Z value Pr(>|z|)
Regressions:                                                  
   bmi.1~ns(age, 1).1    4.16937    0.16669  25.01334   <1e-12
   bmi.1~gendermale.1    1.41160    0.07284  19.37839   <1e-12
Intercepts:                                                   
   bmi.1                22.53618    0.07296 308.87100   <1e-12
Additional Parameters:                                        
   log(var)              2.44580    0.01425 171.68256   <1e-12
   atanh(rhoDZ)          0.29924    0.01848  16.19580   <1e-12

                       Estimate 2.5%    97.5%  
Correlation within MZ: 0.65395  0.62751 0.67889
Correlation within DZ: 0.29061  0.25712 0.32341

'log Lik.' -29020.12 (df=6)
AIC: 58052.24 
BIC: 58093.29
\end{verbatim}

A formal test of genetic effects can be obtained by comparing the MZ and DZ correlation:
\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
estimate(l,contr(5:6,6))
\end{lstlisting}

\begin{verbatim}
                          Estimate Std.Err   2.5%  97.5%   P-value
[atanh(rhoMZ)@1] - [a....   0.4829 0.04176 0.4011 0.5648 6.325e-31

 Null Hypothesis: 
  [atanh(rhoMZ)@1] - [atanh(rhoDZ)@3] = 0
\end{verbatim}


\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
l <- twinlm(bmi ~ ns(age,1)+gender, data=twinbmi, 
	   DZ="DZ", zyg="zyg", id="tvparnr", type="cor", missing=TRUE)
summary(l)
\end{lstlisting}

\begin{verbatim}
____________________________________________________
Group 1
                        Estimate Std. Error   Z value Pr(>|z|)
Regressions:                                                  
   bmi.1~ns(age, 1).1    4.16937    0.16669  25.01334   <1e-12
   bmi.1~gendermale.1    1.41160    0.07284  19.37839   <1e-12
Intercepts:                                                   
   bmi.1                22.53618    0.07296 308.87100   <1e-12
Additional Parameters:                                        
   log(var)              2.44580    0.01425 171.68256   <1e-12
   atanh(rhoMZ)          0.78217    0.02290  34.16186   <1e-12
____________________________________________________
Group 2
                        Estimate Std. Error   Z value Pr(>|z|)
Regressions:                                                  
   bmi.1~ns(age, 1).1    4.16937    0.16669  25.01334   <1e-12
   bmi.1~gendermale.1    1.41160    0.07284  19.37839   <1e-12
Intercepts:                                                   
   bmi.1                22.53618    0.07296 308.87100   <1e-12
Additional Parameters:                                        
   log(var)              2.44580    0.01425 171.68256   <1e-12
   atanh(rhoDZ)          0.29924    0.01848  16.19580   <1e-12

                       Estimate 2.5%    97.5%  
Correlation within MZ: 0.65395  0.62751 0.67889
Correlation within DZ: 0.29061  0.25712 0.32341

'log Lik.' -29020.12 (df=6)
AIC: 58052.24 
BIC: 58093.29
\end{verbatim}

\section*{Twin analysis, censored outcomes}
\label{sec:org4b5f762}

\section*{Twin analysis,  binary traits}
\label{sec:org4c17447}


\section*{Time to event}
\label{sec:org843c812}


\bibliography{mets}
\bibliographystyle{plain}
\end{document}