Mid-parent–offspring heritability estimate

Parent-Offspring-h2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mid-parent–offspring heritability estimate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Demonstrate how heritability can be estimated from regressing offspring on mid-parental phenotype, using symbolic calculations with statistical functions in [SymPy](https://www.sympy.org/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy.stats import P, E, variance, covariance, Variance, Covariance, Normal\n",
    "from sympy import Symbol, symbols, sqrt, simplify"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define additive genetic variance $V_a$ and environment variance $V_e$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "Va = Symbol('V_a', positive = True)\n",
    "Ve = Symbol('V_e', positive = True)\n",
    "Vp = Va + Ve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model of maternal, paternal, and offspring phenotypes ($Y_m, Y_p, Y_o$) based on an additive genetic model with a transmitted ($a_{m_T}, a_{p_T}$) and nontransmitted ($a_{m_{NT}}, a_{p_{NT}}$) effects from parents to offspring, and environmental effects ($e_m, e_p, e_o$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "amT= Normal(\"a_{m_T}\", 0, sqrt(Va/2))\n",
    "amNT = Normal(\"a_{m_{NT}}\", 0, sqrt(Va/2))\n",
    "apT = Normal(\"a_{p_T}\", 0, sqrt(Va/2))\n",
    "apNT = Normal(\"a_{p_{NT}}\", 0, sqrt(Va/2))\n",
    "\n",
    "em = Normal(\"e_m\", 0, sqrt(Ve))\n",
    "ep = Normal(\"e_p\", 0, sqrt(Ve))\n",
    "eo = Normal(\"e_o\", 0, sqrt(Ve))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle a_{m_T} + a_{p_T} + e_{o}$"
      ],
      "text/plain": [
       "a_{m_T} + a_{p_T} + e_o"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ym = amT + amNT + em\n",
    "Yp = apT + apNT + ep\n",
    "Yo = amT + apT + eo\n",
    "Yo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check variance of offspring phenotype (should be $\\mathrm{var}(Y) = V_a + V_e$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle V_{a} + V_{e}$"
      ],
      "text/plain": [
       "V_a + V_e"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "variance(Yo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Covariance between mid-parental phenotype and offspring phenotype, and expand symbolicly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\operatorname{Cov}\\left(a_{m_T} + a_{p_T} + e_{o}, \\frac{a_{m_T}}{2} + \\frac{a_{m_{NT}}}{2} + \\frac{a_{p_T}}{2} + \\frac{a_{p_{NT}}}{2} + \\frac{e_{m}}{2} + \\frac{e_{p}}{2}\\right)$"
      ],
      "text/plain": [
       "Covariance(a_{m_T} + a_{p_T} + e_o, a_{m_T}/2 + a_{m_{NT}}/2 + a_{p_T}/2 + a_{p_{NT}}/2 + e_m/2 + e_p/2)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "covOP = Covariance((Ym + Yp)/2, Yo)\n",
    "covOP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\operatorname{Cov}\\left(a_{m_T}, a_{m_T}\\right)}{2} + \\frac{\\operatorname{Cov}\\left(a_{m_T}, a_{m_{NT}}\\right)}{2} + \\operatorname{Cov}\\left(a_{m_T}, a_{p_T}\\right) + \\frac{\\operatorname{Cov}\\left(a_{m_T}, a_{p_{NT}}\\right)}{2} + \\frac{\\operatorname{Cov}\\left(a_{m_T}, e_{m}\\right)}{2} + \\frac{\\operatorname{Cov}\\left(a_{m_T}, e_{o}\\right)}{2} + \\frac{\\operatorname{Cov}\\left(a_{m_T}, e_{p}\\right)}{2} + \\frac{\\operatorname{Cov}\\left(a_{m_{NT}}, a_{p_T}\\right)}{2} + \\frac{\\operatorname{Cov}\\left(a_{m_{NT}}, e_{o}\\right)}{2} + \\frac{\\operatorname{Cov}\\left(a_{p_T}, a_{p_T}\\right)}{2} + \\frac{\\operatorname{Cov}\\left(a_{p_T}, a_{p_{NT}}\\right)}{2} + \\frac{\\operatorname{Cov}\\left(a_{p_T}, e_{m}\\right)}{2} + \\frac{\\operatorname{Cov}\\left(a_{p_T}, e_{o}\\right)}{2} + \\frac{\\operatorname{Cov}\\left(a_{p_T}, e_{p}\\right)}{2} + \\frac{\\operatorname{Cov}\\left(a_{p_{NT}}, e_{o}\\right)}{2} + \\frac{\\operatorname{Cov}\\left(e_{m}, e_{o}\\right)}{2} + \\frac{\\operatorname{Cov}\\left(e_{o}, e_{p}\\right)}{2}$"
      ],
      "text/plain": [
       "Covariance(a_{m_T}, a_{m_T})/2 + Covariance(a_{m_T}, a_{m_{NT}})/2 + Covariance(a_{m_T}, a_{p_T}) + Covariance(a_{m_T}, a_{p_{NT}})/2 + Covariance(a_{m_T}, e_m)/2 + Covariance(a_{m_T}, e_o)/2 + Covariance(a_{m_T}, e_p)/2 + Covariance(a_{m_{NT}}, a_{p_T})/2 + Covariance(a_{m_{NT}}, e_o)/2 + Covariance(a_{p_T}, a_{p_T})/2 + Covariance(a_{p_T}, a_{p_{NT}})/2 + Covariance(a_{p_T}, e_m)/2 + Covariance(a_{p_T}, e_o)/2 + Covariance(a_{p_T}, e_p)/2 + Covariance(a_{p_{NT}}, e_o)/2 + Covariance(e_m, e_o)/2 + Covariance(e_o, e_p)/2"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "covOP.expand()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute answer for mid-parent–offspring covariance and mid-parent variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{V_{a}}{2}$"
      ],
      "text/plain": [
       "V_a/2"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "covariance((Ym + Yp)/2, Yo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{V_{a}}{2} + \\frac{V_{e}}{2}$"
      ],
      "text/plain": [
       "V_a/2 + V_e/2"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "variance((Ym + Yp)/2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regression of offspring on mid-parent phenotype is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{V_{a}}{V_{a} + V_{e}}$"
      ],
      "text/plain": [
       "V_a/(V_a + V_e)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta = covariance((Ym + Yp)/2, Yo) / variance((Ym + Yp)/2)\n",
    "simplify(beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which yields an estimate of heritability."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}