---
title: "TREC & KREC of twins: Decomposing the covariance matrix"
author:
  - Henrique Aparecido Laureano^[Instituto de Pesquisa
                                 Pelé Pequeno Príncipe,
                                 Curitiba, Paraná, Brazil]
  - Wagner Hugo Bonat^[Laboratório de
                       Estatística e Geoinformação,
                       Universidade Federal do Paraná,
                       Curitiba, Paraná, Brazil]
  - Stéfanne Maria Jeha Bortoletto^[Faculdades Pequeno Príncipe &
                                    Instituto de Pesquisa
                                    Pelé Pequeno Príncipe,
                                    Curitiba, Paraná, Brazil]
  - Carolina Cardoso de Mello Prando^[Faculdades
                                      Pequeno Príncipe &
                                      Instituto de Pesquisa
                                      Pelé Pequeno Príncipe,
                                      Curitiba, Paraná, Brazil]
date: September, 2022

abstract: Multivariate linear *regression* models built to handle
          multiple responses and with genetic and environmental
          interrelations of twins data [@bonat&hjelmborg] were applied
          to understand the dynamics from the TREC and KREC measures.

bibliography: /home/laureano/Git/consu/references.bib
format:
  pdf:
    echo: false
    warning: false
---

**Keywords**:

```{r}
#| label: pkgs

## setwd('papers/twins/')

if(!requireNamespace('pacman', quietly=TRUE)) install.packages('pacman')

pacman::p_load(readxl, dplyr, tidyr, haven, ggplot2, patchwork, tibble,
               kableExtra, Matrix, ggfortify)

## if(!requireNamespace('devtools', quietly=TRUE))
##     install.packages('devtools')
## devtools::install_github('wbonat/mglm4twin')
library(mglm4twin)

```

```{r}
#| label: data

path <- '/home/laureano/Git/carolina-I4P/trec_and_krec/'
file <- 'Dados 50k amostras.xlsx'
dat  <- readxl::read_xlsx(paste0(path, file), sheet='Gemeos')|>
    dplyr::mutate(parto=dplyr::recode(parto,
                                      'cesariana'='Cesariana',
                                      'NI'       ='Não informado'),
                  ID2  =gsub('.$', '', ID))|>
    dplyr::select(
               ID2, TREC, KREC, peso, id_gest, parto, sexo, Zigocidade
           )|>
    dplyr::filter(!Zigocidade %in% 'ERRADO' & !ID2 %in% 'PG74')|>
    dplyr::arrange(Zigocidade)|>
    dplyr::group_by(ID2)|>
    dplyr::mutate(Twin_pair =rep(c(1, 2), times=length(unique(ID2))),
                  Zigocidade=dplyr::recode(Zigocidade,
                                           'diz'='DZ', 'monoz'='MZ'))|>
    dplyr::ungroup()|>
    dplyr::filter(TREC <= 215 & KREC <= 215 & !is.na(peso))

id2drop <- dat|>
    dplyr::count(ID2)|>
    dplyr::filter(n == 1)|>
    dplyr::pull(ID2)

dat <- dat|>
    dplyr::filter(!ID2 %in% id2drop)|>
    dplyr::select(!ID2)

dat$PC1 <- prcomp(~ peso + id_gest,
                  data=dat, center=TRUE, scale.=TRUE)$x[, 'PC1']|>
                                                    as.numeric()
n   <- nrow(dat)
N   <- dat|>
    dplyr::group_by(Zigocidade)|>
    dplyr::summarise(n=dplyr::n()/2)|>
    tidyr::spread(Zigocidade, n)

```

# Introduction

# Data

We have the TREC and KREC measures of `r n/2` pairs of twins. From these
`r n` twins, we have five covariates, as follow:

+ Birth weight (grams);
+ Gestational age (weeks);
+ Type of childbirth (normal or cesarean);
+ Sex (male of female);
+ Zygosity (MZ: monozygotic or DZ: dizygotic).

