{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lista 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rozważmy model: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{y} = \\mathbf{X\\beta+Za+\\epsilon}$, gdzie"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{y}$ - wektor z wartościami fenotypowmi ($n\\times 1$ w naszym przypadku n = ilość zwierząt) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{\\beta}$ - wektor efektów stałych ($p\\times 1$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{a}$ - wektor efektów losowych ($q\\times 1$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{\\epsilon}$ - wektor błędów ($n\\times 1$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{X}$ - macierz wystąpien dla efektów stałych ($n\\times p$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{Z}$ - macierz wystąpien dla efektów losowych ($n\\times q$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zakładamy, że macierze $\\mathbf{X}$ i $\\mathbf{Z}$ są niezależne. Dodatkowo zakładamy, że:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$E(\\mathbf{y}) = \\mathbf{X\\beta}$, $E(\\mathbf{a})=E(\\mathbf{\\epsilon})=0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$V=Var(\\mathbf{y})=Var(\\mathbf{Za})+ Var(\\mathbf{\\epsilon})=\\mathbf{Z}Var(\\mathbf{a})\\mathbf{Z}^{T}+\\mathbf{I}Var(\\mathbf{\\epsilon})\\mathbf{I}^{T}=\\mathbf{ZGZ}^{T}+\\mathbf{R}$, gdzie"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{G}=\\mathbf{A}\\cdot\\sigma^{2}_{a}$, $\\mathbf{R}=\\mathbf{I}\\cdot\\sigma^{2}_{\\epsilon}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ \\left[ \\begin{array}{cc}\n",
    "         \\mathbf{X}^{T}\\mathbf{X} & \\mathbf{X}^{T}\\mathbf{Z} \\\\\n",
    "         \\mathbf{Z}^{T}\\mathbf{X} & \\mathbf{Z}^{T}\\mathbf{Z}+\\mathbf{A}^{-1}\\alpha \\end{array}\\right]\\cdot\n",
    "  \\left[ \\begin{array}{c}\n",
    "         \\widehat{\\mathbf{\\beta}} \\\\\n",
    "         \\widehat{\\mathbf{a}} \\end{array}\\right]=\n",
    "  \\left[ \\begin{array}{c}\n",
    "         \\mathbf{X}^{T}\\mathbf{y} \\\\\n",
    "         \\mathbf{Z}^{T}\\mathbf{y} \\end{array}\\right]$, gdzie"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\alpha=\\frac{\\sigma^{2}_{\\epsilon}}{\\sigma^{2}_{a}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ \\left[ \\begin{array}{c}\n",
    "         \\widehat{\\mathbf{\\beta}} \\\\\n",
    "         \\widehat{\\mathbf{a}} \\end{array}\\right]=\n",
    "  \\left[ \\begin{array}{cc}\n",
    "         \\mathbf{X}^{T}\\mathbf{X} & \\mathbf{X}^{T}\\mathbf{Z} \\\\\n",
    "         \\mathbf{Z}^{T}\\mathbf{X} & \\mathbf{Z}^{T}\\mathbf{Z}+\\mathbf{A}^{-1}\\alpha \\end{array}\\right]^{-1}\\cdot\n",
    "  \\left[ \\begin{array}{c}\n",
    "         \\mathbf{X}^{T}\\mathbf{y} \\\\\n",
    "         \\mathbf{Z}^{T}\\mathbf{y} \\end{array}\\right]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ C = \\left[ \\begin{array}{cc}\n",
    "         \\mathbf{X}^{T}\\mathbf{X} & \\mathbf{X}^{T}\\mathbf{Z} \\\\\n",
    "         \\mathbf{Z}^{T}\\mathbf{X} & \\mathbf{Z}^{T}\\mathbf{Z}+\\mathbf{A}^{-1}\\alpha \\end{array}\\right] = \n",
    "         \\left[ \\begin{array}{cc}\n",
    "         \\mathbf{C}_{11} & \\mathbf{C}_{12} \\\\ \\mathbf{C}_{21} & \\mathbf{C}_{22} \\end{array}\\right]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ C^{-1} = \\left[ \\begin{array}{cc}\n",
    "         \\mathbf{C}^{11} & \\mathbf{C}^{12} \\\\ \\mathbf{C}^{21} & \\mathbf{C}^{22} \\end{array}\\right]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dokładność oszacowania wynosi: $r^2 = diag(1-\\mathbf{C}^{22}\\cdot\\alpha$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Liczymy macierz A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>1</th><th scope=col>2</th><th scope=col>3</th><th scope=col>4</th><th scope=col>5</th><th scope=col>6</th><th scope=col>7</th><th scope=col>8</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>1.00 </td><td>0.00 </td><td>0.00 </td><td>0.500</td><td>0.000</td><td>0.50 </td><td>0.25 </td><td>0.250</td></tr>\n",
       "\t<tr><td>0.00 </td><td>1.00 </td><td>0.00 </td><td>0.000</td><td>0.500</td><td>0.50 </td><td>0.25 </td><td>0.250</td></tr>\n",
       "\t<tr><td>0.00 </td><td>0.00 </td><td>1.00 </td><td>0.000</td><td>0.500</td><td>0.00 </td><td>0.25 </td><td>0.500</td></tr>\n",
       "\t<tr><td>0.50 </td><td>0.00 </td><td>0.00 </td><td>1.000</td><td>0.000</td><td>0.25 </td><td>0.50 </td><td>0.125</td></tr>\n",
       "\t<tr><td>0.00 </td><td>0.50 </td><td>0.50 </td><td>0.000</td><td>1.000</td><td>0.25 </td><td>0.50 </td><td>0.375</td></tr>\n",
       "\t<tr><td>0.50 </td><td>0.50 </td><td>0.00 </td><td>0.250</td><td>0.250</td><td>1.00 </td><td>0.25 </td><td>0.500</td></tr>\n",
       "\t<tr><td>0.25 </td><td>0.25 </td><td>0.25 </td><td>0.500</td><td>0.500</td><td>0.25 </td><td>1.00 </td><td>0.250</td></tr>\n",
       "\t<tr><td>0.25 </td><td>0.25 </td><td>0.50 </td><td>0.125</td><td>0.375</td><td>0.50 </td><td>0.25 </td><td>1.000</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|llllllll}\n",
       " 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8\\\\\n",
       "\\hline\n",
       "\t 1.00  & 0.00  & 0.00  & 0.500 & 0.000 & 0.50  & 0.25  & 0.250\\\\\n",
       "\t 0.00  & 1.00  & 0.00  & 0.000 & 0.500 & 0.50  & 0.25  & 0.250\\\\\n",
       "\t 0.00  & 0.00  & 1.00  & 0.000 & 0.500 & 0.00  & 0.25  & 0.500\\\\\n",
       "\t 0.50  & 0.00  & 0.00  & 1.000 & 0.000 & 0.25  & 0.50  & 0.125\\\\\n",
       "\t 0.00  & 0.50  & 0.50  & 0.000 & 1.000 & 0.25  & 0.50  & 0.375\\\\\n",
       "\t 0.50  & 0.50  & 0.00  & 0.250 & 0.250 & 1.00  & 0.25  & 0.500\\\\\n",
       "\t 0.25  & 0.25  & 0.25  & 0.500 & 0.500 & 0.25  & 1.00  & 0.250\\\\\n",
       "\t 0.25  & 0.25  & 0.50  & 0.125 & 0.375 & 0.50  & 0.25  & 1.000\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |\n",
       "|---|---|---|---|---|---|---|---|\n",
       "| 1.00  | 0.00  | 0.00  | 0.500 | 0.000 | 0.50  | 0.25  | 0.250 |\n",
       "| 0.00  | 1.00  | 0.00  | 0.000 | 0.500 | 0.50  | 0.25  | 0.250 |\n",
       "| 0.00  | 0.00  | 1.00  | 0.000 | 0.500 | 0.00  | 0.25  | 0.500 |\n",
       "| 0.50  | 0.00  | 0.00  | 1.000 | 0.000 | 0.25  | 0.50  | 0.125 |\n",
       "| 0.00  | 0.50  | 0.50  | 0.000 | 1.000 | 0.25  | 0.50  | 0.375 |\n",
       "| 0.50  | 0.50  | 0.00  | 0.250 | 0.250 | 1.00  | 0.25  | 0.500 |\n",
       "| 0.25  | 0.25  | 0.25  | 0.500 | 0.500 | 0.25  | 1.00  | 0.250 |\n",
       "| 0.25  | 0.25  | 0.50  | 0.125 | 0.375 | 0.50  | 0.25  | 1.000 |\n",
       "\n"
      ],
      "text/plain": [
       "  1    2    3    4     5     6    7    8    \n",
       "1 1.00 0.00 0.00 0.500 0.000 0.50 0.25 0.250\n",
       "2 0.00 1.00 0.00 0.000 0.500 0.50 0.25 0.250\n",
       "3 0.00 0.00 1.00 0.000 0.500 0.00 0.25 0.500\n",
       "4 0.50 0.00 0.00 1.000 0.000 0.25 0.50 0.125\n",
       "5 0.00 0.50 0.50 0.000 1.000 0.25 0.50 0.375\n",
       "6 0.50 0.50 0.00 0.250 0.250 1.00 0.25 0.500\n",
       "7 0.25 0.25 0.25 0.500 0.500 0.25 1.00 0.250\n",
       "8 0.25 0.25 0.50 0.125 0.375 0.50 0.25 1.000"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "library(pedigreemm)\n",
    "\n",
    "id = 1:8\n",
    "sire = c(NA, NA, NA, 1, 3, 1, 4, 3)\n",
    "dam = c(NA, NA, NA, NA, 2, 2, 5, 6)\n",
    "\n",
    "ped = pedigree(sire = sire, dam = dam, label = id)\n",
    "(A = as.matrix(getA(ped)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Definiujemy kolejne źródła informacji"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<tbody>\n",
       "\t<tr><td>4.5</td><td>2.9</td><td>3.9</td><td>3.5</td><td>5  </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{lllll}\n",
       "\t 4.5 & 2.9 & 3.9 & 3.5 & 5  \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| 4.5 | 2.9 | 3.9 | 3.5 | 5   |\n",
       "\n"
      ],
      "text/plain": [
       "     [,1] [,2] [,3] [,4] [,5]\n",
       "[1,] 4.5  2.9  3.9  3.5  5   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<tbody>\n",
       "\t<tr><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr>\n",
       "\t<tr><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{lllll}\n",
       "\t 1 & 0 & 0 & 1 & 1\\\\\n",
       "\t 0 & 1 & 1 & 0 & 0\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| 1 | 0 | 0 | 1 | 1 |\n",
       "| 0 | 1 | 1 | 0 | 0 |\n",
       "\n"
      ],
      "text/plain": [
       "     [,1] [,2] [,3] [,4] [,5]\n",
       "[1,] 1    0    0    1    1   \n",
       "[2,] 0    1    1    0    0   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<tbody>\n",
       "\t<tr><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td></tr>\n",
       "\t<tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{llllllll}\n",
       "\t 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\\\\n",
       "\t 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\\\\n",
       "\t 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\\\\n",
       "\t 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\\\\n",
       "\t 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 |\n",
       "| 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |\n",
       "| 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 |\n",
       "| 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 |\n",
       "| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 |\n",
       "\n"
      ],
      "text/plain": [
       "     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]\n",
       "[1,] 0    0    0    1    0    0    0    0   \n",
       "[2,] 0    0    0    0    1    0    0    0   \n",
       "[3,] 0    0    0    0    0    1    0    0   \n",
       "[4,] 0    0    0    0    0    0    1    0   \n",
       "[5,] 0    0    0    0    0    0    0    1   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y = as.matrix(c(4.5, 2.9, 3.9, 3.5, 5.0))\n",
    "t(y)\n",
    "\n",
    "sex = c(1, 0, 0, 1, 1)\n",
    "X = matrix(0, 5, 2)\n",
    "X[,1] = sex\n",
    "X[,2] = 1-sex\n",
    "t(X)\n",
    "\n",
    "I = diag(5)\n",
    "Z = matrix(0, 5, 8)\n",
    "Z[1:5, 4:8] = I\n",
    "Z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Rozwiązujemy układ równan mieszanych"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<dl>\n",
       "\t<dt>$mC</dt>\n",
       "\t\t<dd><table>\n",
       "<tbody>\n",
       "\t<tr><td>3            </td><td>0            </td><td> 0.000000e+00</td><td> 0.000000e+00</td><td> 0.000000e+00</td><td> 1.000000e+00</td><td> 0.000000e+00</td><td> 0.000000e+00</td><td> 1.000000e+00</td><td> 1.000000e+00</td></tr>\n",
       "\t<tr><td>0            </td><td>2            </td><td> 0.000000e+00</td><td> 0.000000e+00</td><td> 0.000000e+00</td><td> 0.000000e+00</td><td> 1.000000e+00</td><td> 1.000000e+00</td><td> 0.000000e+00</td><td> 0.000000e+00</td></tr>\n",
       "\t<tr><td>0            </td><td>0            </td><td> 3.666667e+00</td><td> 1.000000e+00</td><td> 5.107026e-15</td><td>-1.333333e+00</td><td> 1.110223e-15</td><td>-2.000000e+00</td><td>-9.547918e-15</td><td>-4.309053e-15</td></tr>\n",
       "\t<tr><td>0            </td><td>0            </td><td> 1.000000e+00</td><td> 4.000000e+00</td><td> 1.000000e+00</td><td>-2.664535e-15</td><td>-2.000000e+00</td><td>-2.000000e+00</td><td> 1.776357e-15</td><td> 2.914335e-16</td></tr>\n",
       "\t<tr><td>0            </td><td>0            </td><td> 1.665335e-15</td><td> 1.000000e+00</td><td> 4.000000e+00</td><td> 9.992007e-16</td><td>-2.000000e+00</td><td> 1.000000e+00</td><td>-7.771561e-16</td><td>-2.000000e+00</td></tr>\n",
       "\t<tr><td>1            </td><td>0            </td><td>-1.333333e+00</td><td>-2.664535e-15</td><td>-4.329870e-15</td><td> 4.666667e+00</td><td> 1.000000e+00</td><td> 1.998401e-15</td><td>-2.000000e+00</td><td> 1.915135e-15</td></tr>\n",
       "\t<tr><td>0            </td><td>1            </td><td>-2.220446e-16</td><td>-2.000000e+00</td><td>-2.000000e+00</td><td> 1.000000e+00</td><td> 6.000000e+00</td><td> 8.881784e-16</td><td>-2.000000e+00</td><td>-4.232725e-15</td></tr>\n",
       "\t<tr><td>0            </td><td>1            </td><td>-2.000000e+00</td><td>-2.000000e+00</td><td> 1.000000e+00</td><td>-5.551115e-16</td><td>-1.776357e-15</td><td> 6.000000e+00</td><td> 1.221245e-15</td><td>-2.000000e+00</td></tr>\n",
       "\t<tr><td>1            </td><td>0            </td><td> 4.440892e-16</td><td>-1.332268e-15</td><td>-4.440892e-16</td><td>-2.000000e+00</td><td>-2.000000e+00</td><td>-7.771561e-16</td><td> 5.000000e+00</td><td> 1.172673e-15</td></tr>\n",
       "\t<tr><td>1            </td><td>0            </td><td>-1.998401e-15</td><td>-7.709111e-15</td><td>-2.000000e+00</td><td>-5.509482e-15</td><td> 1.970646e-15</td><td>-2.000000e+00</td><td> 5.606626e-15</td><td> 5.000000e+00</td></tr>\n",
       "</tbody>\n",
       "</table>\n",
       "</dd>\n",
       "\t<dt>$est</dt>\n",
       "\t\t<dd><table>\n",
       "<tbody>\n",
       "\t<tr><td> 4.358502330</td></tr>\n",
       "\t<tr><td> 3.404430006</td></tr>\n",
       "\t<tr><td> 0.098444576</td></tr>\n",
       "\t<tr><td>-0.018770099</td></tr>\n",
       "\t<tr><td>-0.041084203</td></tr>\n",
       "\t<tr><td>-0.008663123</td></tr>\n",
       "\t<tr><td>-0.185732099</td></tr>\n",
       "\t<tr><td> 0.176872088</td></tr>\n",
       "\t<tr><td>-0.249458555</td></tr>\n",
       "\t<tr><td> 0.182614688</td></tr>\n",
       "</tbody>\n",
       "</table>\n",
       "</dd>\n",
       "</dl>\n"
      ],
      "text/latex": [
       "\\begin{description}\n",
       "\\item[\\$mC] \\begin{tabular}{llllllllll}\n",
       "\t 3             & 0             &  0.000000e+00 &  0.000000e+00 &  0.000000e+00 &  1.000000e+00 &  0.000000e+00 &  0.000000e+00 &  1.000000e+00 &  1.000000e+00\\\\\n",
       "\t 0             & 2             &  0.000000e+00 &  0.000000e+00 &  0.000000e+00 &  0.000000e+00 &  1.000000e+00 &  1.000000e+00 &  0.000000e+00 &  0.000000e+00\\\\\n",
       "\t 0             & 0             &  3.666667e+00 &  1.000000e+00 &  5.107026e-15 & -1.333333e+00 &  1.110223e-15 & -2.000000e+00 & -9.547918e-15 & -4.309053e-15\\\\\n",
       "\t 0             & 0             &  1.000000e+00 &  4.000000e+00 &  1.000000e+00 & -2.664535e-15 & -2.000000e+00 & -2.000000e+00 &  1.776357e-15 &  2.914335e-16\\\\\n",
       "\t 0             & 0             &  1.665335e-15 &  1.000000e+00 &  4.000000e+00 &  9.992007e-16 & -2.000000e+00 &  1.000000e+00 & -7.771561e-16 & -2.000000e+00\\\\\n",
       "\t 1             & 0             & -1.333333e+00 & -2.664535e-15 & -4.329870e-15 &  4.666667e+00 &  1.000000e+00 &  1.998401e-15 & -2.000000e+00 &  1.915135e-15\\\\\n",
       "\t 0             & 1             & -2.220446e-16 & -2.000000e+00 & -2.000000e+00 &  1.000000e+00 &  6.000000e+00 &  8.881784e-16 & -2.000000e+00 & -4.232725e-15\\\\\n",
       "\t 0             & 1             & -2.000000e+00 & -2.000000e+00 &  1.000000e+00 & -5.551115e-16 & -1.776357e-15 &  6.000000e+00 &  1.221245e-15 & -2.000000e+00\\\\\n",
       "\t 1             & 0             &  4.440892e-16 & -1.332268e-15 & -4.440892e-16 & -2.000000e+00 & -2.000000e+00 & -7.771561e-16 &  5.000000e+00 &  1.172673e-15\\\\\n",
       "\t 1             & 0             & -1.998401e-15 & -7.709111e-15 & -2.000000e+00 & -5.509482e-15 &  1.970646e-15 & -2.000000e+00 &  5.606626e-15 &  5.000000e+00\\\\\n",
       "\\end{tabular}\n",
       "\n",
       "\\item[\\$est] \\begin{tabular}{l}\n",
       "\t  4.358502330\\\\\n",
       "\t  3.404430006\\\\\n",
       "\t  0.098444576\\\\\n",
       "\t -0.018770099\\\\\n",
       "\t -0.041084203\\\\\n",
       "\t -0.008663123\\\\\n",
       "\t -0.185732099\\\\\n",
       "\t  0.176872088\\\\\n",
       "\t -0.249458555\\\\\n",
       "\t  0.182614688\\\\\n",
       "\\end{tabular}\n",
       "\n",
       "\\end{description}\n"
      ],
      "text/markdown": [
       "$mC\n",
       ":   \n",
       "| 3             | 0             |  0.000000e+00 |  0.000000e+00 |  0.000000e+00 |  1.000000e+00 |  0.000000e+00 |  0.000000e+00 |  1.000000e+00 |  1.000000e+00 |\n",
       "| 0             | 2             |  0.000000e+00 |  0.000000e+00 |  0.000000e+00 |  0.000000e+00 |  1.000000e+00 |  1.000000e+00 |  0.000000e+00 |  0.000000e+00 |\n",
       "| 0             | 0             |  3.666667e+00 |  1.000000e+00 |  5.107026e-15 | -1.333333e+00 |  1.110223e-15 | -2.000000e+00 | -9.547918e-15 | -4.309053e-15 |\n",
       "| 0             | 0             |  1.000000e+00 |  4.000000e+00 |  1.000000e+00 | -2.664535e-15 | -2.000000e+00 | -2.000000e+00 |  1.776357e-15 |  2.914335e-16 |\n",
       "| 0             | 0             |  1.665335e-15 |  1.000000e+00 |  4.000000e+00 |  9.992007e-16 | -2.000000e+00 |  1.000000e+00 | -7.771561e-16 | -2.000000e+00 |\n",
       "| 1             | 0             | -1.333333e+00 | -2.664535e-15 | -4.329870e-15 |  4.666667e+00 |  1.000000e+00 |  1.998401e-15 | -2.000000e+00 |  1.915135e-15 |\n",
       "| 0             | 1             | -2.220446e-16 | -2.000000e+00 | -2.000000e+00 |  1.000000e+00 |  6.000000e+00 |  8.881784e-16 | -2.000000e+00 | -4.232725e-15 |\n",
       "| 0             | 1             | -2.000000e+00 | -2.000000e+00 |  1.000000e+00 | -5.551115e-16 | -1.776357e-15 |  6.000000e+00 |  1.221245e-15 | -2.000000e+00 |\n",
       "| 1             | 0             |  4.440892e-16 | -1.332268e-15 | -4.440892e-16 | -2.000000e+00 | -2.000000e+00 | -7.771561e-16 |  5.000000e+00 |  1.172673e-15 |\n",
       "| 1             | 0             | -1.998401e-15 | -7.709111e-15 | -2.000000e+00 | -5.509482e-15 |  1.970646e-15 | -2.000000e+00 |  5.606626e-15 |  5.000000e+00 |\n",
       "\n",
       "\n",
       "$est\n",
       ":   \n",
       "|  4.358502330 |\n",
       "|  3.404430006 |\n",
       "|  0.098444576 |\n",
       "| -0.018770099 |\n",
       "| -0.041084203 |\n",
       "| -0.008663123 |\n",
       "| -0.185732099 |\n",
       "|  0.176872088 |\n",
       "| -0.249458555 |\n",
       "|  0.182614688 |\n",
       "\n",
       "\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "$mC\n",
       "      [,1] [,2]          [,3]          [,4]          [,5]          [,6]\n",
       " [1,]    3    0  0.000000e+00  0.000000e+00  0.000000e+00  1.000000e+00\n",
       " [2,]    0    2  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00\n",
       " [3,]    0    0  3.666667e+00  1.000000e+00  5.107026e-15 -1.333333e+00\n",
       " [4,]    0    0  1.000000e+00  4.000000e+00  1.000000e+00 -2.664535e-15\n",
       " [5,]    0    0  1.665335e-15  1.000000e+00  4.000000e+00  9.992007e-16\n",
       " [6,]    1    0 -1.333333e+00 -2.664535e-15 -4.329870e-15  4.666667e+00\n",
       " [7,]    0    1 -2.220446e-16 -2.000000e+00 -2.000000e+00  1.000000e+00\n",
       " [8,]    0    1 -2.000000e+00 -2.000000e+00  1.000000e+00 -5.551115e-16\n",
       " [9,]    1    0  4.440892e-16 -1.332268e-15 -4.440892e-16 -2.000000e+00\n",
       "[10,]    1    0 -1.998401e-15 -7.709111e-15 -2.000000e+00 -5.509482e-15\n",
       "               [,7]          [,8]          [,9]         [,10]\n",
       " [1,]  0.000000e+00  0.000000e+00  1.000000e+00  1.000000e+00\n",
       " [2,]  1.000000e+00  1.000000e+00  0.000000e+00  0.000000e+00\n",
       " [3,]  1.110223e-15 -2.000000e+00 -9.547918e-15 -4.309053e-15\n",
       " [4,] -2.000000e+00 -2.000000e+00  1.776357e-15  2.914335e-16\n",
       " [5,] -2.000000e+00  1.000000e+00 -7.771561e-16 -2.000000e+00\n",
       " [6,]  1.000000e+00  1.998401e-15 -2.000000e+00  1.915135e-15\n",
       " [7,]  6.000000e+00  8.881784e-16 -2.000000e+00 -4.232725e-15\n",
       " [8,] -1.776357e-15  6.000000e+00  1.221245e-15 -2.000000e+00\n",
       " [9,] -2.000000e+00 -7.771561e-16  5.000000e+00  1.172673e-15\n",
       "[10,]  1.970646e-15 -2.000000e+00  5.606626e-15  5.000000e+00\n",
       "\n",
       "$est\n",
       "              [,1]\n",
       " [1,]  4.358502330\n",
       " [2,]  3.404430006\n",
       " [3,]  0.098444576\n",
       " [4,] -0.018770099\n",
       " [5,] -0.041084203\n",
       " [6,] -0.008663123\n",
       " [7,] -0.185732099\n",
       " [8,]  0.176872088\n",
       " [9,] -0.249458555\n",
       "[10,]  0.182614688\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "library(MASS)\n",
    "\n",
    "mme = function(y, X, Z, A, sigma_a, sigma_e) {\n",
    "    alpha = sigma_e / sigma_a\n",
    "    invA = ginv(A)\n",
    "    C = rbind(cbind(t(X)%*%X, t(X)%*%Z),\n",
    "              cbind(t(Z)%*%X, t(Z)%*%Z+invA*alpha))\n",
    "    rhs = rbind(t(X)%*%y, t(Z)%*%y)\n",
    "    invC = ginv(C)\n",
    "    estimators = invC%*%rhs\n",
    "    list(mC = C, est = estimators)\n",
    "}\n",
    "\n",
    "mme(y, X, Z, A, 20, 40)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Dokładność wartości hodowlanych"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<tbody>\n",
       "\t<tr><td> 0.59556097 </td><td> 0.15730213 </td><td>-0.164119413</td><td>-0.083624565</td><td>-0.13059083 </td><td>-0.26455749 </td><td>-0.14827804 </td><td>-0.16632621 </td><td>-0.2842464  </td><td>-0.2378790  </td></tr>\n",
       "\t<tr><td> 0.15730213 </td><td> 0.80245865 </td><td>-0.132863260</td><td>-0.241250738</td><td>-0.11196840 </td><td>-0.08730803 </td><td>-0.29891465 </td><td>-0.30600266 </td><td>-0.1859495  </td><td>-0.1986488  </td></tr>\n",
       "\t<tr><td>-0.16411941 </td><td>-0.13286326 </td><td> 0.471094211</td><td> 0.006928037</td><td> 0.03264668 </td><td> 0.21954371 </td><td> 0.04495225 </td><td> 0.22077427 </td><td> 0.1386223  </td><td> 0.1341923  </td></tr>\n",
       "\t<tr><td>-0.08362457 </td><td>-0.24125074 </td><td> 0.006928037</td><td> 0.492095721</td><td>-0.01030797 </td><td> 0.02039033 </td><td> 0.23734577 </td><td> 0.24515571 </td><td> 0.1198194  </td><td> 0.1106640  </td></tr>\n",
       "\t<tr><td>-0.13059083 </td><td>-0.11196840 </td><td> 0.032646682</td><td>-0.010307967</td><td> 0.45645878 </td><td> 0.04812709 </td><td> 0.20132326 </td><td> 0.02261354 </td><td> 0.1258983  </td><td> 0.2177471  </td></tr>\n",
       "\t<tr><td>-0.26455749 </td><td>-0.08730803 </td><td> 0.219543709</td><td> 0.020390333</td><td> 0.04812709 </td><td> 0.42768015 </td><td> 0.04704420 </td><td> 0.12757186 </td><td> 0.2428012  </td><td> 0.1231911  </td></tr>\n",
       "\t<tr><td>-0.14827804 </td><td>-0.29891465 </td><td> 0.044952254</td><td> 0.237345770</td><td> 0.20132326 </td><td> 0.04704420 </td><td> 0.42810675 </td><td> 0.16972255 </td><td> 0.2197160  </td><td> 0.1780739  </td></tr>\n",
       "\t<tr><td>-0.16632621 </td><td>-0.30600266 </td><td> 0.220774267</td><td> 0.245155707</td><td> 0.02261354 </td><td> 0.12757186 </td><td> 0.16972255 </td><td> 0.44228277 </td><td> 0.1521830  </td><td> 0.2192238  </td></tr>\n",
       "\t<tr><td>-0.28424641 </td><td>-0.18594950 </td><td> 0.138622268</td><td> 0.119819354</td><td> 0.12589831 </td><td> 0.24280124 </td><td> 0.21971599 </td><td> 0.15218301 </td><td> 0.4418562  </td><td> 0.1680818  </td></tr>\n",
       "\t<tr><td>-0.23787901 </td><td>-0.19864885 </td><td> 0.134192262</td><td> 0.110664009</td><td> 0.21774710 </td><td> 0.12319108 </td><td> 0.17807393 </td><td> 0.21922376 </td><td> 0.1680818  </td><td> 0.4223641  </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{llllllllll}\n",
       "\t  0.59556097  &  0.15730213  & -0.164119413 & -0.083624565 & -0.13059083  & -0.26455749  & -0.14827804  & -0.16632621  & -0.2842464   & -0.2378790  \\\\\n",
       "\t  0.15730213  &  0.80245865  & -0.132863260 & -0.241250738 & -0.11196840  & -0.08730803  & -0.29891465  & -0.30600266  & -0.1859495   & -0.1986488  \\\\\n",
       "\t -0.16411941  & -0.13286326  &  0.471094211 &  0.006928037 &  0.03264668  &  0.21954371  &  0.04495225  &  0.22077427  &  0.1386223   &  0.1341923  \\\\\n",
       "\t -0.08362457  & -0.24125074  &  0.006928037 &  0.492095721 & -0.01030797  &  0.02039033  &  0.23734577  &  0.24515571  &  0.1198194   &  0.1106640  \\\\\n",
       "\t -0.13059083  & -0.11196840  &  0.032646682 & -0.010307967 &  0.45645878  &  0.04812709  &  0.20132326  &  0.02261354  &  0.1258983   &  0.2177471  \\\\\n",
       "\t -0.26455749  & -0.08730803  &  0.219543709 &  0.020390333 &  0.04812709  &  0.42768015  &  0.04704420  &  0.12757186  &  0.2428012   &  0.1231911  \\\\\n",
       "\t -0.14827804  & -0.29891465  &  0.044952254 &  0.237345770 &  0.20132326  &  0.04704420  &  0.42810675  &  0.16972255  &  0.2197160   &  0.1780739  \\\\\n",
       "\t -0.16632621  & -0.30600266  &  0.220774267 &  0.245155707 &  0.02261354  &  0.12757186  &  0.16972255  &  0.44228277  &  0.1521830   &  0.2192238  \\\\\n",
       "\t -0.28424641  & -0.18594950  &  0.138622268 &  0.119819354 &  0.12589831  &  0.24280124  &  0.21971599  &  0.15218301  &  0.4418562   &  0.1680818  \\\\\n",
       "\t -0.23787901  & -0.19864885  &  0.134192262 &  0.110664009 &  0.21774710  &  0.12319108  &  0.17807393  &  0.21922376  &  0.1680818   &  0.4223641  \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "|  0.59556097  |  0.15730213  | -0.164119413 | -0.083624565 | -0.13059083  | -0.26455749  | -0.14827804  | -0.16632621  | -0.2842464   | -0.2378790   |\n",
       "|  0.15730213  |  0.80245865  | -0.132863260 | -0.241250738 | -0.11196840  | -0.08730803  | -0.29891465  | -0.30600266  | -0.1859495   | -0.1986488   |\n",
       "| -0.16411941  | -0.13286326  |  0.471094211 |  0.006928037 |  0.03264668  |  0.21954371  |  0.04495225  |  0.22077427  |  0.1386223   |  0.1341923   |\n",
       "| -0.08362457  | -0.24125074  |  0.006928037 |  0.492095721 | -0.01030797  |  0.02039033  |  0.23734577  |  0.24515571  |  0.1198194   |  0.1106640   |\n",
       "| -0.13059083  | -0.11196840  |  0.032646682 | -0.010307967 |  0.45645878  |  0.04812709  |  0.20132326  |  0.02261354  |  0.1258983   |  0.2177471   |\n",
       "| -0.26455749  | -0.08730803  |  0.219543709 |  0.020390333 |  0.04812709  |  0.42768015  |  0.04704420  |  0.12757186  |  0.2428012   |  0.1231911   |\n",
       "| -0.14827804  | -0.29891465  |  0.044952254 |  0.237345770 |  0.20132326  |  0.04704420  |  0.42810675  |  0.16972255  |  0.2197160   |  0.1780739   |\n",
       "| -0.16632621  | -0.30600266  |  0.220774267 |  0.245155707 |  0.02261354  |  0.12757186  |  0.16972255  |  0.44228277  |  0.1521830   |  0.2192238   |\n",
       "| -0.28424641  | -0.18594950  |  0.138622268 |  0.119819354 |  0.12589831  |  0.24280124  |  0.21971599  |  0.15218301  |  0.4418562   |  0.1680818   |\n",
       "| -0.23787901  | -0.19864885  |  0.134192262 |  0.110664009 |  0.21774710  |  0.12319108  |  0.17807393  |  0.21922376  |  0.1680818   |  0.4223641   |\n",
       "\n"
      ],
      "text/plain": [
       "      [,1]        [,2]        [,3]         [,4]         [,5]        [,6]       \n",
       " [1,]  0.59556097  0.15730213 -0.164119413 -0.083624565 -0.13059083 -0.26455749\n",
       " [2,]  0.15730213  0.80245865 -0.132863260 -0.241250738 -0.11196840 -0.08730803\n",
       " [3,] -0.16411941 -0.13286326  0.471094211  0.006928037  0.03264668  0.21954371\n",
       " [4,] -0.08362457 -0.24125074  0.006928037  0.492095721 -0.01030797  0.02039033\n",
       " [5,] -0.13059083 -0.11196840  0.032646682 -0.010307967  0.45645878  0.04812709\n",
       " [6,] -0.26455749 -0.08730803  0.219543709  0.020390333  0.04812709  0.42768015\n",
       " [7,] -0.14827804 -0.29891465  0.044952254  0.237345770  0.20132326  0.04704420\n",
       " [8,] -0.16632621 -0.30600266  0.220774267  0.245155707  0.02261354  0.12757186\n",
       " [9,] -0.28424641 -0.18594950  0.138622268  0.119819354  0.12589831  0.24280124\n",
       "[10,] -0.23787901 -0.19864885  0.134192262  0.110664009  0.21774710  0.12319108\n",
       "      [,7]        [,8]        [,9]       [,10]     \n",
       " [1,] -0.14827804 -0.16632621 -0.2842464 -0.2378790\n",
       " [2,] -0.29891465 -0.30600266 -0.1859495 -0.1986488\n",
       " [3,]  0.04495225  0.22077427  0.1386223  0.1341923\n",
       " [4,]  0.23734577  0.24515571  0.1198194  0.1106640\n",
       " [5,]  0.20132326  0.02261354  0.1258983  0.2177471\n",
       " [6,]  0.04704420  0.12757186  0.2428012  0.1231911\n",
       " [7,]  0.42810675  0.16972255  0.2197160  0.1780739\n",
       " [8,]  0.16972255  0.44228277  0.1521830  0.2192238\n",
       " [9,]  0.21971599  0.15218301  0.4418562  0.1680818\n",
       "[10,]  0.17807393  0.21922376  0.1680818  0.4223641"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>0.0578115770821027</li>\n",
       "\t<li>0.0158085581151157</li>\n",
       "\t<li>0.0870824309247226</li>\n",
       "\t<li>0.144639692852924</li>\n",
       "\t<li>0.143786506530158</li>\n",
       "\t<li>0.115434468727441</li>\n",
       "\t<li>0.116287655050208</li>\n",
       "\t<li>0.155271707028943</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0.0578115770821027\n",
       "\\item 0.0158085581151157\n",
       "\\item 0.0870824309247226\n",
       "\\item 0.144639692852924\n",
       "\\item 0.143786506530158\n",
       "\\item 0.115434468727441\n",
       "\\item 0.116287655050208\n",
       "\\item 0.155271707028943\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0.0578115770821027\n",
       "2. 0.0158085581151157\n",
       "3. 0.0870824309247226\n",
       "4. 0.144639692852924\n",
       "5. 0.143786506530158\n",
       "6. 0.115434468727441\n",
       "7. 0.116287655050208\n",
       "8. 0.155271707028943\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 0.05781158 0.01580856 0.08708243 0.14463969 0.14378651 0.11543447 0.11628766\n",
       "[8] 0.15527171"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>0.240440381554561</li>\n",
       "\t<li>0.125732088645324</li>\n",
       "\t<li>0.295097324496042</li>\n",
       "\t<li>0.38031525456248</li>\n",
       "\t<li>0.379191912532636</li>\n",
       "\t<li>0.339756484452381</li>\n",
       "\t<li>0.341009757998518</li>\n",
       "\t<li>0.39404531088308</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0.240440381554561\n",
       "\\item 0.125732088645324\n",
       "\\item 0.295097324496042\n",
       "\\item 0.38031525456248\n",
       "\\item 0.379191912532636\n",
       "\\item 0.339756484452381\n",
       "\\item 0.341009757998518\n",
       "\\item 0.39404531088308\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0.240440381554561\n",
       "2. 0.125732088645324\n",
       "3. 0.295097324496042\n",
       "4. 0.38031525456248\n",
       "5. 0.379191912532636\n",
       "6. 0.339756484452381\n",
       "7. 0.341009757998518\n",
       "8. 0.39404531088308\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 0.2404404 0.1257321 0.2950973 0.3803153 0.3791919 0.3397565 0.3410098\n",
       "[8] 0.3940453"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "C = as.matrix(mme(y, X, Z, A, 20, 40)$mC)\n",
    "(invC = ginv(C))\n",
    "\n",
    "invC22 = invC[3:10, 3:10]\n",
    "(r2 = diag(1 - invC22*2))\n",
    "(r = sqrt(r2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Odziedziczalność cechy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.333333333333333"
      ],
      "text/latex": [
       "0.333333333333333"
      ],
      "text/markdown": [
       "0.333333333333333"
      ],
      "text/plain": [
       "[1] 0.3333333"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "(h2 = 20 / (20 + 40))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6. Istotność efektów stałych"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test Walda:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$W = \\frac{\\widehat{\\beta}}{se(\\widehat{\\beta})}\\sim\\mathcal{N}(0,1)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$H_{0}: \\beta=0$ vs. $H_{1}: \\beta\\neq 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$var(\\widehat{\\beta}) = \\left(\\mathbf{X}^{T}\\mathbf{V}^{-1}\\mathbf{X}\\right)^{-1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<tbody>\n",
       "\t<tr><td>60.0</td><td> 0.0</td><td> 5  </td><td>10  </td><td> 2.5</td></tr>\n",
       "\t<tr><td> 0.0</td><td>60.0</td><td> 5  </td><td>10  </td><td> 7.5</td></tr>\n",
       "\t<tr><td> 5.0</td><td> 5.0</td><td>60  </td><td> 5  </td><td>10.0</td></tr>\n",
       "\t<tr><td>10.0</td><td>10.0</td><td> 5  </td><td>60  </td><td> 5.0</td></tr>\n",
       "\t<tr><td> 2.5</td><td> 7.5</td><td>10  </td><td> 5  </td><td>60.0</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{lllll}\n",
       "\t 60.0 &  0.0 &  5   & 10   &  2.5\\\\\n",
       "\t  0.0 & 60.0 &  5   & 10   &  7.5\\\\\n",
       "\t  5.0 &  5.0 & 60   &  5   & 10.0\\\\\n",
       "\t 10.0 & 10.0 &  5   & 60   &  5.0\\\\\n",
       "\t  2.5 &  7.5 & 10   &  5   & 60.0\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| 60.0 |  0.0 |  5   | 10   |  2.5 |\n",
       "|  0.0 | 60.0 |  5   | 10   |  7.5 |\n",
       "|  5.0 |  5.0 | 60   |  5   | 10.0 |\n",
       "| 10.0 | 10.0 |  5   | 60   |  5.0 |\n",
       "|  2.5 |  7.5 | 10   |  5   | 60.0 |\n",
       "\n"
      ],
      "text/plain": [
       "     [,1] [,2] [,3] [,4] [,5]\n",
       "[1,] 60.0  0.0  5   10    2.5\n",
       "[2,]  0.0 60.0  5   10    7.5\n",
       "[3,]  5.0  5.0 60    5   10.0\n",
       "[4,] 10.0 10.0  5   60    5.0\n",
       "[5,]  2.5  7.5 10    5   60.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<tbody>\n",
       "\t<tr><td>23.822439</td><td> 6.292085</td></tr>\n",
       "\t<tr><td> 6.292085</td><td>32.098346</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{ll}\n",
       "\t 23.822439 &  6.292085\\\\\n",
       "\t  6.292085 & 32.098346\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| 23.822439 |  6.292085 |\n",
       "|  6.292085 | 32.098346 |\n",
       "\n"
      ],
      "text/plain": [
       "     [,1]      [,2]     \n",
       "[1,] 23.822439  6.292085\n",
       "[2,]  6.292085 32.098346"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>4.88082357807458</li>\n",
       "\t<li>5.66554023294585</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 4.88082357807458\n",
       "\\item 5.66554023294585\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 4.88082357807458\n",
       "2. 5.66554023294585\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 4.880824 5.665540"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>0.892985017822407</li>\n",
       "\t<li>0.600901214346598</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0.892985017822407\n",
       "\\item 0.600901214346598\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0.892985017822407\n",
       "2. 0.600901214346598\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 0.8929850 0.6009012"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>0.371865195985295</li>\n",
       "\t<li>0.547905784351094</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0.371865195985295\n",
       "\\item 0.547905784351094\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0.371865195985295\n",
       "2. 0.547905784351094\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 0.3718652 0.5479058"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "G = A*20\n",
    "R = diag(5)*40\n",
    "V = Z%*%G%*%t(Z) + R\n",
    "V\n",
    "\n",
    "(varB = ginv(t(X)%*%ginv(V)%*%X))\n",
    "(seB = sqrt(diag(varB)))\n",
    "\n",
    "(testWalda = mme(y, X, Z, A, 20, 40)$est[1:2] / seB)\n",
    "(p_value = 2*pnorm(abs(testWalda), lower.tail = FALSE))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}