<htmllang="en"><head><metacharset="UTF-8"/><metaname="viewport"content="width=device-width, initial-scale=1.0"/><title>How to Calculate EPDs · beefblup</title><scriptdata-outdated-warnersrc="../assets/warner.js"></script><linkrel="canonical"href="https://millironx.com/beefblup/how-to-calculate-epds/"/><linkhref="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css"rel="stylesheet"type="text/css"/><linkhref="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.039/juliamono-regular.css"rel="stylesheet"type="text/css"/><linkhref="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/fontawesome.min.css"rel="stylesheet"type="text/css"/><linkhref="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/solid.min.css"rel="stylesheet"type="text/css"/><linkhref="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/brands.min.css"rel="stylesheet"type="text/css"/><linkhref="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.13.11/katex.min.css"rel="stylesheet"type="text/css"/><script>documenterBaseURL=".."</script><scriptsrc="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js"data-main="../assets/documenter.js"></script><scriptsrc="../siteinfo.js"></script><scriptsrc="../../versions.js"></script><linkclass="docs-theme-link"rel="stylesheet"type="text/css"href="../assets/themes/documenter-dark.css"data-theme-name="documenter-dark"data-theme-primary-dark/><linkclass="docs-theme-link"rel="stylesheet"type="text/css"href="../assets/themes/documenter-light.css"data-theme-name="documenter-light"data-theme-primary/><scriptsrc="../assets/themeswap.js"></script></head><body><divid="documenter"><navclass="docs-sidebar"><divclass="docs-package-name"><spanclass="docs-autofit"><ahref="../">beefblup</a></span></div><formclass="docs-search"action="../search/"><inputclass="docs-search-query"id="documenter-search-query"name="q"type="text"placeholder="Search docs"/></form><ulclass="docs-menu"><li><aclass="tocitem"href="../">Home</a></li><liclass="is-active"><aclass="tocitem"href>How to Calculate EPDs</a><ulclass="internal"><li><aclass="tocitem"href="#The-mathematical-model"><span>The mathematical model</span></a></li><li><aclass="tocitem"href="#The-statistical-model:-the-setup"><span>The statistical model: the setup</span></a></li><li><aclass="tocitem"href="#The-statistical-model:-environment-as-fixed-effects"><span>The statistical model: environment as fixed effects</span></a></li><li><aclass="tocitem"href="#The-statistical-model:-genotype-as-random-effect"><span>The statistical model: genotype as random effect</span></a></li><li><aclass="tocitem"href="#Solving-the-equations"><span>Solving the equations</span></a></li><li><aclass="tocitem"href="#Footnotes"><span>Footnotes</span></a></li></ul></li><li><aclass="tocitem"href="../beefblup-cli/">CLI Reference (WIP)</a></li></ul><divclass="docs-version-selector field has-addons"><divclass="control"><spanclass="docs-label button is-static is-size-7">Version</span></div><divclass="docs-selector control is-expanded"><divclass="select is-fullwidth is-size-7"><selectid="documenter-version-selector"></select></div></div></div></nav><divclass="docs-main"><headerclass="docs-navbar"><navclass="breadcrumb"><ulclass="is-hidden-mobile"><liclass="is-active"><ahref>How to Calculate EPDs</a></li></ul><ulclass="is-hidden-tablet"><liclass="is-active"><ahref>How to Calculate EPDs</a></li></ul></nav><divclass="docs-right"><aclass="docs-edit-link"href="https://github.com/MillironX/beefblup/blob/master/docs/src/how-to-calculate-epds.md#"title="Edit on GitHub"><spanclass="docs-icon fab"></span><spanclass="docs-label is-hidden-touch">Edit on GitHub</span></a><aclass="docs-settings-button fas fa-cog"id="documenter-settings-button"href="#"title="Settings"></a><aclass="docs-sidebar-button fa fa-bars is-hidden-desktop"id="documenter-sidebar-button"href="#"></a></div></header><articleclass="content"id="documenter-page"><h1id="How-to-Calculate-EPDs"><aclass="
\end{aligned}\]</p><p>Drat! Every animal I added brings more variables into the system than it eliminates! In fact, since each cow brings in <em>at least</em> one term (<span>$G_n$</span>), I will never be able to write enough equations to solve for <span>$G$</span> numerically. I will have to use a different approach.</p><h2id="The-statistical-model:-the-setup"><aclass="docs-heading-anchor"href="#The-statistical-model:-the-setup">The statistical model: the setup</a><aid="The-statistical-model:-the-setup-1"></a><aclass="docs-heading-anchor-permalink"href="#The-statistical-model:-the-setup"title="Permalink"></a></h2><p>Since I can never solve for <span>$G$</span> directly, I will have to find some way to estimate it. I can switch to a statistical model and solve for <span>$G$</span> that way. The caveat with a statistical model is that there will be some level of error, but so long as we know and can control the level of error, that will be better than not knowing <span>$G$</span> at all.</p><p>Since we're switching into a statistical space, we should also switch the variables we're using. I'll rewrite the first equation as</p><pclass="math-container">\[y = b + u + e\]</p><p>Where:</p><ul><li><span>$y$</span> = Phenotype</li><li><span>$b$</span> = Environment</li><li><span>$u$</span> = Genotype</li><li><span>$e$</span> = Error</li></ul><p>It's not as easy as simply substituting <span>$b$</span> for every <span>$E$</span> that we had above, however. The reason for that is that we must make the assumption that environment is a <strong>fixed effect</strong> and that genotype is a <strong>random effect</strong>. I'll go over why that is later, but for now, understand that we need to transform the environment terms and genotype terms separately.</p><p>We'll start with the environment terms.</p><h2id="The-statistical-model:-environment-as-fixed-effects"><aclass="docs-heading-anchor"href="#The-statistical-model:-environment-as-fixed-effects">The statistical model: environment as fixed effects</a><aid="The-statistical-model:-environment-as-fixed-effects-1"></a><aclass="docs-heading-anchor-permalink"href="#The-statistical-model:-environment-as-fixed-effects"title="Permalink"></a></h2><p>To properly transform the equations, I will have to introduce <span>$b_{mean}$</span> terms in each animal's equation. This is part of the fixed effect statistical assumption, and it will let us obtain a solution.</p><p>Here are the transformed equations:</p><pclass="math-container">\[\begin{aligned}
\end{aligned}\]</p><p>Statistical methods work best in matrix form, so I'm going to convert the set of equations above to a single matrix equation that means the exact same thing.</p><pclass="math-container">\[\begin{bmatrix}
\end{bmatrix}\]</p><p>That's a nice equation, but now my hand is getting tired writing all those <span>$b$</span> terms over and over again, so I'm going to use <ahref="https://www.khanacademy.org/math/precalculus/x9e81a4f98389efdf:matrices/x9e81a4f98389efdf:multiplying-matrices-by-matrices/v/matrix-multiplication-intro">the dot product</a> to condense this down.</p><pclass="math-container">\[\begin{bmatrix}
\end{bmatrix}\]</p><p>That matrix in the middle with all the zeros and ones is called the <strong>incidence matrix</strong>, and essentially reads like a table with each row corresponding to an animal, and each column corresponding to a fixed effect. For brevity, we'll just call it <span>$X$</span>, though. One indicates that the animal and effect go together, and zero means they don't. For our record, we could write a table to go with <span>$X$</span>, and it would look like this:</p><table><tr><thstyle="text-align: left">Animal</th><thstyle="text-align: left">mean</th><thstyle="text-align: left">1990</th><thstyle="text-align: left">1991</th><thstyle="text-align: left">1992</th><thstyle="text-align: left">male</th><thstyle="text-align: left">female</th></tr><tr><tdstyle="text-align: left">1</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td></tr><tr><tdstyle="text-align: left">2</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td></tr><tr><tdstyle="text-align: left">3</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td></tr><tr><tdstyle="text-align: left">4</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td></tr><tr><tdstyle="text-align: left">5</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td></tr><tr><tdstyle="text-align: left">6</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td></tr><tr><tdstyle="text-align: left">7</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td></tr></table><p>Now that we have <span>$X$</span>, we have the ability to start making changes to allow us to solve for <span>$u$</span>. Immediately, we see that <span>$X$</span> is <strong>singular</strong>, meaning it can't be solved directly. We kind of already knew that, but now we can quantify it. We calculate the <ahref="https://math.stackexchange.com/a/2080577">rank of <span>$X$</span></a>, and find that there is only enough information contained in it to solve for 4 variables, which means we need to eliminate two columns.</p><p>There are several ways to effectively eliminate fixed effects in this type of system, but one of the simplest and the most common methods is to declare a <strong>base population</strong>, and lump the fixed effects of animals within the base population into the mean fixed effect. Note that it is possible to declare a base population that has no animals in it, but that gives weird results. For this example, we'll follow the convention built into <code>beefblup</code> and pick the last occuring form of each variable.</p><h3id="Base-population"><aclass="docs-heading-anchor"href="#Base-population">Base population</a><aid="Base-population-1"></a><aclass="docs-heading-anchor-permalink"href="#Base-population"title="Permalink"></a></h3><ul><li><strong>Year</strong>: 1992</li><li><strong>Sex</strong>: Female</li></ul><p>Now in order to use the base population, we simply drop the columns representing conformity w
\end{bmatrix}\]</p><p>And the table for humans to understand:</p><table><tr><thstyle="text-align: left">Animal</th><thstyle="text-align: left">mean</th><thstyle="text-align: left">1990</th><thstyle="text-align: left">1991</th><thstyle="text-align: left">female</th></tr><tr><tdstyle="text-align: left">1</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">no</td></tr><tr><tdstyle="text-align: left">2</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td></tr><tr><tdstyle="text-align: left">3</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td></tr><tr><tdstyle="text-align: left">4</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">yes</td></tr><tr><tdstyle="text-align: left">5</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td></tr><tr><tdstyle="text-align: left">6</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">yes</td></tr><tr><tdstyle="text-align: left">7</td><tdstyle="text-align: left">yes</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">no</td><tdstyle="text-align: left">no</td></tr></table><p>Even though each animal is said to participate in the mean, the result for the mean will now actually be the average of the base population. Math is weird sometimes.</p><p>Double-checking, the rank of <span>$X$</span> is still 4, so we can solve for the average of the base population, and the effect of being born in 1990, the effect of being born in 1991, and the effect of being male.</p><p>Whew! That was some transformation. We still haven't constrained this model enough to solve it, though. Now on to the genotype.</p><h2id="The-statistical-model:-genotype-as-random-effect"><aclass="docs-heading-anchor"href="#The-statistical-model:-genotype-as-random-effect">The statistical model: genotype as random effect</a><aid="The-statistical-model:-genotype-as-random-effect-1"></a><aclass="docs-heading-anchor-permalink"href="#The-statistical-model:-genotype-as-random-effect"title="Permalink"></a></h2><p>Remember I said above that genotype was a <strong>random effect</strong>? Statisticians say "<em>a random effect is an effect that influences the variance and not the mean of the observation in question.</em>" I'm not sure exactly what that means or how that is applicable to genotype, but it does let us add an additional constraint to our model.</p><p>The basic gist of genetics is that organisms that are related to one another are similar to one another. Based on a pedigree, we can even say how related to one another animals are, and quantify that as the amount that the genotype terms should be allowed to vary between related animals.</p><p>We'll need a pedigree for our animals:</p><h3id="Calf-Records-3"><aclass="docs-heading-anchor"href="#Calf-Records-3">Calf Records</a><aclass="docs-heading-anchor-permalink"href="#Calf-Records-3"title="Permalink"></a></h3><table><tr><thstyle="text-align: left">ID</th><thstyle="text-align: left">Sire</th><thstyle="text-align: left">Dam</th><thstyle="text-align: left">Birth Year</th><thstyle="text-align: left">Sex</th><thstyle="text-align: left">YW (kg)</th></tr><tr><tdstyle="text-align: left">1</td><tdstyle="text-align: left">NA</td><tdstyle="text-align: left">NA</td><tdstyle="text-align: left">1990</td><tdstyle="text-align: left">Male</td><tdstyle="text-align: left">354</td></tr><tr><tdstyle="text-align: left">2</td><tdstyle="text-align: left">NA</td><tdstyle="text-align: left">NA</td><tdstyle="text-align: left">1990</td><tdstyle="text-align: left">Female</t
\end{bmatrix}\]</p><p>Now I'm going to take the matrix with all of the <span>$u$</span> values, and call it <span>$μ$</span>. To quantify the idea of genetic relationship, I will then say that</p><pclass="math-container">\[\textup{var}(μ) = A σ_μ^2\]</p><p>Where:</p><ul><li><span>$A$</span> = the relationship matrix defined above</li><li><span>$σ_μ^2$</span> = the standard deviation of all the genotypes</li></ul><p>To fully constrain the system, I have to make two more assumptions: 1) that the error term in each animal's equation is independent from all other error terms, and 2) that the error term for each animal is independent from the value of the genotype. I will call the matrix holding the <span>$e$</span> values <span>$ε$</span> and then say</p><pclass="math-container">\[\textup{var}(ϵ) = I σ_ϵ^2\]</p><pclass="math-container">\[\textup{cov}(μ, ϵ) = \textup{cov}(ϵ, μ) = 0\]</p><p>Substituting in the matrix names, our equation now looks like</p><pclass="math-container">\[\begin{bmatrix}
+ ϵ\]</p><p>We are going to make three changes to this equation before we are ready to solve it, but they are cosmetic details for this example.</p><ol><li>Call the matrix on the left side of the equation <span>$Y$</span> (sometimes it's called the <strong>matrix of observations</strong>)</li><li>Multiply <span>$μ$</span> by an identity matrix called <span>$Z$</span>. Multiplying by the identity matrix is the matrix form of multiplying by one, so nothing changes, but if we later want to find one animal's genetic effect on another animal's performance (e.g. a <strong>maternal effects model</strong>), we can alter <span>$Z$</span> to allow that</li><li>Call the matrix with all the <span>$b$</span> values <span>$β$</span>.</li></ol><p>With all these changes, we now have</p><pclass="math-container">\[Y = Z μ + X β + ϵ\]</p><p>This is the canonical form of the mixed-model equation, and the form that Charles Henderson used to first predict breeding values of livestock.</p><h2id="Solving-the-equations"><aclass="docs-heading-anchor"href="#Solving-the-equations">Solving the equations</a><aid="Solving-the-equations-1"></a><aclass="docs-heading-anchor-permalink"href="#Solving-the-equations"title="Permalink"></a></h2><p>Henderson proved that the mixed-model equation can be solved by the following:</p><pclass="math-container">\[\begin{bmatrix}
\end{bmatrix}\]</p><p>Where</p><ul><li>The variables with hats are the statistical estimates of their mixed-model counterparts<ul><li>The predicted value of <span>$μ$</span> is called the <em>Best Linear Unbiased Predictor</em> or <em>BLUP</em></li><li>The estimated value of <span>$β$</span> is called the <em>Best Linear Unbiased Estimate</em> or <em>BLUE</em></li></ul></li><li>' is the transpose operator</li><li><span>$λ$</span> is a single real number that is a function of the heritability for the trait being predicted. It can be left out in many cases (<span>$λ = 1$</span>).<ul><li><span>$λ = \frac{1-h^2}{h^2}$</span></li></ul></li></ul><p>What happened to</p><h2id="Footnotes"><aclass="docs-heading-anchor"href="#Footnotes">Footnotes</a><aid="Footnotes-1"></a><aclass="docs-heading-anchor-permalink"href="#Footnotes"title="Permalink"></a></h2><h3id="Exception"><aclass="docs-heading-anchor"href="#Exception">Exception</a><aid="Exception-1"></a><aclass="docs-heading-anchor-permalink"href="#Exception"title="Permalink"></a></h3><p>An animal <strong>can</strong> share its genome with itself by a factor of more than one: that's called inbreeding! We can account for this, and <code>beefblup</code> does as it calculates <span>$A$</span>. This is an area that actually merits a good deal of study: see chapter 2 of <em>Linear Models for the Prediction of Animal Breeding Values</em> by Raphael A. Mrode (ISBN 978 1 78064 391 5).</p></article><navclass="docs-footer"><aclass="docs-footer-prevpage"href="../">« Home</a><aclass="docs-footer-nextpage"href="../beefblup-cli/">CLI Reference (WIP) »</a><divclass="flexbox-break"></div><pclass="footer-message">Powered by <ahref="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <ahref="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><divclass="modal"id="documenter-settings"><divclass="modal-background"></div><divclass="modal-card"><headerclass="modal-card-head"><pclass="modal-card-title">Settings</p><buttonclass="delete"></button></header><sectionclass="modal-card-body"><p><labelclass="label">Theme</label><divclass="select"><selectid="documenter-themepicker"><optionvalue="documenter-light">documenter-light</option><optionvalue="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <ahref="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 0.27.5 on <spanclass="colophon-date"title="Wednesday 1 September 2021 01:08">Wednesday 1 September 2021</span>. Using Julia version 1.5.4.</p></section><footerclass="modal-card-foot"></footer></div></div></div></body></html>