Zygosity is a key term since we need to inform how many twin pairs are
monozygotic and dizygotic. Besides these five covariates, we may also
consider the twin covariate (1 or 2), to verify a randomness
supposition.

```{r}
#| label: fig-TREC-and-KREC-plots
#| fig-width: 10
#| fig-height: 9
#| fig-cap: 'Graph A: TREC and KREC boxplots with their means in
#|           asterisks; Graph B: TREC and KREC scatterplot with tendency
#|           curve in solid black and 95% confidence band in gray; Graph
#|           C: Twins scatterplots per TREC and KREC with tendency
#|           curves in solid black and 95% confidence bands in gray.'
#| fig-pos: 'H'

p1 <- dat|>
    tidyr::pivot_longer(TREC:KREC, names_to='var')|>
    dplyr::mutate(var=haven::as_factor(var))|>
    ggplot()+
    aes(x=var, y=value)+
    scale_y_continuous(breaks=seq(20, 200, by=20))+
    geom_hline(yintercept=seq(20, 200, by=20),
               color     ='gray',
               linetype  ='dashed')+
    geom_boxplot(outlier.shape=21, outlier.size=2)+
    stat_summary(fun='mean', geom='point', size=3, shape=8)+
    labs(x=NULL, y=NULL)+
    theme_classic(base_size=14)+
    theme(axis.text.x=element_text(
              face  ='bold',
              size  =14,
              margin=ggplot2::unit(c(t=3, r=0, b=0, l=0), 'mm')
          ))

p2 <- dat|>
    ggplot()+
    aes(x=TREC, y=KREC)+
    scale_x_continuous(breaks=seq(20, 200, by=20))+
    scale_y_continuous(breaks=seq(20, 200, by=20))+
    geom_vline(xintercept=seq(20, 200, by=20),
               color     ='gray',
               linetype  ='dashed')+
    geom_hline(yintercept=seq(20, 200, by=20),
               color     ='gray',
               linetype  ='dashed')+
    geom_point(size=3, shape=21)+
    geom_smooth(method ='loess',
                formula=y ~ x, se=TRUE, size=1, color='black')+
    labs(x='TREC', y='KREC')+
    theme_classic(base_size=14)+
    theme(axis.title.x=element_text(face='bold'),
          axis.title.y=element_text(face='bold'))

p3 <- dplyr::bind_cols(
                 by=rep(c('TREC', 'KREC'), each=n/2), 
                 dplyr::bind_rows(
                            dat|>
                            dplyr::select(Twin_pair, TREC)|>
                            tidyr::pivot_wider(names_from =Twin_pair,
                                               values_from=TREC,
                                               values_fn  =list)|>
                            tidyr::unnest(c(`1`, `2`)), 
                            dat|>
                            dplyr::select(Twin_pair, KREC)|>
                            tidyr::pivot_wider(names_from =Twin_pair,
                                               values_from=KREC,
                                               values_fn  =list)|>
                            tidyr::unnest(c(`1`, `2`))
                        )
             )|>
    dplyr::mutate(by=haven::as_factor(by))|>
    ggplot()+
    aes(x=`1`, y=`2`)+
    scale_x_continuous(breaks=seq(20, 200, by=20))+
    scale_y_continuous(breaks=seq(0, 200, by=20))+
    geom_vline(xintercept=seq(20, 200, by=20),
               color     ='gray',
               linetype  ='dashed')+
    geom_hline(yintercept=seq(0, 200, by=20),
               color     ='gray',
               linetype  ='dashed')+
    geom_point(size=3, shape=21)+
    geom_smooth(method ='loess',
                formula=y ~ x, se=TRUE, size=1, color='black')+
    facet_wrap(~ by)+
    labs(x='Twin 1', y='Twin 2')+
    theme_classic(base_size=14)+
    theme(axis.title.x=element_text(face  ='bold',
                                    margin=ggplot2::unit(c(
                                                        t=3, r=0,
                                                        b=0, l=0
                                                    ), 'mm')),
          axis.title.y=element_text(face='bold'),
          strip.text.x=element_text(face='bold', size=14))

layout <- '
AABBBBB
CCCCCCC
'
p1 + p2 + p3 +
    patchwork::plot_annotation(tag_levels='A') +
    patchwork::    plot_layout(design=layout)

```

In @fig-TREC-and-KREC-plots and the C panel of @fig-appendix we have
some plots to see how our two responses behave. We can see that the KREC
measures are asymmetric with a stronger left tail. Both are positively
correlated with an almost linear correspondence. By crossing twins, we do
not see any clear leakage. For the TREC we have a mean (standard
deviation, sd) of
`r round(mean(dat$TREC), 2)` (`r round(sd(dat$TREC), 2)`) and for the
KREC a mean (sd) of
`r round(mean(dat$KREC), 2)` (`r round(sd(dat$KREC), 2)`).

In @fig-covariate-plots we have the plots of the covariates. We have a
mean (sd) birth weight of `r round(mean(dat$peso), 2)`
(`r round(sd(dat$peso), 2)`) grams and a mean (sd) gestational age of
`r round(mean(dat$id_gest), 2)` (`r round(sd(dat$id_gest), 2)`) weeks.
88.2% of the twins were born from cesarean childbirth, 74.5% are
zygotes, and 60.9% are from the male sex.

```{r}
#| label: fig-covariate-plots
#| fig-width: 10
#| fig-height: 4
#| fig-cap: Graphs A and B are the boxplots of the numerical variables
#|          with their means in asterisks; Graphs C, D and E are the
#|          barplots of the categorial variables with their levels
#|          frequencies and respective sample percentages.
#| fig-pos: 'H'

p1 <- dat|>
    dplyr::rename(Weight='peso')|>
    ggplot()+
    aes(x=Weight, y=' ')+
    scale_x_continuous(breaks=seq(1e3, 3800, by=400))+
    geom_vline(xintercept=seq(1e3, 3800, by=400),
               color     ='gray',
               linetype  ='dashed')+
    geom_boxplot(outlier.shape=21, outlier.size=3)+
    stat_summary(fun='mean', geom='point', size=4, shape=8)+
    labs(x='Grams', y=NULL, title='Birth weight')+
    theme_classic(base_size=14)+
    theme(plot.title  =element_text(face='bold'),
          axis.title.x=element_text(
              face  ='bold',
              margin=ggplot2::unit(c(t=3, r=0, l=0, b=0), 'mm')
          ))

p2 <- dat|>
    dplyr::rename(`Gestational age`='id_gest')|>
    ggplot()+
    aes(x=`Gestational age`, y=' ')+
    scale_x_continuous(breaks=seq(27, 39, by=1))+
    geom_vline(xintercept=seq(27, 39, by=1),
               color     ='gray',
               linetype  ='dashed')+
    geom_boxplot(outlier.shape=21, outlier.size=3)+
    stat_summary(fun='mean', geom='point', size=4, shape=8)+
    labs(x='Weeks', y=NULL, title='Gestational age')+
    theme_classic(base_size=14)+
    theme(plot.title  =element_text(face='bold'),
          axis.title.x=element_text(
              face  ='bold',
              margin=ggplot2::unit(c(t=3, r=0, l=0, b=0), 'mm')
          ))

p3 <- dat|>
    dplyr::mutate(parto=dplyr::recode(parto, Cesariana='Cesarean'))|>
    ggplot()+
    aes(y=parto)+
    geom_bar(fill='white', color='black')+
    geom_text(stat='count',
              aes(label=paste0(..count..,
                               ' (', 100*round(..count../n, 3), '%)'
                               )),
              hjust=-0.1, fontface='bold')+
    scale_x_continuous(expand=c(0, 0), breaks=NULL, limits=c(0, 155))+
    labs(x=NULL, y=NULL, title='Type of childbirth')+
    theme_classic(base_size=14)+
    theme(axis.text.y=element_text(face='bold', size=14),
          plot.title =element_text(face='bold'))

p4 <- dat|>
    dplyr::mutate(sexo=dplyr::recode(sexo,
                                     Feminino ='Female',
                                     Masculino='Male'))|>
    ggplot()+
    aes(y=sexo)+
    geom_bar(fill='white', color='black')+
    geom_text(stat='count',
              aes(label=paste0(..count..,
                               ' (', 100*round(..count../n, 3), '%)'
                               )),
              hjust=-0.1, fontface='bold')+
    scale_x_continuous(expand=c(0, 0), breaks=NULL, limits=c(0, 110))+
    labs(x=NULL, y=NULL, title='Sex')+
    theme_classic(base_size=14)+
    theme(axis.text.y=element_text(face='bold', size=14),
          plot.title =element_text(face='bold'))

p5 <- dat|>
    ggplot()+
    aes(y=Zigocidade)+
    geom_bar(fill='white', color='black')+
    geom_text(stat='count',
              aes(label=paste0(..count..,
                               ' (', 100*round(..count../n, 3), '%)'
                               )),
              hjust=-0.1, fontface='bold')+
    scale_x_continuous(expand=c(0, 0), breaks=NULL, limits=c(0, 135))+
    labs(x=NULL, y=NULL, title='Zygosity')+
    theme_classic(base_size=14)+
    theme(axis.text.y=element_text(face='bold', size=14),
          plot.title =element_text(face='bold'))

(p1 + p2) / (p3 + p4 + p5) + patchwork::plot_annotation(tag_levels='A')

```

# Modeling framework

The statistical analysis was performed through the R [@R22] language and
environment for statistical computing. Besides the model fitting
performed through the `{mglm4twin}` [@mglm4twin; @bonat&hjelmborg;
@mcglm.seriesc; @mcglm] package, the leading R packages used are:
`{dplyr}` [@dplyr1.0.9], `{tidyr}` [@tidyr1.2.0], `{ggplot2}`
[@ggplot2.16], `patchwork` [@patchwork1.1.1], `{kableExtra}`
[@kableExtra1.3.4], `{Matrix}` [@matrix1.4.1] and `{ggfortify}`
[@ggfortify.16].

# Results

We model the TREC and KREC measures in a bivariate fashion. We use a
model called ACE to study heritability and the genetic and environmental
relationships. Basically, we decompose the covariance matrix into three:

+ A: genetic/heritability effect or component;
+ C: common environment effect or component;
+ E: unique environment effect or component.

```{r}
#| label: bivariate-models
#| fig-height: 10
#| fig-width: 10
#| results: false

## mean linear predictors ----------------------------------------------
## form_TREC <- TREC ~ PC1 + parto + sexo + Zigocidade * Twin_pair
## form_KREC <- KREC ~ PC1 + parto + sexo + Zigocidade * Twin_pair

## ACE model -----------------------------------------------------------

## biv_ACE <- mglm4twin::mt_twin(N_DZ=N$DZ,
##                               N_MZ=N$MZ, n_resp=2, model='ACE')
## GAUSSIAN 
## fitACEbiv <-
##     mglm4twin::mglm4twin(
##                    linear_pred      =c(form_TREC, form_KREC),
##                    matrix_pred      =biv_ACE,
##                    control_algorithm=list(tuning=0.75, max_iter=40),
##                    link             =rep('identity', 2),
##                    variance         =rep('constant', 2),
##                    power_fixed      =rep(TRUE, 2),
##                    data             =dat
##                )
## summary(fitACEbiv, model='ACE', biometric=TRUE)

## TWEEDIE
## control_algorithm=list(tuning=0.01, max_iter=750),
## link             =rep(    'log', 2),
## variance         =rep('tweedie', 2),
## power_fixed      =rep(FALSE, 2),

## GAMMA
## control_algorithm=list(tuning=0.01, max_iter=1e3),
## link             =rep('inverse', 2),
## variance         =rep('tweedie', 2),
## power_fixed      =rep(TRUE, 2),

## cbind(Model=c('Gaussian', 'Tweedie', 'Gamma'),
##       rbind(mglm4twin::gof(fitACEbiv),
##             mglm4twin::gof(fitACEbiv_tw),
##             mglm4twin::gof(fitACEbiv_ga)))

## par(mfrow=c(4, 2))
## 
## hist(      dat$TREC            );hist(      dat$KREC              )
## hist(resid(fitACEbiv   )[1:110]);hist(resid(fitACEbiv   )[111:220])
## hist(resid(fitACEbiv_tw)[1:110]);hist(resid(fitACEbiv_tw)[111:220])
## hist(resid(fitACEbiv_ga)[1:110]);hist(resid(fitACEbiv_ga)[111:220])
## 
## layout(1)

## ADE model -----------------------------------------------------------

## biv_ADE <- mglm4twin::mt_twin(N_DZ=N$DZ,
##                               N_MZ=N$MZ, n_resp=2, model='ADE')
## fitADEbiv <-
##     mglm4twin::mglm4twin(
##                    linear_pred      =c(form_TREC, form_KREC),
##                    matrix_pred      =biv_ADE,
##                    control_algorithm=list(tuning=0.75, max_iter=40),
##                    link             =rep('identity', 2),
##                    variance         =rep('constant', 2),
##                    power_fixed      =rep(TRUE, 2),
##                    data             =dat
##                )
## summary(fitADEbiv, model='ADE', biometric=TRUE)

## ACE vs ADE `gof`s ---------------------------------------------------

## cbind(Model=c('ACE', 'ADE'),
##       rbind(mglm4twin::gof(fitACEbiv), mglm4twin::gof(fitADEbiv)))

## AE model ------------------------------------------------------------

biv_AE <- mglm4twin::mt_twin(N_DZ=N$DZ,
                             N_MZ=N$MZ, n_resp=2, model='AE')
## fitAEbiv <-
##     mglm4twin::mglm4twin(
##                    linear_pred      =c(form_TREC, form_KREC),
##                    matrix_pred      =biv_AE,
##                    control_algorithm=list(tuning=0.75),
##                    link             =rep('identity', 2),
##                    variance         =rep('constant', 2),
##                    power_fixed      =rep(TRUE, 2),
##                    data             =dat
##                )
## summary(fitAEbiv, model='AE', biometric=TRUE)

## ACE vs AE `gof`s ----------------------------------------------------

## cbind(Model=c('ACE', 'AE'), rbind(mglm4twin::gof(fitACEbiv),
##                                   mglm4twin::gof( fitAEbiv)))

## FINAL MODEL ---------------------------------------------------------

## mean linear predictors ----------------------------------------------
form_TREC <- TREC ~ parto
form_KREC <- KREC ~ parto

fitAEbivf <-
    mglm4twin::mglm4twin(linear_pred      =c(form_TREC, form_KREC),
                         matrix_pred      =biv_AE,
                         control_algorithm=list(tuning=0.75), 
                         link             =rep('identity', 2),
                         variance         =rep('constant', 2),
                         power_fixed      =rep(TRUE, 2),
                         data             =dat)

## summary(fitAEbivf, model='AE', biometric=TRUE)

## cbind(Model=c('AE initial', 'AE final'),
##       rbind(mglm4twin::gof(fitAEbiv), mglm4twin::gof(fitAEbivf)))

fitAEbiv <- fitAEbivf
tt       <- summary(fitAEbiv, model='AE', biometric=TRUE)

```

In the mean structure we tested for the effect of six covariates. With a
bivariate model, we are able to handle the TREC and KREC measures
jointly. We started with the ACE model, however, the C component
presented to be non-significant. Thus, we stayed with an AE model and
from this performed a covariates selection procedure.

```{r}
#| label: tbl-AE-bivariate-model-mean-summary
#| tbl-cap: "Summary of the mean structure of the AE bivariate model."

dplyr::bind_rows(
           cbind(
               ' '=rownames(tt[[1]][[1]]), tibble::tibble(tt[[1]][[1]])
           ),
           cbind(
               ' '=rownames(tt[[2]][[1]]), tibble::tibble(tt[[2]][[1]])
           )
       )|>
    dplyr::mutate(` `=c('Cesarean childbirth (intercept)',
                        'Normal childbirth (add to intercept)',
                        'Cesarean childbirth (intercept)',
                        'Normal childbirth (add to intercept)'))|>
    dplyr::rename(Estimate ='Estimates',
                  SD       ='std.error',
                  'z-value'='z value')|>
    kableExtra::            kbl(booktabs=TRUE,
                                digits  =4,
                                caption =NA,
                                position='h')|>
    kableExtra::kable_classic_2(full_width=FALSE)|>
    kableExtra::    column_spec(5, bold=TRUE)|>
    kableExtra::      pack_rows('TREC', 1, 2)|>
    kableExtra::      pack_rows('KREC', 3, 4)

```

<!-- (`r round(var(dat$TREC), 4)`) -->
<!-- (`r round(var(dat$KREC), 4)`) -->

From the six tested covariates we end up with two. For both responses,
only the type of chilbirth presented to be as statistically
significant. In @tbl-AE-bivariate-model-mean-summary we have a summary
of the estimated effects. About the covariance structure, what we do is
decompose the TREC variance (2167.713), KREC variance (1820.538), and
their covariance (`r round(cov(dat$TREC, dat$KREC), 4)`) in the manner
described in @tbl-AE-bivariate-model-covariance.

```{r}
#| label: tbl-AE-bivariate-model-covariance
#| tbl-cap:
#|   "Summary of the covariance structure of the AE bivariate model."

tibble::tibble(tt$Dispersion)|>
    dplyr::mutate(' '=dplyr::recode(Parameters,
                                    'E1' ='TREC',
                                    'E2' ='KREC',
                                    'E12'='TREC & KREC',
                                    'A1' ='TREC',
                                    'A2' ='KREC',
                                    'A12'='TREC & KREC'),
                  Percentage=Estimates/rep(c(var(dat$TREC),
                                             var(dat$KREC),
                                             cov(dat$TREC, dat$KREC)),
                                           2)*100)|>
    select(!Parameters)|>
    dplyr::relocate(       ` `, .before='Estimates')|>
    dplyr::relocate(Percentage, .before='z value')|>
    dplyr::rename(Estimate ='Estimates',
                  SD       ='std.error',
                  'z-value'='z value')|>
    kableExtra::            kbl(booktabs=TRUE,
                                digits  =4,
                                caption =NA,
                                position='h')|>
    kableExtra::kable_classic_2(full_width=FALSE)|>
    kableExtra::    column_spec(6, bold=TRUE)|>
    kableExtra::      pack_rows('Environment\ncomponent (E)', 1, 3)|>
    kableExtra::      pack_rows(    'Genetic\ncomponent (A)', 4, 6)

```

At @tbl-AE-bivariate-model-genetic-components we have the environment
and heritability measures. The heritability is statistically significant
in all scenarios, marginally; jointly; and its correlation.

```{r}
#| label: tbl-AE-bivariate-model-genetic-components
#| tbl-cap:
#|   "Summary of the genetic components of the AE bivariate model."

tibble::tibble(
            cbind(' '=rep(c('TREC',
                            'KREC', 'TREC & KREC', 'Correlation'),
                          2), 
                  rbind(tt$Env_main,
                        tt$Env_cross,
                        c(0.4210331, 0.3306534, 1.273336, 0.2028987), 
                        tt$A_main,
                        tt$A_cross,
                        c(0.3201396, 0.1929394, 1.659276, 0.09706021))
                  )
        )|>
    dplyr::rename(Estimate ='Estimates',
                  SD       ='std.error',
                  'z-value'='z value')|>
    kableExtra::            kbl(booktabs=TRUE,
                                digits  =8,
                                caption =NA,
                                position='h')|>
    kableExtra::kable_classic_2(full_width=FALSE)|>
    kableExtra::column_spec    (5, bold=TRUE)|>
    kableExtra::pack_rows      ('Environmentability', 1, 4)|>
    kableExtra::pack_rows      (      'Heritability', 5, 8)

```

```{r}
#| label: univariate-models
#| results: false
#| eval: false

cov_TREC   <- list(formE=~ sexo + parto,
                   formA=~ sexo,
                   formAE=~ PC1 + sexo + parto)
uni_AEtrec <- mglm4twin::mt_twin(N_DZ   =N$DZ,
                                 N_MZ   =N$MZ,
                                 n_resp =1,
                                 model  ='AE',
                                 formula=cov_TREC,
                                 dat    =dat)
fitAEtrec  <-
    mglm4twin::mglm4twin(linear_pred      =c(form_TREC),
                         matrix_pred      =uni_AEtrec,
                         control_algorithm=list(tuning=0.25, max_iter=1e3), 
                         link             ='identity',
                         variance         ='constant',
                         power_fixed      =TRUE,
                         data             =dat)

fitAEtrec$Covariance
sqrt(diag(fitAEtrec$vcov))

tt         <- rbind(cbind(fitAEtrec$Covariance, 
                          sqrt(diag(fitAEtrec$vcov))[-c(1:3)])[1, ], 
                    cbind(fitAEtrec$Covariance, 
                          sqrt(diag(fitAEtrec$vcov))[-c(1:3)])[2:3, ])

#| label: tbl-TREC-AE-model-covariance
#| tbl-cap: "Summary of the covariance structure of the AE TREC model."
#| eval: false

tibble::tibble(' '       =c('Environment component (E)',
                            'Cesarean childbirth',
                            'Normal childbirth'),
               Estimate  =tt[, 1],
               SD        =tt[, 2],
               Percentage=Estimate/rep(var(dat$TREC), 3)*100)|>
    kableExtra::            kbl(booktabs=TRUE,
                                digits  =4,
                                caption =NA,
                                position='h')|>
    kableExtra::kable_classic_2(full_width=FALSE)|>
    kableExtra::pack_rows      ('Genetic component (A)', 2, 3)

#| label: KREC-model
#| results: false
#| eval: false

cov_KREC   <- list(formE=~ parto + Zigocidade, formA=~ parto)
uni_AEkrec <- mglm4twin::mt_twin(N_DZ   =N$DZ,
                                 N_MZ   =N$MZ,
                                 n_resp =1,
                                 model  ='AE',
                                 formula=cov_KREC,
                                 dat    =dat)
fitAEkrec  <-
    mglm4twin::mglm4twin(
                   linear_pred      =c(form_KREC),
                   matrix_pred      =uni_AEkrec,
                   control_algorithm=list(tuning=0.5, max_iter=80),
                   link             ='identity',
                   variance         ='constant',
                   power_fixed      =TRUE,
                   data             =dat
               )
tt <- rbind(cbind(fitAEkrec$Covariance, 
                  sqrt(diag(fitAEkrec$vcov))[-c(1:2)])[1:3, ], 
            cbind(fitAEkrec$Covariance, 
                  sqrt(diag(fitAEkrec$vcov))[-c(1:2)])[4:5, ])

#| label: tbl-KREC-AE-model-covariance
#| tbl-cap: "Summary of the covariance structure of the AE KREC model."
#| eval: false

tibble::tibble(' '       =c('Cesarean childbirth & DZ zygosity',
                            'Normal childbirth',
                            'MZ zygosity',
                            'Cesarean childbirth', 
                            'Normal childbirth'),
               Estimate  =tt[, 1],
               SD        =tt[, 2],
               Percentage=Estimate/rep(var(dat$KREC), 5)*100)|>
    kableExtra::            kbl(booktabs=TRUE,
                                digits  =4,
                                caption =NA,
                                position='h')|>
    kableExtra::kable_classic_2(full_width=FALSE)|>
    kableExtra::pack_rows      ('Environment component (E)', 1, 3)|>
    kableExtra::pack_rows      (    'Genetic component (A)', 4, 5)

```

# References

::: {#refs}
:::

# Appendix

```{r}
#| label: fig-appendix
#| fig-height: 9
#| fig-width: 10
#| fig-cap: 'Graph A: Gestational age and weight scatterplot with
#|           tendency curve in solid black and 95% confidence band in
#|           gray; Graph B: Gestational age and weight biplot from a
#|           Principal Components Analysis (PCA); Graph C: TREC and KREC
#|           histograms with density curve estimations in solid black.'
#| fig-pos: 'H'

p1 <- dat|>
    ggplot()+
    aes(x=id_gest, y=peso)+
    scale_x_continuous(breaks=seq(27, 39, by=1))+
    scale_y_continuous(breaks=seq(750, 3750, by=250))+
    geom_vline(xintercept=seq(27, 39, by=1),
               color     ='gray',
               linetype  ='dashed')+
    geom_hline(yintercept=seq(750, 3750, by=250),
               color     ='gray',
               linetype  ='dashed')+
    geom_point(pch=21, size=3)+
    geom_smooth(method ='loess',
                formula=y ~ x, se=TRUE, color='black', size=1)+
    labs(x='Gestational age\n(Weeks)', y='Weight\n(Grams)')+
    theme_classic(base_size=14)+
    theme(axis.title.x=element_text(face  ='bold',
                                    margin=ggplot2::unit(c(
                                                        t=3, r=0,
                                                        b=0, l=0
                                                    ), 'mm')),
          axis.title.y=element_text(face='bold'))

dat.pca <- dat|>
    dplyr::select(peso, id_gest)|>
    dplyr::rename(Weight='peso', Age='id_gest')

p2 <- autoplot(prcomp(dat.pca, center=TRUE, scale.=TRUE),
               data=dat.pca, shape=21, size=3, 
               loadings             =TRUE,
               loadings.colour      ='black', 
               loadings.label       =TRUE,
               loadings.label.size  =4.5,
               loadings.label.colour='black')+
    geom_vline(xintercept=seq(-0.2, 0.2, by=0.1),
               color     ='gray',
               linetype  ='dashed')+
    geom_hline(yintercept=seq(-0.2, 0.2, by=0.1),
               color     ='gray',
               linetype  ='dashed')+
    theme_classic(base_size=14)+
    theme(axis.title.x=element_text(face='bold'),
          axis.title.y=element_text(face='bold'))

p3 <- dat|>
    tidyr::pivot_longer(c(TREC, KREC))|>
    dplyr::mutate(name=haven::as_factor(name))|>
    ggplot()+
    aes(x=value)+
    scale_x_continuous(breaks=seq(20, 220, by=20))+
    scale_y_continuous(breaks=NULL)+
    geom_vline(xintercept=seq(20, 220, by=10),
               color     ='gray',
               linetype  ='dashed')+
    geom_histogram(aes(y=..density..),
                   binwidth=10, color='black', fill='white')+
    geom_density(size=1)+
    facet_wrap(~ name)+
    labs(x=NULL, y=NULL)+
    theme_classic(base_size=14)+
    theme(strip.text.x=element_text(face='bold', size=14))

layout <- '
AB
CC
'
p1 + p2 + p3 +
    patchwork::plot_annotation(tag_levels='A') +
    patchwork::    plot_layout(design=layout)

```
