{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "#### Model mieszany + estymacja parametrów wariancji + genetyka" ] }, { "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$), $\\mathbf{a}\\sim\\mathcal{N}\\left(0, \\mathbf{A}\\cdot\\sigma^{2}_{a}\\right)$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathbf{\\epsilon}$ - wektor błędów ($n\\times 1$), $\\mathbf{\\epsilon}\\sim\\mathcal{N}\\left(0, \\mathbf{I}\\cdot\\sigma^{2}_{\\epsilon}\\right)$" ] }, { "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 = \\sqrt{diag(1-\\mathbf{C}^{22}\\cdot\\alpha}$)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Liczymy macierz A" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
idsiredam
100
200
300
410
532
612
745
836
\n" ], "text/latex": [ "\\begin{tabular}{lll}\n", " id & sire & dam\\\\\n", "\\hline\n", "\t 1 & 0 & 0\\\\\n", "\t 2 & 0 & 0\\\\\n", "\t 3 & 0 & 0\\\\\n", "\t 4 & 1 & 0\\\\\n", "\t 5 & 3 & 2\\\\\n", "\t 6 & 1 & 2\\\\\n", "\t 7 & 4 & 5\\\\\n", "\t 8 & 3 & 6\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| id | sire | dam |\n", "|---|---|---|\n", "| 1 | 0 | 0 |\n", "| 2 | 0 | 0 |\n", "| 3 | 0 | 0 |\n", "| 4 | 1 | 0 |\n", "| 5 | 3 | 2 |\n", "| 6 | 1 | 2 |\n", "| 7 | 4 | 5 |\n", "| 8 | 3 | 6 |\n", "\n" ], "text/plain": [ " id sire dam\n", "[1,] 1 0 0 \n", "[2,] 2 0 0 \n", "[3,] 3 0 0 \n", "[4,] 4 1 0 \n", "[5,] 5 3 2 \n", "[6,] 6 1 2 \n", "[7,] 7 4 5 \n", "[8,] 8 3 6 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "library(AGHmatrix)\n", "\n", "id = 1:8\n", "sire = c(0, 0, 0, 1, 3, 1, 4, 3)\n", "dam = c(0, 0, 0, 0, 2, 2, 5, 6)\n", "\n", "(ped = cbind(id, sire, dam))" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Verifying conflicting data... \n", "Organizing data... \n", "Your data was chronologically organized with success. \n", "Constructing matrix A using ploidy = 2 \n", "Completed! Time = 0 minutes \n" ] } ], "source": [ "A = Amatrix(ped)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
12345678
1.00 0.00 0.00 0.5000.0000.50 0.25 0.250
0.00 1.00 0.00 0.0000.5000.50 0.25 0.250
0.00 0.00 1.00 0.0000.5000.00 0.25 0.500
0.50 0.00 0.00 1.0000.0000.25 0.50 0.125
0.00 0.50 0.50 0.0001.0000.25 0.50 0.375
0.50 0.50 0.00 0.2500.2501.00 0.25 0.500
0.25 0.25 0.25 0.5000.5000.25 1.00 0.250
0.25 0.25 0.50 0.1250.3750.50 0.25 1.000
\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": [ "(A = as.matrix(A))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. Definiujemy kolejne źródła informacji" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dla osobników z przynajmniej jednym znanym rodzicem posiadamy wartość fenotypową" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\n", "
4.52.93.93.55
\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": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
10
01
01
10
10
\n" ], "text/latex": [ "\\begin{tabular}{ll}\n", "\t 1 & 0\\\\\n", "\t 0 & 1\\\\\n", "\t 0 & 1\\\\\n", "\t 1 & 0\\\\\n", "\t 1 & 0\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 1 | 0 |\n", "| 0 | 1 |\n", "| 0 | 1 |\n", "| 1 | 0 |\n", "| 1 | 0 |\n", "\n" ], "text/plain": [ " [,1] [,2]\n", "[1,] 1 0 \n", "[2,] 0 1 \n", "[3,] 0 1 \n", "[4,] 1 0 \n", "[5,] 1 0 " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
factor(sex)0factor(sex)1
01
10
10
01
01
\n" ], "text/latex": [ "\\begin{tabular}{r|ll}\n", " factor(sex)0 & factor(sex)1\\\\\n", "\\hline\n", "\t 0 & 1\\\\\n", "\t 1 & 0\\\\\n", "\t 1 & 0\\\\\n", "\t 0 & 1\\\\\n", "\t 0 & 1\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| factor(sex)0 | factor(sex)1 |\n", "|---|---|\n", "| 0 | 1 |\n", "| 1 | 0 |\n", "| 1 | 0 |\n", "| 0 | 1 |\n", "| 0 | 1 |\n", "\n" ], "text/plain": [ " factor(sex)0 factor(sex)1\n", "1 0 1 \n", "2 1 0 \n", "3 1 0 \n", "4 0 1 \n", "5 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", "X\n", "\n", "model.matrix(~factor(sex) - 1)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 2
    2. \n", "\t
  2. 2
    3. \n", "\t
  3. 1
    4. \n", "\t
  4. 1
    5. \n", "\t
  5. 2
    6. \n", "\t
  6. 3
    7. \n", "\t
  7. 3
    8. \n", "\t
  8. 3
    9. \n", "\t
  9. 2
    10. \n", "\t
  10. 2
    11. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 2\n", "\\item 2\n", "\\item 1\n", "\\item 1\n", "\\item 2\n", "\\item 3\n", "\\item 3\n", "\\item 3\n", "\\item 2\n", "\\item 2\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 2\n", "2. 2\n", "3. 1\n", "4. 1\n", "5. 2\n", "6. 3\n", "7. 3\n", "8. 3\n", "9. 2\n", "10. 2\n", "\n", "\n" ], "text/plain": [ " [1] 2 2 1 1 2 3 3 3 2 2" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\t
1
\n", "\t\t
1
\n", "\t
2
\n", "\t\t
1
\n", "\t
3
\n", "\t\t
1
\n", "\t
4
\n", "\t\t
1
\n", "\t
5
\n", "\t\t
1
\n", "\t
6
\n", "\t\t
1
\n", "\t
7
\n", "\t\t
1
\n", "\t
8
\n", "\t\t
1
\n", "\t
9
\n", "\t\t
1
\n", "\t
10
\n", "\t\t
1
\n", "
\n" ], "text/latex": [ "\\begin{description*}\n", "\\item[1] 1\n", "\\item[2] 1\n", "\\item[3] 1\n", "\\item[4] 1\n", "\\item[5] 1\n", "\\item[6] 1\n", "\\item[7] 1\n", "\\item[8] 1\n", "\\item[9] 1\n", "\\item[10] 1\n", "\\end{description*}\n" ], "text/markdown": [ "1\n", ": 12\n", ": 13\n", ": 14\n", ": 15\n", ": 16\n", ": 17\n", ": 18\n", ": 19\n", ": 110\n", ": 1\n", "\n" ], "text/plain": [ " 1 2 3 4 5 6 7 8 9 10 \n", " 1 1 1 1 1 1 1 1 1 1 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "k = sample(1:3, 10, replace = T)\n", "k\n", "D = model.matrix(~factor(k) - 1)\n", "rowSums(D)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
00010000
00001000
00000100
00000010
00000001
\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": [ "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": 7, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\t
$C
\n", "\t\t
\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\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
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
0 0 3.666667e+00 1.000000e+00 5.551115e-15-1.333333e+00 8.881784e-16-2.000000e+00-1.065814e-14-3.712308e-15
0 0 1.000000e+00 4.000000e+00 1.000000e+00-1.554312e-15-2.000000e+00-2.000000e+00 1.776357e-15-2.289835e-16
0 0 2.220446e-15 1.000000e+00 4.000000e+00 4.440892e-16-2.000000e+00 1.000000e+00-7.771561e-16-2.000000e+00
1 0 -1.333333e+00-3.330669e-15-2.553513e-15 4.666667e+00 1.000000e+00 1.332268e-15-2.000000e+00 5.065393e-16
0 1 2.220446e-16-2.000000e+00-2.000000e+00 1.000000e+00 6.000000e+00-1.776357e-15-2.000000e+00-1.942890e-15
0 1 -2.000000e+00-2.000000e+00 1.000000e+00-2.331468e-15-1.554312e-15 6.000000e+00 2.331468e-15-2.000000e+00
1 0 1.110223e-16-8.881784e-16-2.220446e-16-2.000000e+00-2.000000e+00-1.110223e-16 5.000000e+00 2.706169e-16
1 0 -2.428613e-15-7.667478e-15-2.000000e+00-5.197232e-15 2.317591e-15-2.000000e+00 4.857226e-15 5.000000e+00
\n", "
\n", "\t
$est
\n", "\t\t
\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
4.358502330
3.404430006
0.098444576
-0.018770099
-0.041084203
-0.008663123
-0.185732099
0.176872088
-0.249458555
0.182614688
\n", "
\n", "
\n" ], "text/latex": [ "\\begin{description}\n", "\\item[\\$C] \\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.551115e-15 & -1.333333e+00 & 8.881784e-16 & -2.000000e+00 & -1.065814e-14 & -3.712308e-15\\\\\n", "\t 0 & 0 & 1.000000e+00 & 4.000000e+00 & 1.000000e+00 & -1.554312e-15 & -2.000000e+00 & -2.000000e+00 & 1.776357e-15 & -2.289835e-16\\\\\n", "\t 0 & 0 & 2.220446e-15 & 1.000000e+00 & 4.000000e+00 & 4.440892e-16 & -2.000000e+00 & 1.000000e+00 & -7.771561e-16 & -2.000000e+00\\\\\n", "\t 1 & 0 & -1.333333e+00 & -3.330669e-15 & -2.553513e-15 & 4.666667e+00 & 1.000000e+00 & 1.332268e-15 & -2.000000e+00 & 5.065393e-16\\\\\n", "\t 0 & 1 & 2.220446e-16 & -2.000000e+00 & -2.000000e+00 & 1.000000e+00 & 6.000000e+00 & -1.776357e-15 & -2.000000e+00 & -1.942890e-15\\\\\n", "\t 0 & 1 & -2.000000e+00 & -2.000000e+00 & 1.000000e+00 & -2.331468e-15 & -1.554312e-15 & 6.000000e+00 & 2.331468e-15 & -2.000000e+00\\\\\n", "\t 1 & 0 & 1.110223e-16 & -8.881784e-16 & -2.220446e-16 & -2.000000e+00 & -2.000000e+00 & -1.110223e-16 & 5.000000e+00 & 2.706169e-16\\\\\n", "\t 1 & 0 & -2.428613e-15 & -7.667478e-15 & -2.000000e+00 & -5.197232e-15 & 2.317591e-15 & -2.000000e+00 & 4.857226e-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": [ "$C\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.551115e-15 | -1.333333e+00 | 8.881784e-16 | -2.000000e+00 | -1.065814e-14 | -3.712308e-15 |\n", "| 0 | 0 | 1.000000e+00 | 4.000000e+00 | 1.000000e+00 | -1.554312e-15 | -2.000000e+00 | -2.000000e+00 | 1.776357e-15 | -2.289835e-16 |\n", "| 0 | 0 | 2.220446e-15 | 1.000000e+00 | 4.000000e+00 | 4.440892e-16 | -2.000000e+00 | 1.000000e+00 | -7.771561e-16 | -2.000000e+00 |\n", "| 1 | 0 | -1.333333e+00 | -3.330669e-15 | -2.553513e-15 | 4.666667e+00 | 1.000000e+00 | 1.332268e-15 | -2.000000e+00 | 5.065393e-16 |\n", "| 0 | 1 | 2.220446e-16 | -2.000000e+00 | -2.000000e+00 | 1.000000e+00 | 6.000000e+00 | -1.776357e-15 | -2.000000e+00 | -1.942890e-15 |\n", "| 0 | 1 | -2.000000e+00 | -2.000000e+00 | 1.000000e+00 | -2.331468e-15 | -1.554312e-15 | 6.000000e+00 | 2.331468e-15 | -2.000000e+00 |\n", "| 1 | 0 | 1.110223e-16 | -8.881784e-16 | -2.220446e-16 | -2.000000e+00 | -2.000000e+00 | -1.110223e-16 | 5.000000e+00 | 2.706169e-16 |\n", "| 1 | 0 | -2.428613e-15 | -7.667478e-15 | -2.000000e+00 | -5.197232e-15 | 2.317591e-15 | -2.000000e+00 | 4.857226e-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": [ "$C\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.551115e-15 -1.333333e+00\n", " [4,] 0 0 1.000000e+00 4.000000e+00 1.000000e+00 -1.554312e-15\n", " [5,] 0 0 2.220446e-15 1.000000e+00 4.000000e+00 4.440892e-16\n", " [6,] 1 0 -1.333333e+00 -3.330669e-15 -2.553513e-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 -2.331468e-15\n", " [9,] 1 0 1.110223e-16 -8.881784e-16 -2.220446e-16 -2.000000e+00\n", "[10,] 1 0 -2.428613e-15 -7.667478e-15 -2.000000e+00 -5.197232e-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,] 8.881784e-16 -2.000000e+00 -1.065814e-14 -3.712308e-15\n", " [4,] -2.000000e+00 -2.000000e+00 1.776357e-15 -2.289835e-16\n", " [5,] -2.000000e+00 1.000000e+00 -7.771561e-16 -2.000000e+00\n", " [6,] 1.000000e+00 1.332268e-15 -2.000000e+00 5.065393e-16\n", " [7,] 6.000000e+00 -1.776357e-15 -2.000000e+00 -1.942890e-15\n", " [8,] -1.554312e-15 6.000000e+00 2.331468e-15 -2.000000e+00\n", " [9,] -2.000000e+00 -1.110223e-16 5.000000e+00 2.706169e-16\n", "[10,] 2.317591e-15 -2.000000e+00 4.857226e-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, sigma2_a, sigma2_e) {\n", " alpha = sigma2_e / sigma2_a\n", " invA = ginv(A)\n", " C = rbind(\n", " cbind(t(X)%*%X, t(X)%*%Z),\n", " cbind(t(Z)%*%X, t(Z)%*%Z+invA*c(alpha))\n", " )\n", " rhs = rbind(t(X)%*%y, t(Z)%*%y)\n", " invC = ginv(C)\n", " estimators = invC%*%rhs\n", " list(C = C, est = estimators)\n", "}\n", "\n", "mme(y, X, Z, A, 20, 40)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "-0.00572209096278661" ], "text/latex": [ "-0.00572209096278661" ], "text/markdown": [ "-0.00572209096278661" ], "text/plain": [ "[1] -0.005722091" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "0.157163480373342" ], "text/latex": [ "0.157163480373342" ], "text/markdown": [ "0.157163480373342" ], "text/plain": [ "[1] 0.1571635" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "est = mme(y, X, Z, A, 20, 40)$est[3:10]\n", "(m = mean(est))\n", "(s = sd(est))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "100" ], "text/latex": [ "100" ], "text/markdown": [ "100" ], "text/plain": [ "[1] 100" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "10" ], "text/latex": [ "10" ], "text/markdown": [ "10" ], "text/plain": [ "[1] 10" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 106.627918039179
    2. \n", "\t
  2. 99.1697811662678
    3. \n", "\t
  3. 97.7499790740002
    4. \n", "\t
  4. 99.8128679962948
    5. \n", "\t
  5. 88.5463208053011
    6. \n", "\t
  6. 111.618104804649
    7. \n", "\t
  7. 84.491533065325
    8. \n", "\t
  8. 111.983495048983
    9. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 106.627918039179\n", "\\item 99.1697811662678\n", "\\item 97.7499790740002\n", "\\item 99.8128679962948\n", "\\item 88.5463208053011\n", "\\item 111.618104804649\n", "\\item 84.491533065325\n", "\\item 111.983495048983\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 106.627918039179\n", "2. 99.1697811662678\n", "3. 97.7499790740002\n", "4. 99.8128679962948\n", "5. 88.5463208053011\n", "6. 111.618104804649\n", "7. 84.491533065325\n", "8. 111.983495048983\n", "\n", "\n" ], "text/plain": [ "[1] 106.62792 99.16978 97.74998 99.81287 88.54632 111.61810 84.49153\n", "[8] 111.98350" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "est = 10*(est - m) / s + 100\n", "mean(est)\n", "sd(est)\n", "est" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 0.240440381554564
    2. \n", "\t
  2. 0.125732088645319
    3. \n", "\t
  3. 0.295097324496044
    4. \n", "\t
  4. 0.38031525456248
    5. \n", "\t
  5. 0.379191912532635
    6. \n", "\t
  6. 0.33975648445238
    7. \n", "\t
  7. 0.341009757998517
    8. \n", "\t
  8. 0.394045310883079
    9. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.240440381554564\n", "\\item 0.125732088645319\n", "\\item 0.295097324496044\n", "\\item 0.38031525456248\n", "\\item 0.379191912532635\n", "\\item 0.33975648445238\n", "\\item 0.341009757998517\n", "\\item 0.394045310883079\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.240440381554564\n", "2. 0.125732088645319\n", "3. 0.295097324496044\n", "4. 0.38031525456248\n", "5. 0.379191912532635\n", "6. 0.33975648445238\n", "7. 0.341009757998517\n", "8. 0.394045310883079\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)$C)\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": [ "4. Estymacja parametrów wariancji" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\n", "
0.678
\n" ], "text/latex": [ "\\begin{tabular}{l}\n", "\t 0.678\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 0.678 |\n", "\n" ], "text/plain": [ " [,1] \n", "[1,] 0.678" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "var(y)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "sigma_a = 121.01 #starting value for random effect\n", "sigma_e = 10.51 #starting value for error variance" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "EM = function(y, X, Z, A, sigma_a, sigma_e) {\n", " n = nrow(X)\n", " p = ncol(X) \n", " q = nrow(A) \n", " \n", " t = 1 #iteration number 1\n", " tmp = 0.1 #test for convergance\n", " \n", " while (tmp > 0.00001) {\n", " mme_new = mme(y, X, Z, A, sigma_a, sigma_e)\n", " C_new = ginv(mme_new$C)\n", " Ck = C_new[(p+1):(p+q), (p+1):(p+q)]\n", " mme2 = mme_new$est\n", " \n", " a = as.matrix(mme2[(p+1):(p+q)])\n", " sigma_a_new = (t(a)%*%ginv(A)%*%a + sum(diag(ginv(A)%*%Ck))*c(sigma_e))/q\n", " \n", " res = as.matrix(y-X%*%as.matrix(mme2[1:p]) - Z%*%as.matrix(mme2[(p+1):(p+q)]))\n", " X.tmp1 = cbind(X,Z) %*% C_new\n", " X.tmp2 = t(cbind(X,Z))\n", " sigma_e_new = (t(res)%*%res + sum(diag(X.tmp1%*%X.tmp2))*c(sigma_e))/n\n", " \n", " tmp = max(abs(sigma_a - sigma_a_new), abs(sigma_e - sigma_e_new))\n", " sigma_a = sigma_a_new\n", " sigma_e = sigma_e_new\n", " \n", " t = t + 1\n", " }\n", " list(t = t, sigma_a = sigma_a, sigma_e = sigma_e)\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$y=X\\beta+Za+\\epsilon$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\epsilon = y - X\\beta - Za$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $\\sigma^{2}_{\\epsilon[t+1]} = \\frac{\\widehat{\\epsilon}^{'}_{[t]}\\widehat{\\epsilon}_{[t]} + tr([X, Z]C^{22}_{[t]}[X, Z]^{'})\\cdot\\sigma^{2}_{\\epsilon[t]}}{n}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $\\sigma^{2}_{a[t+1]} = \\frac{\\widehat{a}^{'}_{[t]}A^{-1}\\widehat{a}_{[t]} + tr(A^{-1}C^{22}_{[t]})\\cdot\\sigma^{2}_{\\epsilon[t]}}{q}$" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
\n", "\t
$t
\n", "\t\t
658
\n", "\t
$sigma_a
\n", "\t\t
\n", "\n", "\t\n", "\n", "
0.6694982
\n", "
\n", "\t
$sigma_e
\n", "\t\t
\n", "\n", "\t\n", "\n", "
0.005737211
\n", "
\n", "
\n" ], "text/latex": [ "\\begin{description}\n", "\\item[\\$t] 658\n", "\\item[\\$sigma\\_a] \\begin{tabular}{l}\n", "\t 0.6694982\\\\\n", "\\end{tabular}\n", "\n", "\\item[\\$sigma\\_e] \\begin{tabular}{l}\n", "\t 0.005737211\\\\\n", "\\end{tabular}\n", "\n", "\\end{description}\n" ], "text/markdown": [ "$t\n", ": 658\n", "$sigma_a\n", ": \n", "| 0.6694982 |\n", "\n", "\n", "$sigma_e\n", ": \n", "| 0.005737211 |\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "$t\n", "[1] 658\n", "\n", "$sigma_a\n", " [,1]\n", "[1,] 0.6694982\n", "\n", "$sigma_e\n", " [,1]\n", "[1,] 0.005737211\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "\t\n", "\n", "
0.6752354
\n" ], "text/latex": [ "\\begin{tabular}{l}\n", "\t 0.6752354\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 0.6752354 |\n", "\n" ], "text/plain": [ " [,1] \n", "[1,] 0.6752354" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "\t\n", "\n", "
0.678
\n" ], "text/latex": [ "\\begin{tabular}{l}\n", "\t 0.678\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 0.678 |\n", "\n" ], "text/plain": [ " [,1] \n", "[1,] 0.678" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "(wyniki = EM(y, X, Z, A, sigma_a, sigma_e))\n", "wyniki$sigma_a+wyniki$sigma_e\n", "var(y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "5. Odziedziczalność" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "99.1503391400188" ], "text/latex": [ "99.1503391400188" ], "text/markdown": [ "99.1503391400188" ], "text/plain": [ "[1] 99.15034" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "(h2 = 0.6694982 / (0.6694982 + 0.005737211)) * 100" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "-0.0555680479925584" ], "text/latex": [ "-0.0555680479925584" ], "text/markdown": [ "-0.0555680479925584" ], "text/plain": [ "[1] -0.05556805" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "0.501887852905857" ], "text/latex": [ "0.501887852905857" ], "text/markdown": [ "0.501887852905857" ], "text/plain": [ "[1] 0.5018879" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "100" ], "text/latex": [ "100" ], "text/markdown": [ "100" ], "text/plain": [ "[1] 100" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "10" ], "text/latex": [ "10" ], "text/markdown": [ "10" ], "text/plain": [ "[1] 10" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
estest2
106.62792106.68655
99.16978 98.37324
97.74998 99.40987
99.81287102.17470
88.54632 89.43013
111.61810109.25737
84.49153 82.57042
111.98350112.09773
\n" ], "text/latex": [ "\\begin{tabular}{ll}\n", " est & est2\\\\\n", "\\hline\n", "\t 106.62792 & 106.68655\\\\\n", "\t 99.16978 & 98.37324\\\\\n", "\t 97.74998 & 99.40987\\\\\n", "\t 99.81287 & 102.17470\\\\\n", "\t 88.54632 & 89.43013\\\\\n", "\t 111.61810 & 109.25737\\\\\n", "\t 84.49153 & 82.57042\\\\\n", "\t 111.98350 & 112.09773\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| est | est2 |\n", "|---|---|\n", "| 106.62792 | 106.68655 |\n", "| 99.16978 | 98.37324 |\n", "| 97.74998 | 99.40987 |\n", "| 99.81287 | 102.17470 |\n", "| 88.54632 | 89.43013 |\n", "| 111.61810 | 109.25737 |\n", "| 84.49153 | 82.57042 |\n", "| 111.98350 | 112.09773 |\n", "\n" ], "text/plain": [ " est est2 \n", "[1,] 106.62792 106.68655\n", "[2,] 99.16978 98.37324\n", "[3,] 97.74998 99.40987\n", "[4,] 99.81287 102.17470\n", "[5,] 88.54632 89.43013\n", "[6,] 111.61810 109.25737\n", "[7,] 84.49153 82.57042\n", "[8,] 111.98350 112.09773" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "est2 = mme(y, X, Z, A, 0.6694982, 0.005737211)$est[3:10]\n", "(m = mean(est2))\n", "(s = sd(est2))\n", "\n", "est2 = 10*(est2 - m) / s + 100\n", "mean(est2)\n", "sd(est2)\n", "cbind(est, est2)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 0.189936800414011
    2. \n", "\t
  2. 0.0776489403630648
    3. \n", "\t
  3. 0.291999858362966
    4. \n", "\t
  4. 0.482421097371713
    5. \n", "\t
  5. 0.459810563736023
    6. \n", "\t
  6. 0.457015682570914
    7. \n", "\t
  7. 0.479626216206605
    8. \n", "\t
  8. 0.483239239421287
    9. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.189936800414011\n", "\\item 0.0776489403630648\n", "\\item 0.291999858362966\n", "\\item 0.482421097371713\n", "\\item 0.459810563736023\n", "\\item 0.457015682570914\n", "\\item 0.479626216206605\n", "\\item 0.483239239421287\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.189936800414011\n", "2. 0.0776489403630648\n", "3. 0.291999858362966\n", "4. 0.482421097371713\n", "5. 0.459810563736023\n", "6. 0.457015682570914\n", "7. 0.479626216206605\n", "8. 0.483239239421287\n", "\n", "\n" ], "text/plain": [ "[1] 0.18993680 0.07764894 0.29199986 0.48242110 0.45981056 0.45701568 0.47962622\n", "[8] 0.48323924" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 0.435817393427581
    2. \n", "\t
  2. 0.278655594530354
    3. \n", "\t
  3. 0.54037011238869
    4. \n", "\t
  4. 0.694565401795765
    5. \n", "\t
  5. 0.678093329664894
    6. \n", "\t
  6. 0.676029350376826
    7. \n", "\t
  7. 0.692550515274232
    8. \n", "\t
  8. 0.695154111993368
    9. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.435817393427581\n", "\\item 0.278655594530354\n", "\\item 0.54037011238869\n", "\\item 0.694565401795765\n", "\\item 0.678093329664894\n", "\\item 0.676029350376826\n", "\\item 0.692550515274232\n", "\\item 0.695154111993368\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.435817393427581\n", "2. 0.278655594530354\n", "3. 0.54037011238869\n", "4. 0.694565401795765\n", "5. 0.678093329664894\n", "6. 0.676029350376826\n", "7. 0.692550515274232\n", "8. 0.695154111993368\n", "\n", "\n" ], "text/plain": [ "[1] 0.4358174 0.2786556 0.5403701 0.6945654 0.6780933 0.6760294 0.6925505\n", "[8] 0.6951541" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
rr_2
0.24044040.4358174
0.12573210.2786556
0.29509730.5403701
0.38031530.6945654
0.37919190.6780933
0.33975650.6760294
0.34100980.6925505
0.39404530.6951541
\n" ], "text/latex": [ "\\begin{tabular}{ll}\n", " r & r\\_2\\\\\n", "\\hline\n", "\t 0.2404404 & 0.4358174\\\\\n", "\t 0.1257321 & 0.2786556\\\\\n", "\t 0.2950973 & 0.5403701\\\\\n", "\t 0.3803153 & 0.6945654\\\\\n", "\t 0.3791919 & 0.6780933\\\\\n", "\t 0.3397565 & 0.6760294\\\\\n", "\t 0.3410098 & 0.6925505\\\\\n", "\t 0.3940453 & 0.6951541\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| r | r_2 |\n", "|---|---|\n", "| 0.2404404 | 0.4358174 |\n", "| 0.1257321 | 0.2786556 |\n", "| 0.2950973 | 0.5403701 |\n", "| 0.3803153 | 0.6945654 |\n", "| 0.3791919 | 0.6780933 |\n", "| 0.3397565 | 0.6760294 |\n", "| 0.3410098 | 0.6925505 |\n", "| 0.3940453 | 0.6951541 |\n", "\n" ], "text/plain": [ " r r_2 \n", "[1,] 0.2404404 0.4358174\n", "[2,] 0.1257321 0.2786556\n", "[3,] 0.2950973 0.5403701\n", "[4,] 0.3803153 0.6945654\n", "[5,] 0.3791919 0.6780933\n", "[6,] 0.3397565 0.6760294\n", "[7,] 0.3410098 0.6925505\n", "[8,] 0.3940453 0.6951541" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "C = as.matrix(mme(y, X, Z, A, 0.6694982, 0.005737211)$C)\n", "invC = ginv(C)\n", "\n", "invC22 = invC[3:10, 3:10]\n", "(r2 = diag(1 - invC22*(0.005737211/0.6694982)))\n", "(r_2 = sqrt(r2))\n", "\n", "cbind(r, r_2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "6. Test Walda z nowymi $\\sigma^{2}_{a}$ i $\\sigma^{2}_{\\epsilon}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$W=\\frac{\\hat{\\beta}}{se(\\hat{\\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(\\beta) = \\left(X^{T}V^{-1}X \\right)^{-1}$, \n", "$se(\\beta) = \\sqrt{diag(var(\\beta))}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "6.1 stare parametry wariancji" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
60.0 0.0 5 10 2.5
0.060.0 5 10 7.5
5.0 5.060 5 10.0
10.010.0 5 60 5.0
2.5 7.510 5 60.0
\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": [ "\n", "\n", "\t\n", "\t\n", "\n", "
23.822439 6.292085
6.29208532.098346
\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": [ "
    \n", "\t
  1. 4.88082357807458
    2. \n", "\t
  2. 5.66554023294586
    3. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 4.88082357807458\n", "\\item 5.66554023294586\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 4.88082357807458\n", "2. 5.66554023294586\n", "\n", "\n" ], "text/plain": [ "[1] 4.880824 5.665540" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 0.892985017822407
    2. \n", "\t
  2. 0.600901214346598
    3. \n", "
\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": [ "
    \n", "\t
  1. 0.371865196
    2. \n", "\t
  2. 0.5479057844
    3. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.371865196\n", "\\item 0.5479057844\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.371865196\n", "2. 0.5479057844\n", "\n", "\n" ], "text/plain": [ "[1] 0.3718652 0.5479058" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "G = A*c(20)\n", "R = diag(5)*c(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 = round(2*pnorm(abs(testWalda), lower.tail = FALSE), digits = 10))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "6.2 nowe parametry wariancji" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
0.67523530.00000000.16737420.33474840.0836871
0.00000000.67523530.16737420.33474840.2510613
0.16737420.16737420.67523530.16737420.3347484
0.33474840.33474840.16737420.67523530.1673742
0.08368710.25106130.33474840.16737420.6752353
\n" ], "text/latex": [ "\\begin{tabular}{lllll}\n", "\t 0.6752353 & 0.0000000 & 0.1673742 & 0.3347484 & 0.0836871\\\\\n", "\t 0.0000000 & 0.6752353 & 0.1673742 & 0.3347484 & 0.2510613\\\\\n", "\t 0.1673742 & 0.1673742 & 0.6752353 & 0.1673742 & 0.3347484\\\\\n", "\t 0.3347484 & 0.3347484 & 0.1673742 & 0.6752353 & 0.1673742\\\\\n", "\t 0.0836871 & 0.2510613 & 0.3347484 & 0.1673742 & 0.6752353\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 0.6752353 | 0.0000000 | 0.1673742 | 0.3347484 | 0.0836871 |\n", "| 0.0000000 | 0.6752353 | 0.1673742 | 0.3347484 | 0.2510613 |\n", "| 0.1673742 | 0.1673742 | 0.6752353 | 0.1673742 | 0.3347484 |\n", "| 0.3347484 | 0.3347484 | 0.1673742 | 0.6752353 | 0.1673742 |\n", "| 0.0836871 | 0.2510613 | 0.3347484 | 0.1673742 | 0.6752353 |\n", "\n" ], "text/plain": [ " [,1] [,2] [,3] [,4] [,5] \n", "[1,] 0.6752353 0.0000000 0.1673742 0.3347484 0.0836871\n", "[2,] 0.0000000 0.6752353 0.1673742 0.3347484 0.2510613\n", "[3,] 0.1673742 0.1673742 0.6752353 0.1673742 0.3347484\n", "[4,] 0.3347484 0.3347484 0.1673742 0.6752353 0.1673742\n", "[5,] 0.0836871 0.2510613 0.3347484 0.1673742 0.6752353" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\n", "
0.34510430.2065161
0.20651610.3626310
\n" ], "text/latex": [ "\\begin{tabular}{ll}\n", "\t 0.3451043 & 0.2065161\\\\\n", "\t 0.2065161 & 0.3626310\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 0.3451043 | 0.2065161 |\n", "| 0.2065161 | 0.3626310 |\n", "\n" ], "text/plain": [ " [,1] [,2] \n", "[1,] 0.3451043 0.2065161\n", "[2,] 0.2065161 0.3626310" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 0.58745576587621
    2. \n", "\t
  2. 0.602188504574106
    3. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.58745576587621\n", "\\item 0.602188504574106\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.58745576587621\n", "2. 0.602188504574106\n", "\n", "\n" ], "text/plain": [ "[1] 0.5874558 0.6021885" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 7.56094185962023
    2. \n", "\t
  2. 5.79304271276703
    3. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 7.56094185962023\n", "\\item 5.79304271276703\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 7.56094185962023\n", "2. 5.79304271276703\n", "\n", "\n" ], "text/plain": [ "[1] 7.560942 5.793043" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 0
    2. \n", "\t
  2. 6.9e-09
    3. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0\n", "\\item 6.9e-09\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0\n", "2. 6.9e-09\n", "\n", "\n" ], "text/plain": [ "[1] 0.0e+00 6.9e-09" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "G = A*0.6694968\n", "R = diag(5)*0.005738462\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, 0.6694968, 0.005738462)$est[1:2] / seB)\n", "(p_value = round(2*pnorm(abs(testWalda), lower.tail = FALSE), digits = 10))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "7. Korelacja oceny " ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "0.98640763758191" ], "text/latex": [ "0.98640763758191" ], "text/markdown": [ "0.98640763758191" ], "text/plain": [ "[1] 0.9864076" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "0.977061470732081" ], "text/latex": [ "0.977061470732081" ], "text/markdown": [ "0.977061470732081" ], "text/plain": [ "[1] 0.9770615" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cor(est, est2)\n", "cor(r, r_2)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAMFBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////QFLu4AAAACXBIWXMAABJ0AAAS\ndAHeZh94AAAT2ElEQVR4nO3dbUPayhqG0QlgQAzw///thuALWEvdcGdCYK0Plt2e04cIVzGT\nEcsOuFkZ+w7AIxASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAg\nJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkC\nhAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFB\ngJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAk\nCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKE\nBAFCggAhQYCQIEBIECAkCKgQUoGJueJZng9nhBGQJCQIEBIECAkChAQBQoIAIUGAkCBASBAg\nJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBAS/NaF7ygXEvxOX9HfUhIS/E45+fiXP7zi7xuU\nkLg75duvP//pFX/hgITE3RESBAgJEpwjQYBVO4hwHQmGJSQIEBIECAkChAQBQoIAIUGAkCBA\nSBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAVVDelsuysGifRtqBIyiYkjb\nWfkyH2QEjKRiSG1pXrv+1mbdlHaIETCSiiE1pfu83ZVmiBEwkoohnb1L5d/fsvKWETASr0gQ\nUPccab3pbzlH4tHUXP6en6zazbaDjIBx1L2O1PbXkZrF0nUkHoudDRAgJAgQEgQICQKEBAFV\ndzacGWIEjKRiSKvLIf26Mrg/Nb+065rL3zwRGAHjqHqO1F3eGJQYAaOou9iwOtm3OtAIGINV\nOwgQEgQICQKEBAFCggAhQYAtQvAvv9hrcz9bhBIjIK9/rv4rJVuE4LJy8vEf/6Mr/t5r2CLE\n9JRvv17+X13xF1/BFiEm5x5DuqMR8DtCgoT7O0e6pxHwS3e3andXI+DX7us60n2NgCQhQYCQ\nIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQ\nEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQB\nQoIAITEtpdzls0NITElf0T2mJCSmpJx8vCtCYkLKt1/vh5CYECHd3QimSEh3N4JJco50byOY\nJKt29zaCiXId6b5GQJKQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIE1A9pNStl\nsR50BNRWMaTj9vd56bWDjICR1A6pLe12t9u0ZTXECBhJ7ZCasj3c3pbZECNgJLVD+vj2xsvf\n5igkJqZ2SC8fITVDjICRVA1psVyty+v+5ra9vNogJCamakhH/c1mO8QIGEnN60hdt1otFv2S\nQ3uxIyExNXY2QICQIKBqSG/LRX+WtGjfhhoBo6gY0nZWvswHGQEjqRhSW5rXrr+1WTeWv3ko\nFUNqSvd5u3NBlodSfff3T//x/jsnrhwBI/GKBAF1z5HWm/6WcyQeTc3l7/nJ124zW4R4JHWv\nI7X9daRmsXQdicdiZwMECAkChAQBQoIAIUHACN8h+4vNC0JiYiqGtBISD6vqt5o3l795IjAC\nxlH1HKn7xzsVB0bAKOouNqxO9q0ONALGYNUOAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQB\nQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAg\nQEgQICQIEBIECAkChAQBQoIAIUGAkPimFJ/+/09InOkrktL/JiTOlJOP/J6QOFW+/covCYlT\nQrqSkDglpCsJiTPOka4jJM5YtbuOkPjGdaRrCAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKE\nBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQUDWkt+WiHCzat6FGwCgqhrSd\nlS/zQUbASCqG1JbmtetvbdZNaYcYASOpGFJTus/bXWmGGAEjqRjS2TvhXn5bXCExMV6RIKDu\nOdJ6099yjsSjqbn8PT9ZtZttBxkB46h7HantryM1i6XrSDwWOxsgQEgQYIsQBNgiBAG2CEGA\nC7IQcD9bhMqpK0fASLwiQYAtQhBgixAE2CIEAXY2QICQIEBID8/VhBqE9OD6iqQ0OCE9uHLy\nkeFU3dnw680LHveU8u1XBlIxpJWQ6hNSJTW/tOuay988ERjBN0KqpOo5Und5Y1BiBN84R6qj\n7mLD6mTf6kAjOGfVrg6rdg/PdaQahAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQ\nsBkvQEjPzvbwCCE9O9+wFCGkJ+dbaDOE9OSElCGkJyekjFtC2r6UMl+//2b0gfCo1uMcKeKG\nkLbN8UeUH39TSBNl1S7ihpDastrXtDq+x5aQpst1pIAbQmqONzbNbCMkntwNIX20s53PhcST\nuyGkWfn48ZWzuZB4bjeEtCov77c2ZS4kntoty9/tZz3r8PmqkJiYmy7IdouPW5sXIfHM7GyA\nACFBwI0hfZ4aNU3i3vw0AiYgFNLGYgNP7YaQ1mc/yXI28r2CMd3yijQ77eht5HsFY0qdI2UJ\niYmxagcBt4a02p8bbWbhr+yExNTcGNL68LVd/w1+zpF4ZjeGNC+vu67Mdq9lHrtLOyExOYHF\nhq60vh+JJxcIaVHWQuLJ3fylXbcuzc6Xdjy52xcbSlkeXpDWsbu0ExKTc/Pyd3M4Q9rNXkP3\n54cRcP9ckIUAIUHAzSGtF/3K3SZ0f34aAXfv1pDmpX/fk9JESxISE3NjSKsy3x5C+nprrggh\nMTE3htSU7RDvwi4kJiaws0FIcGNIs/dXpM63mvPUMudI6+bwI15yhMTE3Lpqt3h/z4boVjsh\nMTWR60hlkd0hJCSmxs4GCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGA\nkCBASBAgJAgQEgTUD2k1K2Xxjx8CIyQmpmJIx7e+mx/fLKUdZASMpHZIbWm3u92mvfz2XUJi\nYmqHdHiP473t5TeUFBITUzukj7c2vvwWx0JiYmqH9PIRUjPECBhJ1ZAWy9W6HN5LctteXm0Q\nEhNTNaSj/mazHWIEjKTmdaSuW60Wi37Job3YkZCYGjsbIOB+QiqnhhkBQ6kaUtce9zXM/vXT\nK4TExNQMaXnykrMYZgSMo2JI6/Ky2e3e5otdt5qVi9tWhcTEVAxpftwd1JXlPqfLL0lCYmKq\n7/5+39RgixAPpWJIn/tVT/fcZUfASCqG1Jb52263WZSX3fZl/2GAETCSmqt279/T12wPW4Q2\ng4yAcVS9jrTapzRb7mwR4uHcz86GyiMgSUgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQ\nEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQB\nQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAg\nQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBAS\nBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkCqob0tlyUg0X7NtQIGEXFkLaz8mU+yAgYScWQ\n2tK8dv2tzbop7RAjYCQVQ2pK93m7K80QI2AkFUMq5W//ERsBI/GKBAF1z5HWm/6WcyQeTc3l\n7/nJqt1sO8gIGEfd60htfx2pWSxdR+Kx2NkAAfcTUjk1zAgYStWQuvZ4mjRbvA41AkZRM6Tl\nyUvOYpgRMI6KIa3Ly2a3e5svdt1qVtZDjICRVAxpXvol764s9zldfkkSEhMzwhahflODLUI8\nlKpbhPpXpG3fkJB4KFW3CM3fdrvNorzsti/7DwOMuJF1d641whahZrt/xjabQUbc4hevlPAX\nVa8jrfYpzZb7G017cavdSCGNN5rJu5+dDZVH/HWmkriCkL7PFBJXENL3mULiCkL6NlRHXENI\nX0Ot2nE1IZ2OlRFXEhIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBI\nECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQI\nCQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIA\nIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQ\nUD+k1ayUxXrQEVBbxZBK/3+cl147yAgYSe2Q2tJud7tNW1ZDjICR1A6pKdvD7W2ZDTECRlI7\npFJO/iM+AkZSO6SXj5CaIUbASKqGtFiu1uV1f3PbXl5tEBITUzWko/5msx1iBIyk5nWkrlut\nFot+yaG92JGQmBo7GyDgfkIqp4YZAUOpGdLmpTTLfo9Qc3ljg1ckpqZiSNvm8FqzWvYvOfNB\nRsBIKobUL3m3TXnZWv7m0VQMqen/j+W4R8gFWR5K9d3f7wsJtgjxUEZ4RTp83HpF4qGMcI50\nuBjrHInHYtUOAlxHgoD72dlQeQQkCQkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAh\nQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAg\nJAgQEgRMK6Ry+Yehw1imFFJfkZS4R5MKqdZ4+L8mFFK59IcwKiFBgJAgYEIhOUfifk0qJKt2\n3KspheQ6EndrWiHBnRISBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBA\nSBAgJAgQEgTcaUgwMVc8y/PhXGWU+zHGUAf6YDPHH33G8+vRhj7LzPFHn/H8erShzzJz/NFn\nPL8ebeizzBx/9BnPr0cb+iwzxx99xvPr0YY+y8zxR5/x/Hq0oc8yc/zRZzy/Hm3os8wcf/QZ\nz69HG/osM8cffcbz69GGPsvM8Uef8fx6tKHPMnP80Wc8vx5t6LPMHH80PA4hQYCQIEBIECAk\nCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBI4e0bZvStNv+nlz79uX/337ofP15\n831+tZmVDnT1MeHkEIc+2h9mDn60q/LnzWqP6pdxQ9o0/We52ex2Xb2Q5v2g5dfNWc2ZlQ60\n+5hwcohDH+0PMwc/2u7r7/7pkKsZN6SX0u4/tuXl8FlYVBq6KvPtbvtSut3urTTdrmvKW8WZ\ndQ50f1DHh/bkEIc+2p9mDn20nzN/Hl/PuCG9H/nhl9XxJaKCef8Z3hwabsvhq63X4UefzKxy\noPtw3z+3J4c48NH+OHPgo/2a+fP4esYNqXkPqTl8HlaVhn7UO9/tFmX/RWWN14iTmVUOdF/s\n+8iTQxz4aH+cOfDRfs38eXw944a0fP/Sbnk4+vXL/gyxwtCTl8GTm9VmVjnQ7vuxVTjaH2cO\nfLRfM38eX8/Iq3arw2pDc/hHa3E8KZ0PP3PW/4P1VjWkk5m1DrR6SD/OHP5oTw7niUNafi5m\nlfJ6WA2v8HXPsiy2u25eNaSzmXUO9D5CGv5ohbQ7fAW9f9Hfvnx9mrc1Fi37NfdF1ZBOZh4N\nf6D3EdLRkEcrpN3hK57DVbPTT3ONo9+X2yz7SU21T/nXzHeDz3wfcHKIwx/tnzPP/2DAmRfH\nV3A3y9/nvzO87lDvcX1nU2t9p6v5L8bZEtbma9VuyKP9c+b5Hww48+L4Cu5h+Xt7WP5u+hen\nGkd/nLQ6TFr2VxzWZfDFwpOZtQ70/Vl1cojDH+2fM4c/2j9Dqvaont2NmsP+0JbDjqj2eG20\n7c9K1xWGvux2b7PDWXC1a+AnM2sdaP2dDT/NHP5o/wzpCXc2vO+KOiyObo/b7ir8K/I+qf83\nclZpzf1kZq0D/XiCnRzi4Ef758zhj/aH84Jaj+rZ3ag67U/9Pt3+1mEj+KzK7obNy/4pvf4c\nWuUq8LeZFQ7041m1Pf8UD3q0f5k56NH+EFK1R/X0btQdB49JSBAgJAgQEgQICQKEBAFCggAh\nQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAg\nJAgQEgQICQKEBAFCggAhQYCQIEBID2D4n7vLvwhp+mYexPF5DKaveBDH5zGYPiHdAY/BxKxm\npTn+kPD1vJT5+tBRkdLoPALTsuizme9vrfpbZSWku+ARmJR1mW9323nZvxA1pdvtXsvMl3Z3\nwWMwKYuy3X/clsUhn49VbyHdAY/BpJQPu11byqLrjr859t1CSNNyEtJu2ex/bTZCugseg0k5\nb2bdzpwj3QmPwaQsyrftQIeIhHQHPAaT8lqa7rDyvThsDHr9XLXbjH2/ENK0zPszpMOZ0evx\nZOntkFRpxr5fT09IE7PaZ/PSvwL1Oxv2He3eZkIanZAgQEgQICQIEBIECAkChAQBQoIAIUGA\nkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQI\nEBIECAkChAQBQoKA/wB3/WKFfAkbVwAAAABJRU5ErkJggg==", "text/plain": [ "plot without title" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot(est, est2)" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\n", "Call:\n", "lm(formula = est ~ est2)\n", "\n", "Residuals:\n", " Min 1Q Median 3Q Max \n", "-2.33226 -1.18758 0.04123 1.00187 2.48655 \n", "\n", "Coefficients:\n", " Estimate Std. Error t value Pr(>|t|) \n", "(Intercept) 1.35922 6.73744 0.202 0.847 \n", "est2 0.98641 0.06708 14.705 6.21e-06 ***\n", "---\n", "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n", "\n", "Residual standard error: 1.775 on 6 degrees of freedom\n", "Multiple R-squared: 0.973,\tAdjusted R-squared: 0.9685 \n", "F-statistic: 216.2 on 1 and 6 DF, p-value: 6.214e-06\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAM1BMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD/AAD///89ODILAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAaHElEQVR4nO3d60LiyBaA0XARaUSG93/akeAFFBHIrkpVZa0fPUz3ObON\n8jWppMRuDwzWjf0BQAuEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGE\nBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGE\nBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGE\nBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGE\nBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGE\nBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAEyhNRBZR54lseHM8IIiCQkCCAkCCAk\nCCAkCCAkCCAkGO6//+7//wgJzv0nJBjqLaP/nNrBMIeOrJFgmONpnZBgiPflkZDgcf99XGYQ\nEjzssyMhwcNOrnoLCR50evdISPCYs7uwQoJH/He+m0FIcKuT7yj/1pGQ4EZ9Re8p/dhcJyS4\nTff1689NqkKCm3Rf/7yw2VtIcJPPkL4vj87+9IH/YEJCojgfT8qLHQkJbvTLZYbTP3zgv5eU\nkChPf8Hut++FFRLc6vLy6PhHD/zXhn0whYyAu/3ekZDgVtfe4kRIcJurbxUkJLjJ9bfcEhLc\n4MryqCck+NtfHQkJ/vb3O6kKCf5ywzsSCwn+cMs7ewsJrvpzedQTElxzW0dCgmtu/YEtQoLf\n3fyDj4QEv7nxtO5ASPCLOzoSEvzirp9nKSS46L6fCyskuOTOn68sJPjpnuVRT0jww90dCQl+\nuDujzCG9PC+7g+XqJdUIGOyBjnKGtJt3XxZJRsBwj3SUM6RVN/u37R+9bmbdKsUIGOr+5VEv\nY0izbvv5eNvNUoyAgR7sKGdIXffbv4SNgGEezMgrEpx4uKPMa6TNa//IGokSPXpad5Dz8vfi\n5KrdfJdkBDxsSEeZ7yOt+vtIs+Wz+0iUZkhGdjbA0bCOhAQHAzsSEgxcHvWEBMM7EhIMzyjz\nzoYzKUbA/SI6yhnS+npIN1cGkUI6ynpqt51d/+aJgBFwn4DlUS/rGml7fWNQxAi4R1RHmS82\nrE/2rSYaAbeLyshVO6YsriMhMVlhp3UHQmKiQjsSEhMVmpGQmKjgjoTEJEV3ZIsQE3Tv8uiG\nvTblbBGKGAE3uLOj/rn6V0q2CDE1j/2giYJCskWIAjz6A1uuPyltEWJaHv6BLSWFVNAIJumB\nu7BCgm8e2s1Q3hqppBFM0IM/aKK0q3ZFjWB6Hn+H/KLuI5U1gqmJ3aT6nZCYhrQdCYlpSJuR\nkJiG1B0JiSlI3pGQaF/i5VFPSLQuR0dConU5MhISrcvTkZBoW6aOhETLsiyPekKiXfk6EhLt\nypeRkGhXzo6ERKMyntYdCIkmZe5ISDQpc0ZCoknZOxISDcrfkZBoTu7lUU9INGaUjoREY0bJ\nSEg0ZqSOhERTxupISDRknOVRT0g0Y8SOhEQzRsxISDRj1I6ERBvGPK07EBItGLsjIdGCsTMS\nEi0YvyMhUb8COhIStRt9edQTEnUroyMhUbcyMhISdSulIyFRs2I6EhL1KmR51BMStSqpIyFR\nq5IyEhK16brjs6OsjoREVfqK3n4p6rTuQEjU5PjM6IrrSEjU5P2JUVxGQqIqRS6PekKiIv0T\n462j8p4gQqImx8sMBT4/hERN+ssMXYHPDyFRk0IzEhJVKfEyw5GQqEe5HQmJapR3F/aEkKhE\n0R0JiUoUnZGQqEThHQmJGpR9WncgJMpXfkdConzlZyQkyldDR0KidFV0JCTKVsHyqCckSlZL\nRyOEtJ533XKTdAStqCWjrCEdt78vut4qyQjaUk9H2UNadavdfv+66tYpRtCUijrKHtKs2x0e\n77p5ihE0pJrlUS93SB/f3nj92xyFRF0dZQ/p6SOkWYoRNKOujDKHtHxeb7p/bw93q+tXG4Q0\ndbV1lDeko/7hbJdiBG2o7LTuIOd9pO12vV4u+0sOq6sdCWnaKuzIzgaKU2FGQqI4VXaUN6SX\n52W/SlquXlKNoHZ1dpQzpN28+7JIMoLa1bg86mUMadXN/m37R6+bmcvfXFBtRzlDmnXbz8db\nN2T5qdqMRtj9felf3n/nxIMjqFrFHXlFohg1d5R5jbR57R9ZI/FDvcujXs7L34uTc7e5LUKc\nqryjzPeRVv19pNny2X0kzlSekZ0NFKH6joTE+Go/rTsQEmNroSMhMbYWMhISY2ujozG+Q/aG\nzQtCmoxGOsoZ0lpIfNPE8qiX9VvNZ9e/eSJgBFVpp6O8a6TtH+9UHDCCirSTUe6LDeuTfauJ\nRlCNljpy1Y6xNNWRkBhHQ8ujnpAYQ2sdCYkxtJaRkBhDex0JieyaO607EBKZNdmRkMisyYyE\nRGaNdiQksmq1IyGRUZvLo56QyKbhjoRENg1nJCSyabojIZFJ2x0JiSxaXh71hEQGzXckJDJo\nPiMhkcEEOhISqbV/WncgJNKaRkdCIq1pZCQk0ppKR0Iipcl0JCTSmcjyqCckUplSR0IilSll\nJCRSmVZHQiKNiXUkJFKY1PKoJyTiTa8jIRFvehkJiXhT7EhIBJvgad2BkAg10Y6ERKiJZiQk\nfui6xz/9k+1ISJzrK3o0pel2JCTOdSe/3mmqy6OekDjVffvnHSbdkZA483hIk85ISJx7OKSJ\ndyQkzj24Rpp6R0Li3ENX7aa9POoJiW/uv4+kIyExnIz2QmIwHR0IiUGc1h0JiSF09E5IDCCj\nD0LicTr6JCQepqMvQuJBlkenhMRjdHRGSDxERueExCN09I2QeICOvhMSd7M8+klI3EtHFwiJ\nO8noEiFxHx1dJCTu4bTuF0LiDjr6jZC4nYx+JSRupqPfCYlb6egKIXEby6Orsob08rzsDpar\nl1QjSERH12UMaTfvviySjCAVGf0hY0irbvZv2z963cy6VYoRJKKjv2QMadZtPx9vu1mKEaSh\noz9lDOnsnXCvvy2ukEpieXQDr0j8QUe3yLtG2rz2j6yRKiKjm+S8/L04uWo33yUZQTQd3Sbv\nfaRVfx9ptnx2H6kOTutuZWcDv9PRzYTEr2R0O1uE+I2O7mCLEL/Q0T1sEeIiy6P7uCHLJTq6\nUzlbhLpTD44giIzu5RWJn3R0N1uE+EFH97NFiG8sjx5hixDndPQQOxs4I6PHCIlTOnqQkJp3\nx90Ep3UPE1Lj+opuTElHjxNS47qTX/8gowGy7my4efOCkKJ03/55hY6GyBjSWkj53R6SjgbJ\neWq3nV3/5omAEXxza0iWRwNlXSNtr28MihjBN7etkXQ0VN6LDeuTfauJRnDupqt2MhrMVbvm\n/X0fSUfDCQkdBRDS1FkehRDSxOkohpCmTUZBhDRpOooipAlzWhdHSNOlo0BCmiwZRRLSVOko\nlJAmSkexhDRJlkfRhDRF5x15i+gAQpqgbxl9/MIAQpqe89O6O97Ugd8JaXIudeRTPpSQJub7\nZQYhxRDStPy4XCekGENC2j113WLz/puhXwhf1UQuXPW2RgoxIKTd7Pgjyo+/KaQKXLp75Kpd\niAEhrbr1W03r43tsCal8v92FdR8pwICQZscHr7P5q5AqYDdDSgNC+mhnt1gIqXwySmpASPPu\n48dXzhdCKp2O0hoQ0rp7en/02i2EVDYdJTbk8vfqs55N8HpVSLEsj5IbdEN2u/x49PokpHLp\nKD07G9onowyE1Dwd5TAwpM+l0WwW8dFcGsEwOsoiKKRXFxvKZHmUyYCQNmc/yXI+8kfFJTrK\nZcgr0vy0o5eRPyoukFE2UWukWEIKoaN8XLVrltO6nIaGtH5bG73Og8/shBRAR1kNDGlzOLfr\nv8HPGqksMsprYEiL7t9+2833/7pF2Ie0F9JwOsos4GLDtlv5fqTC6Ci3gJCW3UZIRbE8ym/w\nqd120832Tu1KoqMRDL/Y0HXPhxekTdiHtBfSIDIaw+DL37PDCmk//xf08VwYwV10NAo3ZBuj\no3EIqSmWR2MZHNJm2V+5ew36eC6N4FY6Gs3QkBZd/74n3Sy0JCE9REbjGRjSulvsDiF9vTVX\nCCE9QkcjGhjSrNuleBd2Id3Pad2oAnY2CKkEOhrXwJDm769IW99qPi4ZjSxmjbSZHX7ESxwh\n3UlHYxt61W75/p4NoVvthHQnHY0u5D5St4zdISSku1geFcDOhurpqARCqp2MiiCkyumoDEKq\nm44KIaSaWR4VQ0gV01E5hFQvGRVESNXSUUmEVCmndWURUp10VBghVUlGpRFSjXRUHCFVSEfl\nEVJ1LI9KJKTa6KhIQqqMjMokpLroqFBCqoqOSpU/pPW865Z//BAYIV1keVSujCEd3/pucXyz\nlFWSEW3TUcFyh7TqVrv9/nV1/e27hHSBjEqWO6TDexy/2V1/Q0kh/aSjouUO6eOtja+/xbGQ\nvnNaV7jcIT19hDRLMaJZOipd1pCWz+tNd3gvyd3q+tUGIZ2TUfGyhnTUP5ztUoxolI7Kl/M+\n0na7Xi+X/SWH1dWOhHRGRxWws6F0lkdVKCek7lSaETXSUR2yhrRdHfc1zP/66RVC+iCjSuQM\n6fnkJWeZZkRrdFSLjCFtuqfX/f5lsdxv1/Pu6rZVIR3pqBoZQ1ocdwdtu+e3nK6/JAnpwPKo\nItl3f79varBF6E86qknGkD73q57uuYsd0RIZVSVjSKtu8bLfvy67p/3u6e2XBCMaoqO65Lxq\n9/49fbPdYYvQa5IRrXBaV5us95HWbynNn/e2CP1FR9UpZ2dD5hElk1F9hFQeHVVISMXRUY2E\nVBjLozoJqSw6qpSQiiKjWgmpJDqqlpAKoqN6CakYlkc1E1IpdFQ1IRVCRnUTUhl0VDkhlcBp\nXfWEVAAd1U9I45NRA4Q0Oh21QEhj01EThDQuy6NGCGlUOmqFkMYko2YIaUQ6aoeQxqOjhghp\nLJZHTRHSSHTUFiGNQ0aNEdIodNQaIY3AaV17hJSfjhokpOxk1CIh5aajJgkpMx21SUhZWR61\nSkg56ahZQspIRu0SUj46apiQstFRy4SUieVR24SUh44aJ6QsZNQ6IeWgo+YJKT2ndRMgpOR0\nNAVCSk1GkyCkxHQ0DUJKS0cTIaSULI8mQ0gJ6Wg6hJSOjCZESMnoaEqElIqOJkVIaVgeTYyQ\nktDR1AgpBRlNjpAS0NH0CCmc07opElI0HU2SkILJaJqEFEtHEyWkUDqaKiEFsjyaLiHF0dGE\nCSmMjKZMSFF0NGlCCqKjaRNSCMujqRNSBB1NnpACyAghDacjhDSY0zr2QhpMRxwIaRgZ0RPS\nIDriSEhD6Ih3Qnqc5RGfhPQwHfFFSI+SESeE9CAdcSprSC/Py+5guXpJNSIXHXEmY0i7efdl\nkWRELpZHfJMxpFU3+7ftH71uZt0qxYhMdMR3GUOaddvPx9tulmJEHjLih4whdd1v/xI2Igsd\n8ZNXpDs5reOSvGukzWv/qOI1ko64KOfl78XJVbv5LsmI1GTEZXnvI636+0iz5XOl95F0xC/s\nbLiDjvhNOSF1p9KMGMbyiN9lDWm7Oi6T5st/qUakoyOuyBnS88lLzjLNiHRkxDUZQ9p0T6/7\n/ctiud+u590mxYh0dMRVGUNadP0l7233/JbT9Zek4kLSEdeNsEWo39RQ1RYhyyP+knWLUP+K\ntOsbqikkHfGnrFuEFi/7/euye9rvnt5+STBioMvX3WXE30bYIjTbvT1jZ69JRgzxyyuljrhB\n1vtI67eU5s9vD2arq1vtRgrp0minddyknJ0NmUf8OvNsto64jZC+zzydLSNuJKTvM09m64hb\nCenbUB3xCCF9DT2/amd5xB2EdDr25OK3jriHkC6TEXcR0kU64j5CukRH3ElIP1kecTch/aAj\n7iek72TEA4T0jY54hJDOOK3jMUI6pSMeJKQTMuJRQvqiIx4mpE864nFCemd5xBBCOtIRgwip\nJyOGEdKBjhhISHsdMZyQLI8IICQdEWDyIcmICFMPSUeEmHZITusIMumQdESUKYckI8JMOCQd\nEWe6IemIQFMNyfKIUBMNSUfEmmZIMiLYJEPSEdGmGJKOCDe9kCyPSGByIemIFKYWkoxIYmIh\n6Yg0JhWS0zpSmVJIOiKZCYUkI9KZTkg6IqHJhKQjUppISJZHpDWNkHREYpMISUakNoWQdERy\nEwhJR6TXfEiWR+TQekg6IovGQ5IRebQdko7IpOWQnNaRTcMh6Yh82g1JRmTUbEg6IqdWQ9IR\nWbUZkuURmTUZko7IrcWQZER2DYakI/JrLyQdMYLWQrI8YhSNhaQjxtFWSDJiJE2FpCPG0lBI\nTusYTzsh6YgRNROSjBhTKyHpiFE1EpKOGFcTIVkeMbYWQtIRo2sgJBkxvvwhreddt9zEjdAR\nBcgYUtf/HxddbxU1QkeUIHdIq2612+9fV906ZITlEWXIHdKs2x0e77p5xAgdUYjcIXXdyb8M\nHSEjSpE7pKePkGbDR+iIYmQNafm83nT/3h7uVtevNtwywmkdBcka0lH/cLYbOEJHlCTnfaTt\ndr1eLvtLDqurHd0wQkYUpdKdDTqiLOWE1J3643+rIwqTM6TXp2723O8Rml3f2PDHCMsjipMx\npN3s8Fqzfu5fchaPj9AR5ckYUn/JezXrnnaDLn/LiAJlDGnW/x+74x6hh2/I6ogSZd/9/X4h\n4dEtQjqiSCO8Ih1+3T32imR5RKFGWCMdbsY+tkbSEaWq6aqdjChWRfeRdES5ytnZ8McIp3WU\nrJaQdETRKglJRpStjpB0ROGqCElHlK6CkCyPKF/5IemIChQfkoyoQekh6YgqFB6SjqhD0SFZ\nHlGLkkPSEdUoOCQZUY9yQ9IRFSk1JKd1VKXQkHREXcoMSUZUptSQ0s+AQGWGlGEERBISBBAS\nBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBKgrpO76D0OHsdQUUl+RlChR\nVSHlGg/3qiik7tofwqiEBAGEBAEqCskaiXJVFZKrdpSqppDcR6JYdYUEhRISBBASBBASBBAS\nBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBCg0JKjMA8/y+HAeMsrH\nMcZQB9rYzPFHn/H8am3oVGaOP/qM51drQ6cyc/zRZzy/Whs6lZnjjz7j+dXa0KnMHH/0Gc+v\n1oZOZeb4o894frU2dCozxx99xvOrtaFTmTn+6DOeX60NncrM8Uef8fxqbehUZo4/+oznV2tD\npzJz/NFnPL9aGzqVmeOPhnYICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQII\nCQIICQIICQKMHNJuNetmq13/kTz69uX3exu62Hw+fJ+fbWamA11/TDg5xNRHe2Fm8qNddz8f\nZvuqfhk3pNdZ/1meve7323whLfpBz18P5zlnZjrQ7ceEk0NMfbQXZiY/2u3Xf/vSIWczbkhP\n3ert11X3dPgsLDMNXXeL3X731G33+5dutt1vZ91Lxpl5DvTtoI5f2pNDTH20l2amPtrPmZfH\n5zNuSO9HfvjH+vgSkcGi/wy/HhpedYezrX/pR5/MzHKgb+G+f25PDjHx0V6cmfhov2ZeHp/P\nuCHN3kOaHT4P60xDP+pd7PfL7u2kMsdrxMnMLAf6Vuz7yJNDTHy0F2cmPtqvmZfH5zNuSM/v\np3bPh6PfPL2tEDMMPXkZPHmYbWaWA91+P7YMR3txZuKj/Zp5eXw+I1+1Wx+uNswOf2ktj4vS\nRfqZ8/4vrJesIZ3MzHWg2UO6ODP90Z4czoRDev68mNV1/w5XwzOc9zx3y91+u8ga0tnMPAda\nRkjpj1ZI+8MZ9NuL/u7p69O8y3HRsr/mvswa0snMo/QHWkZIRymPVkj7wxnP4a7Z6ac5x9G/\nlTt77ifNsn3Kv2a+Sz7zfcDJIaY/2p8zz/8g4cyr4zMo5vL3+e+ktz3Ue7y+85rr+s42598Y\nZ5ewXr+u2qU82p8zz/8g4cyr4zMo4fL37nD5e9a/OOU4+uOk9WHSc3/HYdMlv1h4MjPXgb4/\nq04OMf3R/pyZ/mh/hpTtq3r2YeQc9sOqO+yIWh3vja76Vekmw9Cn/f5lflgFZ7sHfjIz14Hm\n39lwaWb6o/0Z0gR3NrzvijpcHN0dt91l+FvkfVL/d+Q80zX3k5m5DvTjCXZyiMmP9ufM9Ed7\nYV2Q66t69mFknfZTv0+3f3TYCD7Psrvh9entKb35HJrlLvC3mRkO9ONZtTv/FCc92l9mJj3a\nCyFl+6qefhh5x0GbhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQB\nhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQB\nhNSA9D93l78IqX5zX8Tx+RrUr/NFHJ+vQf2EVABfg8qs593s+EPCN4uuW2wOHXVSGp2vQF2W\nfTaLt0fr/lG3FlIRfAWqsukWu/1u0b29EM267X7/r5s7tSuCr0FVlt3u7dddtzzk83HVW0gF\n8DWoSvdhv1913XK7Pf7m2B8WQqrLSUj759nbP2evQiqCr0FVzpvZrObWSIXwNajKsvu2HegQ\nkZAK4GtQlX/dbHu48r08bAz693nV7nXsjwsh1WXRr5AOK6N/x8XSyyGpbjb2xzV5QqrM+i2b\np/4VqN/Z8NbR/mUupNEJCQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQII\nCQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQL8\nD57NCuHFe9kMAAAAAElFTkSuQmCC", "text/plain": [ "plot without title" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot(est, est2)\n", "model = lm(est~est2)\n", "summary(model)\n", "abline(model, lwd= 2, col = \"red\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "8. Genetyka" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 'AB'
    2. \n", "\t
  2. 'BB'
    3. \n", "\t
  3. 'AA'
    4. \n", "\t
  4. 'AA'
    5. \n", "\t
  5. 'BB'
    6. \n", "\t
  6. 'BB'
    7. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 'AB'\n", "\\item 'BB'\n", "\\item 'AA'\n", "\\item 'AA'\n", "\\item 'BB'\n", "\\item 'BB'\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 'AB'\n", "2. 'BB'\n", "3. 'AA'\n", "4. 'AA'\n", "5. 'BB'\n", "6. 'BB'\n", "\n", "\n" ], "text/plain": [ "[1] \"AB\" \"BB\" \"AA\" \"AA\" \"BB\" \"BB\"" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Genotypes = sample(c(\"AA\", \"AB\", \"BB\"), 5000, replace = TRUE)\n", "head(Genotypes)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
V1V2V3V4V5V6V7V8V9V10...V991V992V993V994V995V996V997V998V999V1000
AB BB AA BB BB AA AB BB AA AA ...BB BB BB AB AB BB AB AB BB AA
BB BB AA AB AB AA BB BB BB BB ...BB AA AA AA AB AB BB AB AB AB
AA AA AA AA BB AB AB AB BB AB ...BB AA BB AA AB AB AA AA BB AB
AA AB BB AB BB AA BB AA AA BB ...BB BB AA BB AA AA AB AB BB AB
BB AA AB AA AB AB AA AB AA BB ...AA AB BB AB AA AB BB BB AA BB
\n" ], "text/latex": [ "\\begin{tabular}{r|llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll}\n", " V1 & V2 & V3 & V4 & V5 & V6 & V7 & V8 & V9 & V10 & ... & V991 & V992 & V993 & V994 & V995 & V996 & V997 & V998 & V999 & V1000\\\\\n", "\\hline\n", "\t AB & BB & AA & BB & BB & AA & AB & BB & AA & AA & ... & BB & BB & BB & AB & AB & BB & AB & AB & BB & AA \\\\\n", "\t BB & BB & AA & AB & AB & AA & BB & BB & BB & BB & ... & BB & AA & AA & AA & AB & AB & BB & AB & AB & AB \\\\\n", "\t AA & AA & AA & AA & BB & AB & AB & AB & BB & AB & ... & BB & AA & BB & AA & AB & AB & AA & AA & BB & AB \\\\\n", "\t AA & AB & BB & AB & BB & AA & BB & AA & AA & BB & ... & BB & BB & AA & BB & AA & AA & AB & AB & BB & AB \\\\\n", "\t BB & AA & AB & AA & AB & AB & AA & AB & AA & BB & ... & AA & AB & BB & AB & AA & AB & BB & BB & AA & BB \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| V1 | V2 | V3 | V4 | V5 | V6 | V7 | V8 | V9 | V10 | ... | V991 | V992 | V993 | V994 | V995 | V996 | V997 | V998 | V999 | V1000 |\n", "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n", "| AB | BB | AA | BB | BB | AA | AB | BB | AA | AA | ... | BB | BB | BB | AB | AB | BB | AB | AB | BB | AA |\n", "| BB | BB | AA | AB | AB | AA | BB | BB | BB | BB | ... | BB | AA | AA | AA | AB | AB | BB | AB | AB | AB |\n", "| AA | AA | AA | AA | BB | AB | AB | AB | BB | AB | ... | BB | AA | BB | AA | AB | AB | AA | AA | BB | AB |\n", "| AA | AB | BB | AB | BB | AA | BB | AA | AA | BB | ... | BB | BB | AA | BB | AA | AA | AB | AB | BB | AB |\n", "| BB | AA | AB | AA | AB | AB | AA | AB | AA | BB | ... | AA | AB | BB | AB | AA | AB | BB | BB | AA | BB |\n", "\n" ], "text/plain": [ " V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 ... V991 V992 V993 V994 V995 V996 V997 V998\n", "1 AB BB AA BB BB AA AB BB AA AA ... BB BB BB AB AB BB AB AB \n", "2 BB BB AA AB AB AA BB BB BB BB ... BB AA AA AA AB AB BB AB \n", "3 AA AA AA AA BB AB AB AB BB AB ... BB AA BB AA AB AB AA AA \n", "4 AA AB BB AB BB AA BB AA AA BB ... BB BB AA BB AA AA AB AB \n", "5 BB AA AB AA AB AB AA AB AA BB ... AA AB BB AB AA AB BB BB \n", " V999 V1000\n", "1 BB AA \n", "2 AB AB \n", "3 BB AB \n", "4 BB AB \n", "5 AA BB " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Genotypes = as.data.frame(matrix(Genotypes, 5, 1000))\n", "head(Genotypes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathbf{G} = \\frac{\\textbf{MM}^{T}}{2\\sum_{i=1}^{N_{SNP}}p_{i}(1-p_{i})}$, gdzie $\\textbf{M}_{ij}\\in\\{2-2p_{i}, 1-2p_{i}, -2p_{i}\\}$ dla genotypu $SNP_{i} = \\{AA, AB, BB\\}$ i osobnika $j$." ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "data": { "text/html": [ "5000" ], "text/latex": [ "5000" ], "text/markdown": [ "5000" ], "text/plain": [ "[1] 5000" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
2 1 2 2 0 0 2 2 0 2 ...0 1 0 1 2 0 1 1 1 0
2 0 1 2 2 0 2 1 1 1 ...1 1 1 1 0 1 2 0 1 2
2 1 1 1 1 0 2 2 0 2 ...1 0 0 1 0 0 2 1 1 1
0 2 2 0 1 1 2 0 0 0 ...2 2 0 1 2 0 0 1 1 1
1 0 2 2 2 2 0 0 1 1 ...2 2 0 1 1 1 2 0 0 2
\n" ], "text/latex": [ "\\begin{tabular}{llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll}\n", "\t 2 & 1 & 2 & 2 & 0 & 0 & 2 & 2 & 0 & 2 & ... & 0 & 1 & 0 & 1 & 2 & 0 & 1 & 1 & 1 & 0 \\\\\n", "\t 2 & 0 & 1 & 2 & 2 & 0 & 2 & 1 & 1 & 1 & ... & 1 & 1 & 1 & 1 & 0 & 1 & 2 & 0 & 1 & 2 \\\\\n", "\t 2 & 1 & 1 & 1 & 1 & 0 & 2 & 2 & 0 & 2 & ... & 1 & 0 & 0 & 1 & 0 & 0 & 2 & 1 & 1 & 1 \\\\\n", "\t 0 & 2 & 2 & 0 & 1 & 1 & 2 & 0 & 0 & 0 & ... & 2 & 2 & 0 & 1 & 2 & 0 & 0 & 1 & 1 & 1 \\\\\n", "\t 1 & 0 & 2 & 2 & 2 & 2 & 0 & 0 & 1 & 1 & ... & 2 & 2 & 0 & 1 & 1 & 1 & 2 & 0 & 0 & 2 \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 2 | 1 | 2 | 2 | 0 | 0 | 2 | 2 | 0 | 2 | ... | 0 | 1 | 0 | 1 | 2 | 0 | 1 | 1 | 1 | 0 |\n", "| 2 | 0 | 1 | 2 | 2 | 0 | 2 | 1 | 1 | 1 | ... | 1 | 1 | 1 | 1 | 0 | 1 | 2 | 0 | 1 | 2 |\n", "| 2 | 1 | 1 | 1 | 1 | 0 | 2 | 2 | 0 | 2 | ... | 1 | 0 | 0 | 1 | 0 | 0 | 2 | 1 | 1 | 1 |\n", "| 0 | 2 | 2 | 0 | 1 | 1 | 2 | 0 | 0 | 0 | ... | 2 | 2 | 0 | 1 | 2 | 0 | 0 | 1 | 1 | 1 |\n", "| 1 | 0 | 2 | 2 | 2 | 2 | 0 | 0 | 1 | 1 | ... | 2 | 2 | 0 | 1 | 1 | 1 | 2 | 0 | 0 | 2 |\n", "\n" ], "text/plain": [ " [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]\n", "[1,] 2 1 2 2 0 0 2 2 0 2 ... 0 1 0 \n", "[2,] 2 0 1 2 2 0 2 1 1 1 ... 1 1 1 \n", "[3,] 2 1 1 1 1 0 2 2 0 2 ... 1 0 0 \n", "[4,] 0 2 2 0 1 1 2 0 0 0 ... 2 2 0 \n", "[5,] 1 0 2 2 2 2 0 0 1 1 ... 2 2 0 \n", " [,15] [,16] [,17] [,18] [,19] [,20] [,21]\n", "[1,] 1 2 0 1 1 1 0 \n", "[2,] 1 0 1 2 0 1 2 \n", "[3,] 1 0 0 2 1 1 1 \n", "[4,] 1 2 0 0 1 1 1 \n", "[5,] 1 1 1 2 0 0 2 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Genotypes2 = sample(c(2, 1, 0), 5000, replace = TRUE)\n", "length(Genotypes2)\n", "Genotypes2 = matrix(Genotypes2, 5, 1000)\n", "head(Genotypes2)" ] }, { "cell_type": "code", "execution_count": 70, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Initial data: \n", "\tNumber of Individuals: 5 \n", "\tNumber of Markers: 988 \n", "\n", "Missing data check: \n", "\tTotal SNPs: 988 \n", "\t 0 SNPs dropped due to missing data threshold of 1 \n", "\tTotal of: 988 SNPs \n", "MAF check: \n", "\t 5 SNPs dropped with MAF below 0.01 \n", "\tTotal: 983 SNPs \n", "Monomorphic check: \n", "\t 6 monomorphic SNPs \n", "\tTotal: 977 SNPs \n", "Summary check: \n", "\tInitial: 988 SNPs \n", "\tFinal: 977 SNPs ( 11 SNPs removed) \n", " \n", "Completed! Time = 0.06 seconds \n" ] }, { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
1.2357915-0.3175200-0.3494598-0.2747769-0.2940348
-0.3175200 1.2245186-0.2611555-0.3297323-0.3161109
-0.3494598-0.2611555 1.2287459-0.3006106-0.3175200
-0.2747769-0.3297323-0.3006106 1.2325035-0.3273837
-0.2940348-0.3161109-0.3175200-0.3273837 1.2550493
\n" ], "text/latex": [ "\\begin{tabular}{lllll}\n", "\t 1.2357915 & -0.3175200 & -0.3494598 & -0.2747769 & -0.2940348\\\\\n", "\t -0.3175200 & 1.2245186 & -0.2611555 & -0.3297323 & -0.3161109\\\\\n", "\t -0.3494598 & -0.2611555 & 1.2287459 & -0.3006106 & -0.3175200\\\\\n", "\t -0.2747769 & -0.3297323 & -0.3006106 & 1.2325035 & -0.3273837\\\\\n", "\t -0.2940348 & -0.3161109 & -0.3175200 & -0.3273837 & 1.2550493\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 1.2357915 | -0.3175200 | -0.3494598 | -0.2747769 | -0.2940348 |\n", "| -0.3175200 | 1.2245186 | -0.2611555 | -0.3297323 | -0.3161109 |\n", "| -0.3494598 | -0.2611555 | 1.2287459 | -0.3006106 | -0.3175200 |\n", "| -0.2747769 | -0.3297323 | -0.3006106 | 1.2325035 | -0.3273837 |\n", "| -0.2940348 | -0.3161109 | -0.3175200 | -0.3273837 | 1.2550493 |\n", "\n" ], "text/plain": [ " [,1] [,2] [,3] [,4] [,5] \n", "[1,] 1.2357915 -0.3175200 -0.3494598 -0.2747769 -0.2940348\n", "[2,] -0.3175200 1.2245186 -0.2611555 -0.3297323 -0.3161109\n", "[3,] -0.3494598 -0.2611555 1.2287459 -0.3006106 -0.3175200\n", "[4,] -0.2747769 -0.3297323 -0.3006106 1.2325035 -0.3273837\n", "[5,] -0.2940348 -0.3161109 -0.3175200 -0.3273837 1.2550493" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Gmatrix(\n", " SNPmatrix = Genotypes2,\n", " method = \"VanRaden\",\n", " missingValue = -9,\n", " maf = 0.01)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "9. Wrzucenie macierzy G do modelu" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathbf{y} = \\mathbf{X\\beta+Z_{1}a+Z_{2}g+\\epsilon}$, $\\mathbf{g}\\sim\\mathcal{N}\\left(0,\\mathbf{G}\\sigma^{2}_{g}\\right)$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$ \\left[ \\begin{array}{c}\n", " \\widehat{\\mathbf{\\beta}} \\\\\n", " \\widehat{\\mathbf{a}} \\\\ \n", " \\widehat{\\mathbf{g}} \\end{array} \\right]=\n", " \\left[ \\begin{array}{cc}\n", " \\mathbf{X}^{T}\\mathbf{X} & \\mathbf{X}^{T}\\mathbf{Z}_{1} & \\mathbf{X}^{T}\\mathbf{Z}_{2} \\\\\n", " \\mathbf{Z}_{1}^{T}\\mathbf{X} & \\mathbf{Z}_{1}^{T}\\mathbf{Z}_{1}+\\mathbf{A}^{-1}\\alpha_{1} & \\mathbf{Z}_{1}^{T}\\mathbf{Z}_{2} \\\\ \n", " \\mathbf{Z}_{2}^{T}\\mathbf{X} & \\mathbf{Z}_{2}^{T}\\mathbf{Z}_{1} & \\mathbf{Z}_{2}^{T}\\mathbf{Z}_{2}+\\mathbf{G}^{-1}\\alpha_{2}\\end{array}\\right]^{-1}\\cdot\n", " \\left[ \\begin{array}{c}\n", " \\mathbf{X}^{T}\\mathbf{y} \\\\\n", " \\mathbf{Z}_{1}^{T}\\mathbf{y} \\\\\n", " \\mathbf{Z}_{2}^{T}\\mathbf{y}\n", " \\end{array}\\right]$," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\alpha_{1} = \\frac{\\sigma^{2}_{\\epsilon}}{\\sigma^{2}_{a}}$ i $\\alpha_{2} = \\frac{\\sigma^{2}_{\\epsilon}}{\\sigma^{2}_{g}}$" ] }, { "cell_type": "code", "execution_count": 72, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
4.3584837470
3.4044246438
-0.0086571956
-0.1856393783
0.1767871430
-0.2493226696
0.1825171246
0.0001067897
-0.0002298413
0.0002327889
-0.0004635470
0.0003682570
\n" ], "text/latex": [ "\\begin{tabular}{l}\n", "\t 4.3584837470\\\\\n", "\t 3.4044246438\\\\\n", "\t -0.0086571956\\\\\n", "\t -0.1856393783\\\\\n", "\t 0.1767871430\\\\\n", "\t -0.2493226696\\\\\n", "\t 0.1825171246\\\\\n", "\t 0.0001067897\\\\\n", "\t -0.0002298413\\\\\n", "\t 0.0002327889\\\\\n", "\t -0.0004635470\\\\\n", "\t 0.0003682570\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 4.3584837470 |\n", "| 3.4044246438 |\n", "| -0.0086571956 |\n", "| -0.1856393783 |\n", "| 0.1767871430 |\n", "| -0.2493226696 |\n", "| 0.1825171246 |\n", "| 0.0001067897 |\n", "| -0.0002298413 |\n", "| 0.0002327889 |\n", "| -0.0004635470 |\n", "| 0.0003682570 |\n", "\n" ], "text/plain": [ " [,1] \n", " [1,] 4.3584837470\n", " [2,] 3.4044246438\n", " [3,] -0.0086571956\n", " [4,] -0.1856393783\n", " [5,] 0.1767871430\n", " [6,] -0.2493226696\n", " [7,] 0.1825171246\n", " [8,] 0.0001067897\n", " [9,] -0.0002298413\n", "[10,] 0.0002327889\n", "[11,] -0.0004635470\n", "[12,] 0.0003682570" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "library(MASS)\n", "\n", "mme3 = function(y, X, Z1, Z2, A, G, sigma_a, sigma_e, n) {\n", " sigma_g = sigma_a / n\n", " alpha1 = sigma_e / sigma_a\n", " alpha2 = sigma_e / sigma_g\n", " invA = ginv(A)\n", " invG = ginv(G)\n", " C = rbind(cbind(t(X)%*%X, t(X)%*%Z1, t(X)%*%Z2),\n", " cbind(t(Z1)%*%X, t(Z1)%*%Z1+invA*c(alpha1), t(Z1)%*%Z2),\n", " cbind(t(Z2)%*%X, t(Z2)%*%Z1, t(Z2)%*%Z2 + invG*c(alpha2)))\n", " rhs = rbind(t(X)%*%y, t(Z1)%*%y, t(Z2)%*%y)\n", " invC = ginv(C)\n", " estimators = invC%*%rhs\n", " list(C = C, est = estimators)\n", "}\n", "\n", "(est = mme3(y, X, diag(5), diag(5), A[4:8, 4:8], G, 20, 40, 988)$est)" ] }, { "cell_type": "code", "execution_count": 73, "metadata": {}, "outputs": [ { "data": { "text/html": [ "12" ], "text/latex": [ "12" ], "text/markdown": [ "12" ], "text/plain": [ "[1] 12" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "length(est)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "10. Estymacja efektów SNP" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\widehat{SNP} = \\left(\\left(\\frac{1-k}{N_{SNP}}\\right)^{-1} + \\frac{1}{k}\\textbf{MA}^{-1}\\textbf{M}^{T} \\right)^{-1}\\frac{1}{k}\\textbf{MA}^{-1}\\widehat{g}$" ] }, { "cell_type": "code", "execution_count": 74, "metadata": {}, "outputs": [], "source": [ "k = 0.2\n", "snp = ginv((t(M)%*%ginv(A)[4:8, 4:8]%*%M)/k + 1 / ((1-k)/length(p2)) ) / k " ] }, { "cell_type": "code", "execution_count": 75, "metadata": {}, "outputs": [], "source": [ "snp = snp%*%t(M)%*%ginv(A)[4:8, 4:8]" ] }, { "cell_type": "code", "execution_count": 76, "metadata": {}, "outputs": [], "source": [ "snp = snp %*% est[8:12]" ] }, { "cell_type": "code", "execution_count": 77, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " V1 \n", " Min. :-2.154e-06 \n", " 1st Qu.:-6.785e-07 \n", " Median : 1.567e-08 \n", " Mean : 0.000e+00 \n", " 3rd Qu.: 6.004e-07 \n", " Max. : 2.207e-06 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "summary(snp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "11. Istotność markerów SNP" ] }, { "cell_type": "code", "execution_count": 78, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "8.55663583183359e-07" ], "text/latex": [ "8.55663583183359e-07" ], "text/markdown": [ "8.55663583183359e-07" ], "text/plain": [ "[1] 8.556636e-07" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "-6.0666215783731e-22" ], "text/latex": [ "-6.0666215783731e-22" ], "text/markdown": [ "-6.0666215783731e-22" ], "text/plain": [ "[1] -6.066622e-22" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "(s = sd(snp))\n", "(m= mean(snp))" ] }, { "cell_type": "code", "execution_count": 79, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "1" ], "text/latex": [ "1" ], "text/markdown": [ "1" ], "text/plain": [ "[1] 1" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "8.28603979405942e-18" ], "text/latex": [ "8.28603979405942e-18" ], "text/markdown": [ "8.28603979405942e-18" ], "text/plain": [ "[1] 8.28604e-18" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "snp2 = (snp - m) / s\n", "sd(snp2)\n", "mean(snp2)" ] }, { "cell_type": "code", "execution_count": 80, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ " V1 \n", " Min. :-2.51783 \n", " 1st Qu.:-0.79293 \n", " Median : 0.01831 \n", " Mean : 0.00000 \n", " 3rd Qu.: 0.70167 \n", " Max. : 2.57984 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "summary(snp2)" ] }, { "cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [], "source": [ "W = snp2\n", "p.value = 2*pnorm(abs(W), lower.tail = FALSE)" ] }, { "cell_type": "code", "execution_count": 84, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAMFBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////QFLu4AAAACXBIWXMAABJ0AAAS\ndAHeZh94AAAgAElEQVR4nO2diXqrug6FnTZNuzvx/m+7mwGs0dggICTr/+49ScGWZFnLJmTY\nqQMAzCZtHQAAjwCEBEAAEBIAAUBIAAQAIQEQAIQEQAAQEgABQEgABAAhARAAhARAABASAAFA\nSAAEACEBEACEBEAAEBIAAUBIAAQAIQEQAIQEQAAQEgABQEgABAAhARAAhARAABASAAFASAAE\nACEBEACEBEAAEBIAAUBIAAQAIQEQAIQEQAAQEgABQEgABAAhARAAhARAABASAAFASAAEACEB\nEACEBEAAEBIAAUBIAAQAIQEQAIQEQAAQEgABQEgABAAhARAAhARAABASAAFASAAEACEBEACE\nBEAAEBIAAUBIAAQAIQEQAIQEQAAQEgABQEgABAAhARAAhARAABASAAFASAAEACEBEACEBEAA\nEBIAAUBIAAQAIQEQAIQEQAAQEgABQEgABAAhARAAhARAABASAAFASAAEACEBEACEBEAAEBIA\nAUBIAAQAIQEQAIQEQAAQEgABQEgABAAhARAAhARAABASAAFASAAEACEBEACEBEAAEBIAAUBI\nAAQAIQEQAIQEQAAQEgABQEgABAAhARAAhARAABASAAFASAAEACEBEMAKQkoA7IwJVR4vnA1c\nABAJhARAABASAAFASAAEACEBEACEBEAAEBIAAUBIAAQAIQEQwKpC+no/Xt4EPp6+lnIBwCas\nKKTfF/KBitdFXACwESsK6ZQO/74vz34+D+m0hAsANmJFIR3S9/D8Ox2WcAHARqwoJPYB2fKn\nZSEksDOwIwEQwLqvkT5/Ls/wGgk8Gmve/n4ld+1efhdxAcA2rPs+0unyPtLh+I73kcBjgU82\nABAAhARAABASAAFsJSS8jwQeivsR0szfNgJgS3BpB0AAEBIAAUBIAASwgZA+DunlY1kXAKzM\nmkL6PqbDR/eOL/aBx2NFIX1fFHRKb7/dzzEV9yQICeyMFYX0dv7E9+n6/Ynf9LKECwA2YvUv\n9qUj+SPaBQAbsbqQ/l2v6fDFPvBQrHpp99Z/Cen3DV/sAw/Fmj/HdRiu51J5Q4KQwN5Y9X2k\nUy+fQ3E/gpDA7sAnGwAIAEICIAAICYAAICQAAoCQAAgAQgIgAAgJgAAgJAACgJAACABCAiAA\nCAmAACAkAAKAkAAIAEICIAAICYAAICQAAoCQAAgAQgIgAAgJgAAgJAACgJAACABCAiAACAmA\nACAkAAKAkAAIAEICIAAICYAAICQAAoCQAAgAQgIgAAgJgAAgJAACgJAACABCAiAACAmAACAk\nAAKAkAAIAEICIAAICYAAICQAAoCQAAgAQgIgAAgJgAAgJAACgJAACABCAiAACAmAACAkAAKA\nkAAIAEICIAAICYAAICQAAoCQAAgAQgIgAAgJgAAgJAACgJAACABCAiAACAmAACAkAAKAkAAI\nAEICIAAICYAAICQAAoCQAAgAQgIgAAgJgAAgJAACgJAACGB9IX28pHT8XNQFAGuzopDSpeNr\nunBaxAUAG7G2kE7p9Nt1P6f0sYQLADZibSEd0u/5+W96WcIFABuxtpBSIn+EuwBgI9YW0lsv\npMMSLgDYiFWFdHz/+Ez//p7+nsp3GyAksDNWFdKVy9PD7xIuANiINd9H+v7++DgeL7ccTkUd\nQUhgb+CTDQAEACEBEACEBEAAWwkJ7yOBh+J+hJQoES4AWA9c2gEQAIQEQAAQEgABrCmk37eU\nXm9f6cPNBvBQrCik38PlPsLxagRCAo/EikK6fJnv9+PwejECIYFHYkUhHa4dfw4vPxASeDBW\n/82Gv03p9RVCAg/GikJ6Sf1Hvl9eISTwWKwopI/0dnv2k14hJPBQrHn7+zSo53PkU0AQEtgZ\nq74h+33sn/28QUjgkcAnGwAIAEICIAAICYAAICQAAoCQAAgAQgIgAAgJgAAgJAACgJAACABC\nAiAACAmAACAkAAKAkAAIAEICIAAICYAAICQAAoCQAAgAQgIgAAgJgAAgJAACgJAACABCAiAA\nCAmAACAkAAKAkAAIAEICIAAICYAAICQAAoCQAAgAQgIgAAgJgAAgJAACgJAACABCAiAACAmA\nACAkAAKAkAAIAEICIAAICYAAICQAAoCQAAgAQgIgAAgJgAAgJAACgJAACABCAiAACAmAACAk\nAAKAkAAIAEICIAAICYAAICQAAoCQAAgAQgIgAAgJgAAgJAACgJAACABCAiAACAmAACAkAAKA\nkAAIAEICIAAICYAAICQAAoCQAAgAQgIgAAgJgAAgJAACgJAACABCAiCAVYX09X5MZ46nr6Vc\nALAJKwrp9yVlXhdxAcBGrCikUzr8+748+/k8pNMSLgDYiBWFdEjfw/PvdFjCBQAbsaKQUvL+\nCHMBwEZgRwIggHVfI33+XJ7hNRJ4NNa8/f1K7tq9/C7iAoBtWPd9pNPlfaTD8X3D95HOESxo\nHjwnT/fJhtuOuJwD8JQ8m5BuGoKSQCzPJyT6AEAQWwlpo/eRercQEojlfoSUKBEuSn4hJBDL\n813a4TUSWIBnExLu2oFFeDoh4X0ksARrCun3LaXXz5sRfGgVPBJrfrHvcP167NUIhAQeiVU/\ntPrxp6aPw+XLsRASeChW/RrF5eHn8PIDIYEHY4Mv9v2+vkJI4MFYUUgvqf/qxMsrhAQeixWF\n9JHebs9+0iuEtBDq7j49sM2t/6d4w2HN29+nIZ+fI6l9/LwvhXq/mR7Y5s3o53gLfNU3ZL+P\n/bOfNwhpCdQnoOiBbT4e9SQfynq+TzY8NOozufTANh/YfZKPCUNIj4T6lgg9sM1XSJ7liysQ\n0kOBHWkrIKSHAq+RtgJCeixw124jIKQHA+8jbQOEBEAAEBIAAUBIAAQAIQEQAIQEQABzhfR1\nOv8bE69j/7ryHBcA3D/zhPQv//vKL59xQUFIYG/MEdLPa3r9+D5/W+/36/3v+c+WUQGwJTOE\n9JlO9F8L+zmlsE0JQgI7Y4aQjvIf3ft9k02nAiGBnYG7dgAEACEBEMAsIX29H68/nhp89xtC\nAntjhpB+873vlF63jgqALZkhpFM6/Pu+PPv5PKRTXEwQEtgdM4R0SN/Dse90iImHuwBgJ8wQ\nEvu2VuxXtyAksDOwIwEQwLzXSJ/XTwXhNRJ4dubc/n4ld+1e5Mcc1o4KgC2Z9z7S6fI+0uH4\njveRwHODTzYAEACEBEAAEBIAATyMkCp+hXDaDxU+xc8bzgVJehQhVfwu7rSfzn2OH9ydCZI0\n75MNjE2jqvil9mk/5v4kPwE/DyRplpA+7kpI9GFqk7BeTwaSNO/S7vsQ++UJw0Vt+8QfpzUJ\n6/VkIEndzNdI37EfDLJc1HbAjrQhSNLcmw0f5HOrkeA10q5AknDXLsIwQJIeRkh4H2lTkKSH\nERIAWwIhARAAhARAABASAAFASAAEACHFg1tYMVzyKJI5/GknebvUQ0jh4E2VGKwPcg5/2kne\nMPUQUjR4mz+GQS8d/TTfLbd2krdMfZSQ/sYQ9q8j7V1I9AFM5aYgIaTbg53kLVMfJ6Tu33Fu\nMGUX+wAfhY6hFxAVktiYZJI3TT0u7aLBjhTDk+5Isey5CPEaKYYnfY0Uy66LEHftYniyu3af\nx3Pgx5+geCwXewMyiuGp3kd6vS0Qh1AloQ7BzpgppI/0+nsW0kfkzW8ICeyOmUI6pF92gzII\nCAnsjJlCur7u6yAk8OTMFNLLbUf6Ti9hIXUQEtgdMa+RPg/pIyykDkICu2PuXbvj7dZ97E9F\nQkhgZ4S8j5SO/4LCMV0AcPfgkw0ABAAhARAAhARAALPfR1KfK4wAQgI7A0ICIICYS7uv17hv\nxzouALhngl4j/eJDq+CpCfzNhkAgJLAzgoT0kQ613T9eUjp+trsA4I4Ju9nwPt7v0vH12rz8\nb2ZCSGBnBAnppeIzqxchndLpt+t+TuUPuUJIYGes+IbsRUjnbwJ255sTxa9dQEhgZ6wtpLof\n8YOQwM6YIaTEGe93bvLWC6l4c6LGmOfxDn7DZ5UQ7mCc7SwY9Mb5WFVIx/ePz3T+xsXvqXy3\nYdSY73PDnzZbNYQ7GGc7Cwa9dT5WvbQbxprS4XeOi+T+pKZ/ZjVWCeEOxtnOgkFvno81P/39\n/f3xcTxebjmcijqqEBJ9qDuzGquEcAfjbGfBoDfPR5SQvkI/bDcSlX/H4g7+LYhVQriDcbaz\nYNDb52OukE6bfPobO9IdjLMd7Ehul6yjkQ/9THdhnp/8GmmFF6QyhEVcNr8mmBhFaPC3BZdY\nDDPfmo/wSZkppEP6172mn5/X9NVoZN77SFPv2m1wQ20hl41mJ0YRG7y6eAk032YqflJmCukc\ny/vfbvTd+ntc1ssbQkX/Ke8jrXRvh4awmMumQpgYRWzwanZDzTcuK3GObybndTmH8nn+3Nw+\nvkaxwZX05hfvc6KIDV5Z2yw3CzieKaTj36XdT3rpvnYhpA3u7Wx/O2lGFLHBK2ub5WYJxzOF\n9HmO5fLNiF18QxY70irdaq1hR8q8n/96G/t+0SwXkXbXeY20scu4KMJfI3Frm+Xm/l4jLcRi\nLpa5hVY0ucqNwlHu6a5djfkgr56Zu7trV/6kz2SWK7sN7kTfgYymRxEbvHn/u9AwyJ9zZq51\nbnBel/Ra/0Zs4gRHtR33cfG2C2pTFZTSNWdmppBe/mI9Vb4X+/GwQqIPoEBtqoJSuubMzH2N\n9PP+p6WX96pLvO9D7bu2e6rJ+7jBvQtqUxWU0lVnJuBmw8/pkOou8b5rb+7tqiaxI1WDHWmk\ny0flS7eP9D3VxWwWe83ffCVuR7LVPYk1/KYhR1aqdATP9xrpzPXqLvTf7Ftg5AvehW40bTdf\nML4J0Szlw/JWe2yW48UJeY10OP1ExWO4iLG45OLUNFnuqmwcXZ41/FIftbtPUP2vtzYF3LV7\na/wGRaOLIIt380LGjmSr+NbwW/ZxPzMzi9nvIwX/M8zaRYzBu7m1ZkeyVXxr+C37uJ+ZmcdM\nIe3mkw33s+5hR1o/ghWIuWsXzc5eIzXx3K+RtolgBaKEFJuHfd21a+SZ79ptFcHyPI2QNnuf\nRvO87yNtF8HiPMulHQCLAiEBEACEBEAAs4X073h+Tzb05yEhJLA75grp9k/CptCf/oaQwN6Y\nKaRTOpw3o89D+d+EneMCgB0wU0iH2/civsv/JuwcFwDsgLmftUvySQgQEtgZsy/t+h1pzX8f\nCYB7Y+7NhvfLa6Sv6l9jmOACgPtn9qUdY8OoANgSCAmAAPDJBgACgJAACGCGkI7y67G/Yf+0\nC4QEdsYMIX2mE5XSzynuH2SGkMDOmHNp9/OaXj++z2L6/Xr/ex73m1yrCMn8MbU1HG+JP8Zn\nGP2VRUY67zXSv5fhht1L5AfA15hRdaPxMb7yXMYf4zOM/soyI517s+HrdP7892vtv0gxxcUy\nqN/ceJAf4Sjij/EZRn9loZE+7V079StQ1T8LtdTCvcKG4I9xhR/FupMN7xZGMZYJoT6rkNTv\nElb/UOFS10ArXFv5Y1zhZxrv5NqxD2Ps11haQ437ZMNr4L/HfMc70lLXQKtcW224I93JteOw\nHxUCmRRq5EeEDu224qJqdjHxNdJSFbfKD45u+BrpTn5Q9bbVjAiJPtQanhAL/eOt/4bsV3es\n/WfEGl0sxLS7dktdA61wbXU1v9Fdu5XGVxPG2EinhTpTSPn7SK/db9y3ZO/4faRd70gbvo90\nNzvSsCsV23QrC4m9Vn+KT3/v+jXShtzJ+GrC2OI1Uv7NhsOTCGnPd+025U7GVxPGBnft8q8I\nnbp/Kexrspunu8RS5XAPZbYkdzK+mjAmhDr3ZkP/u3avZ+9hP8m1xEgBWJC5Quo+z7+0ejxv\nS+k9JiTlwjh/H5cJAPTMFtIibPSmKABT2amQ6MM9spv9cqNAXbe7SZxktpD+nV8lHYP/Seax\nmyr88f7YzaXnRoG6bneTOEXgzYZAdr4j7ebSc6NAXbe7SZxmppA+NvkR/XvP970LfWCjQF23\nu0mcZqaQXrb5Ef37vgK4/0vPGxsF6rrdTeIMIj8iFMfO30fazcKKHSmMsB0p7jsU3aSo7ol7\nv/QcwGukMHb5Gunuue9LTwLu2kWxy7t2989uqgHvIwUx/32k4+rvIwFwd+zykw0A3BsQEgAB\nzBBS4mwcFQBbAiEBEAAu7QAIAEICIIAIIcXf+YeQwM6AkAAIAEICIAAICYAAICQAAoCQAAgA\nt78BCABCAiAACGmHLPmlnV1+IWjBoGtNQ0j7Y8mvke7yK6oLBl1tGkLqXV7TdQ9VRGKwwpn/\nwwb+ICs/glxqsUYGuQ87IVVxjDWqzzWEdPO42CfZp4YintLz9GGeA+fUjH8/aI0MCh9mQqoX\nhGKj+lxDSFeHQ+FeH9b2b4XihDP7x98Kg+zLqlyBpSytkUHhw0xIVRyjjRpyDSFdHSbjYRtI\nDHY4c4Ms9ecLyiQD84KrQfqwfFbFMd7oSXekydcUQ7543uZco0ztSxZBZz2ct+gnOUjb9vjl\njrOl+eeiUD6MhFTF4TXir1G7urncrZCMwc24PDd3pDmX+3bfqhkZ25GC4pr+ImmzHYlHTl8k\nqZBvV8ajl3aJLFiWtZRqMtKtLKSv9+MlpOPpa64LM3fTV2rrNdJMe3aEVS+Ae7deBNN3ymvd\nJKJSL/C6CJvOzaWPTPvQAVfVvzdJ1PpwfmRIKwrp94XIu/yDkhXLth7crMUwMWbbG4qVX29U\nFRmJoUp5bXGVrdYtv6U24SFTy9eHKh9ZdaNt9F7W0Vmqq4MVhXRKh3/XXwr/ufwr6HNcGIOb\neXl+y2dO6xx7/QRKIdUZJL2Ca5Jcs7oNaoRQarGQjMQ1b80lcuEKNjeRWlPTXlkHKwrpcPvB\n/TMjP7o/XmqdnnBap9Nmc1rde5bkajdT6UVf9W2LISwVYgSNsfHm5upgaCbLil7a1fhdUUim\n8Ce6sDZlcuU07fpC9Kq8EivHF6XMUVfVjYshLBNiDI2xyXXVyJKhmaEV2bzHt7ZurzsSuRAT\nB9NNR9e/q65TemuqyLSw6oduC6lFmda8G10b5T4iuzmLx9K0j7RvnssieU06dnFPtqHaZXnd\n10ifP5dn818jea8TrdVkNAlJ5E7Zo61GIxuJr9aEbmv3bt1DRgJoG+W6NMaWm/MHx2L/nDnJ\nF4aj3mrDmtXlSv9PwJx5+Z3loh+tOUC+VdfUDu9V02o0PvO6vL4OtDs7gPBXNXcro645tryu\nunrIh8w2DdldU0jd1+nyPtLh+B7wPhJ98E56K5HTumByzhX6FKzt0bR4z69q7oWaOrDb1Gd3\nVSEFuihuEOMrEWs8dCqa5I+z4lNN3YWSvkKzA2jbKX2Xj0zFlUlKndWmPrt7FVL5gplcQ3YN\nO1LBJDdVUYjC2Gis9qHxHWnKq5oJXdqMlw9UdQwNsDy3tyZmWqpTtVshjW/UQw5GhdSvOlXK\nZBqtjc/vYax5KVmf4zEaGo5qcC1FoEZaW4q8XbDWx62RpOugKiLZSkgjo5ptfbimq0jDeJsh\ny/1+dD3YHo91qpOnhk2PBxVWW0u+rFIjrU0Wb7eo1t0AHB1d19ix7hM8tncxjJgLa7X+66xX\nrybjUqOt2gvR72G89CHBN4bZFMxSQhK2a5Nl5nbVV3JGem8KumchLeoiulTkMtk6y6V4vB1p\nwSpa0L4aae1U8HaLar2eNHxk8kmFFF4qc5fLQg/nNVKb/TaWtL/jHcng2Xek8FLhe1HoayTr\n2jz4lXaNy0DT1wf3QFXHDV4jWezgNdLCLqJLJVkExWNemy9bQgvaVyOtTRZvt/RaUkn9ZD+o\nkMJL5WaPP2wYz/2iRlo7dK2/0LimUS3oRxUSAKuyopAaLowgJLAzVhTSB4QEHpY1L+2+D+Wf\nPAlwAcA2rPoa6bv8db4IFwBswro3Gz7It80XcgHAFuCuHQABQEgABAAhUbf38SbgwoSMcuVU\nUXfruW7xBCERr/fxsZSFCRnlyqmi7tZz3eQJQspO7+ODkgsTMsqVU0Xdree6zROElJ3ex0f3\nFyZklCunirpbz3Wbpx0KaaGNvfBlsgiPlTaWvmohoxxxlYYV2f0Qqt1fdrj+XfUJarsNdVd2\nHUmjp/0JabFrZHcFivBYaWP5FwD8lYbvqj+tmpEDVd+vSpSR0Jw22JGmUnCx3DWyZ9k83lju\nlVHXNJspNLHPeLbcZnR7sYued+hbVg1MuUrqDF4jBblYcEVyVkTLI12v6+xqG37LUrPZexbf\nIXwh9Y3pn8MzPwrZYWg//OGFr3uSrY/mcOlNOwek91a/8QT7U4KKcrHoNbKZKctj4oxbrVvd\nKiwGLMjDClBwJQctdDGqBvFIx++51T2TNlByHQ/3VJ6avQlpxWvkgsdh53CvcFjbypbkeqol\nmGmMKHJkR/IDGNuRfL+i5/oTXWYsYRMszglntoupS3JpMWndBOiKPhoLvZyp2JGyQ9eYMDqN\nUSGpraSqnzrPhpTUP4jn9qy79GAXfPaZwpEyov2IsHcnpInXyIVeo9dT+nS+MBqVR8OVfar4\nZ4bJdE5LxGBnZNT9WdVqzK08nyhFffCeNTvS0KMiytZ0ifZjwt6fkCYtxHy1U+ZS2ao6x/tU\n7Eh1Uafxf2IwD6QwpBpEYdsN6KM6MdaR/514vKYJKUC3oWyi2urOrekyTPqBd7sU0iSDfhYG\nNbSUJFvQR+vKcz2lKVmEu1q7M8MKpd7v+BYi16hS1bemyzBZvhautDuvyzou/LQX9mV34R11\nxcQ03rYizNrLJqreyUqojD6c6ryVp4QYSCoXt61dK6u0lIpjqn056AcSUmmg0TtSvjKqXFrr\nwhyxx7rO3JHErroiQ+LmeKb9zR1JL1+iiWmqK7UvC7st/oldlnTRD294sJuMnXSy5CWvYHA8\n3NT/pHRbd+5zcgQ8lpkmDHO1TV0l1W1WHUmByoxOjpsu50RrencvpGFCRpec0rzJ1YtcNRQU\npnaa+ninLMpihKX+VYZzIPUxFDyMjoecLc1ITV6YZJLeoAzz3oRZL6ekyQr2LiSyKPUH7CaF\npPgpH6w7vZSR0cATCbhxRzFeCZSuEEeDUalrQXsYHQ/pQuvfDavonZgcnuX+Sb1G6rwJ03XT\nH25LzP6F1D+4O1Lriwkyk/Vdq2XhX9dX+KjsURdMXizaZ9TwMBYd61LwXTXKQiNquGBELiN0\neRr3b9hraz6xy2IuyHriJcDfqso2U/ndQ7fTWLthAhvDmiLWkUburlARivQwOh7aJckyVs1G\ng/dSIYQ6PgJpqn2Bu7Zvaz6xy3IuyLC965nGzNCZbNgD+GOpJb+maLp8qL18rAomTb93ZnkY\nGQ/vIrOgG1Zsp3bkzEHNCHjTCQvctX1b84ldlnNB1xM3s217NddmZdeWy67+Arz9EqKy6qtr\nsStXW5OHsfGwLoWXIZV5cQP3NzsvHN70SXekmkW1cd2V2qzq2iK5m8mp20FYMDMCsDyMmONd\n3A1pfl7q+rs5mrLAdQ8gpJpVunFe6ExUd62f/nbbzVQGMyMAy8OIOd6lEOHcvFQvfZFC3r+Q\nlnA/ZSaXk8UElg9mSq3dU4bKF4cTpn9CAO1d7tAFAJFASAAEACEBEACEBEAAEBIAAUBIAATw\nDEK6s9uu4BF5AiEt9wkCAHoeX0gTP/KxJq2fEVgpigf3G8szCIk+3CNqy9xkD91q436QC4aH\nFhKdIm+2Np3Ei3O1ZSbnY9m3Q8tETKPQHtqqvfEDioPfBlvicKXHa7O2BNaarrc4o8s2LhKd\nfqcUNl0PRYD5sPuZ0Ex8LMOD9tDmtropaeg1r5u2So9TEljb+uGElEedF/Dhv9bn/sVhkjWz\nlAM3BL48Erf5DA+MPcTEkK1z5yon6fZLCHIEOo7qCAezfqkOc1h0wRNZ8jcSnrZRP5ixBiFd\n1nPBFrlEjuhXSuZ6yPvzxE5Yz0ZipU54DDJw+VovfEvKdrWH1P98mHBr56I2wiTwmxhCYi7s\nbdSwdrPphWfYqE73gwmJZClrh2RHrPD9HNFNSPRn9e0ukdNi7YO4WVaDoK7yaORAYuBOOx7P\nkCW5oxuhGKkec+jtJyom20VOZNEpN1a3ldYO5gGFlB/YxuQss1IXbGpEj2Gb6GoSWxcsWUvN\nra8zqyYqgN4Xddo/qdiRnFBoYscdDyO1K5vk3Pfs7ZiGuYZ/oMk9ZhsfaxDSJcyFNz9W0eVp\nGP7Ll31RrHzjUmsRW5T9WR8J1AzaWHJzbeVAEnto81WKQSz8NDEyHHpcKZwrMpVfTib5SlA3\npU1U1NkztTAupJR0As3xGJ5K7EtIcobVcbaAsIIQXelc03qevCNx+16g5bBNa9ImKYf83zp/\nZgR8kNSPcYidsZJttHZdywkqNDHiNryVC340LDOQ2szuSkjWtPPjvAVLAUvHLTtk1Uxy41LO\nSO2WQigFagzIXR71ucG3Vde1Ey7d0Afm2haXN0rx3EygGmFnGDWbmL1pEAUjdKwlGXkBV2Z1\nZ0KiD+ZxmqfUT0Sn0qGTqjculXJZu9Ro7lIK1ArEr3/3nBBUPlRxCSIM8UfjGBu0GFf2TxNj\nblil4fk5kAn3h5AkVWP1ItKH3T5Do9EWEV1iXHipYMfzqFVu6dpOp5HboOKwFiflIHcX9eYs\nstysP0ul+mLDuzV2PBaxdyTyYCZQZYi3YBEVQkq5SUlJI4tNbjgipdF4nMkouSetRs7HdAly\nMSxfxvFOpohPrSgzkm65I9UHaWjopkxj5R4a2ZGOmGetiI66wedUIemLJ2qKFLIxLt6APiOt\nxkJykmLHo5NFvbBJ1cvt9BRVdNuZkOzVxi3PLKSONxiqj1Rna5qNKewGIcm5J1UmvZhDUjHT\nVmRByRYnVQnLkhFQYq78BakbqpdREZKTFC8RKlrihfjVQvJqZ5Ta9aDR7LQuUS5ovVrHzabD\n/Bd3pModnDrIj2yWpEYTLYH+NAuUx6bNy1Z5ZLRyLSs1w7Clnb1mV2zDyVElpqAhvqqMOlfb\n1H0AACAASURBVEnxEmEMc/AyCMgS0hBQXVYc96V2rYY3FZI57eS4OMTmiXfLJ4aDrWlm24O8\nRCHGVFHKQpCxWcflg6R5HRhz35+VO1JyYpdLSR9mhXsjKV6EVrQ57zlCQ0ilYY5FWNFzX0Kq\nXR6ubWgNyIUsT/mE3BIH1IS96skVVdQMHxLTH9uDZGsqoU50bRoEf9Rj9IuTRtPRYOoDcZIi\nY+ibetEay4oUmzvMsbzp2rGblU8HdQlz0bKw6GVbrmPF+at2QMMyHDGZGRsHX3DJWfEHffCd\nTRiENCxOewlknXmLhpCcpFgx+NEaMy3NucMcz2FdlncnpPqZug5/KHVrFau25XvgYVk2yUkr\nENJF9KZNDcMRMhrPaEFGfNxThWTPjoqhEG3K1x4FyduR1URcleWdCampesbaxlQiM2UufPmk\n6S8fLu0OgcEGGpad24019rCa013NNeUmnz7MYG9CanopMNY2sDJTKlzD55Xa1QjtZQe1jIzm\nGpad241NEJ48cv3PiCTNc+WUt7A7Id038xa4sOXxuaBXh+0ypg9zglilyx26WIaaS+6lej8t\nqX+N1E0RUlDKIaRY5r2UWeyF0GOTCBP7zo9hlS536GIp5s0KZDSJ4V7O1L4BIazS5Q5dABAJ\nhARAABASAAFASPf0suSeYtmcFZMR4OrphXRPN8ruKZbNWTEZEa4eWEhVuUn03bzFqLNP3w/J\nnZpiGxrP/rxAe2fnowPT31OjnxUZ8S0dtSbt+jCx/7VHW/OJXbZwUbfM0M9pBTidE0o3fPZS\nfHisPrahceOIZiWgFOh0w3Ujdwbc5lZ/umFC2A8rJGOZMVt11YvfwqH0Our6/6WGvsJRU6+Z\nb++XAp1umFsdb8UctbnVn7ebEvYDC4k+lJolW0hxe1Tt57kS/eVd3qlWSLcHe42dHaDVk/nK\nXlhKFxQScaYuMc2+RiJU8ylhry+kj5eUjp+TXLRv11XX2KaQZl3tTAolFwK78h/rS4t3OCI9\nFkcjXpQ1kM2K6yt+XcbMVnlhZguthtbEJwmiFLEKiFwYGmGPRtzQdnqXa79Lx9frSE4TXDTV\ndv020JuV013Vfdy4ueLZdqmma3ckmhR3RxoZDblAqlIsOSLkP2ynbD+h/Uwv9pHx2aYDTtl1\n3tqNMUh/Sbu68x3pEtcpnX677ueUPppdtNV2beusIiEk+tCKWJfFTPkl0leEqkU3DnaelIm5\nxvpCYoLwxyTilhdUwxIgFM1q1BiPkRBfSEIErP+goP4v3Vn6o2Mnrjvbuc/aQjqk3/Pz3/TS\n7KKxtmtWNFILvOmk3V265tOcss2CNnJbU4N2/GxpHiqDKndkNCxat4ljVgbKRiqjtyaxVlzG\ncTXg5Ewo9Z76RmZAVWnXrC2kugq1V+uKnnTwdYmgq6r2NklIWUBDHImskGXTeZaHZ4UgZFIS\ndcqNjuhxLL1qpU+JD5DuP3mIwimrZDEI4dkOl4su8SC4B8sA30M7ncCOLQYN07+2kN76wA/N\nLipqu3kdcfcH53ClSfpAlzh+uKZzoy+vUcUaO6bvbpCIXPT1AP3s6X2tbfdnUTJDfPlyxnsV\nPNWQO1nVId06Vrec0+XaLx3fPz7Tv7+nv6fy3QZHSGODm1T97trVrMmhJ4uGrJUVO1JTVXV1\nQ6bFUdyUSpaIiG5ioJY78YebvdvxRF+SNqwdLD1qZRJYRtUp3o4MpTak3nB1yzldrv1ong+/\n7S5Ga7tlNWdLmfWtsEYdkdZaQEJIepJoNM4YnHB48Zrt+vGVhWSml2VJiqkbcteHLJJqmTIk\n3VK1ND0yVSwBxk5DRslMJRX2fQup+/7++Dgez90Pp6KOxt5HcmuKPxY98KoximhUtq49EuVA\nPmxZJs2cKUy6Ezmj3AnjfZWXc2MXf1IFR7YTarKcKzpAOcKGVJPORb9qXbp6yCFkf+xp58U4\nGlh1yzldolxYeWCn6UPRgboqkB3bUslbi6JR1W1ols91EpGkfHI0AnO7GzbYlk3WHtUgpuGU\nnXY1iL6N0VxLvzTOJOOy25GHnByqJ7tXonhRGD3rm87oEuSCjc8YZX3xq42f/ekcarFnKqHU\nOxe6UVVD9ZZtOGFnfTYLiT+QUdBUm2lXawcdqW5e6GqcHZ75hthJI2jXeDeiNa9nS+PJXWJc\njC9EteuIIyB2taEO1dtrDKtf1s1CH7aTrigkXip2UTeusWJUogqFqqRp1ZoNdSQd9cthwZCx\n5lYsju1Z6juu0sUwUrTi7rudVytDk6bKX2pHagwrqRfx3HJKuQLrIhDtpshIGpOW+X4jTevW\nzFS5lq0hOE3Luwv3XqycKpOlfqt0MYwoK3TZdDZr9jhtvNkZN6kWwfp1sb211dvZkXLlpfIk\nJ1Yss5KjbXYq+6M9ReuWbabFUR3j1zLzXazSJcYFv3yZlw4hWPvipGF1mnxJIHqbOxJZZEZt\nzAvEjWv2Fl0f1xLlzvICIZHleXatCAuGwTYf8yLyFZCokmoiiJMRNzZ3i66Oa5Fyz1t6ZHqo\ng1W6BLno87BUMrbFGVX0LjOZ1bboJct9sTSuKqSv9+MlRcfT1zQXd1BO63M3g15ti76bETew\nopB+X8g1/OsiLgDYiBWFdEqHf9+XZz+fhwkfWgXgfllRSIf0PTz/nvA1ijgWunSwLu3nunL7\nr3j549wDqb51MKXGqnq12J4SR5P9VuPTq1zcJVvCRW0gi7yYtW4KzHXl9l/x/oPlqtr9tDjr\nerXYnhJHU58H2pEaZna08YSN4GqYpz4V7+SOB2z3z/O7hpIsV9Xuiw3N8Sd7dEbblhxMyVdb\nn3VfI33+XJ4t8hqpaRUbaTxlI8g6ou9C0of2gM3+TqkthRVD9XumpYZ6HsQB8SaUzNXQqiIL\nU97kbXwHusHy9C5XXknqXiZ8se92zk5cXWnl/ahL/pcSXFsFJ1RH/Xn52Bqw2Z/GHiWkQjFa\nMRTHNdaZ+Oxdi6d6mbB3xU70nxZHUJ81hdR9nS7vIx2O7xPfR+rMxN2yf3vu9upPJ/qDpgUf\nppCUE7J+Ei357R1bySwIepIdGxaExotQ66DcF+wOKdEsDtkcKWLeXri8PSMxdP3/kmxFrDHb\ntxaVShrb/71td1yq+/pkQ2euTeYmoHpladAi9HzYpWU4yQ0Hs0xIhR2MPibbpygzs72jAUeY\n+mCllbzf5oOFjnQAqhE5kMR4ekf0oLM5DLbVGSeOkWYyzNxnbJSXZuUAgrrEudBLC826u/jn\nlHQdLQTPiZN2ueOwvSsXgzzviDs/DOVpFAupOXbIOJX72EPwWid5k4R1IPqh5m330gzpnJhV\npgOyCnXEJZ82NSFkHstatsRsj9VcJN1yYN1LJ8O6hLkw94Tbg5UrVpxDTsz6473s08qJmPHU\n2xcB+E5ygEMNsX7EHd+dcikL9bEQZVXrg6QW9QZMpEMjY16sfIhgxHG6F/Fhcpty2nLg1OiI\nSGi+zFM0C51OpZM27ad0MqxLnAsxKJ4JlSx+kK1147m3zidp0Zhh1c/xxAOh/+d2WP2x2KyR\nDPXI+rIotOxMIVE/XV/PiaVe+s9hcRO3fmz3oZOTJHLcfJSqOY1LzoxMkHXK8Erb67Rp9iek\nPFyzpszGpKpI+5ITbSyfYcKhD7Zdz1Yu3t7mYJpPqgyf9mMDvLUyy8KrlaEMByPkbxqTWF9y\nAbLljHkiUzOMj+8ReWc10m5MW9dRszwqPSr1a/7avVoPzEmnDw57ExKfJpEEo0b4Y617Y2o6\nasHdkWxjbpNhiujEdnLqiAVVzLSmbo1obdNFXB9kY00kGLO4eDa5Qbm89dGyDuK88Dc81ZkT\nM6HCUbMlJp7KUNSOHG7HYhufwdyodDKsS6QLUUWytIY2tA5kpvMfcpJEb32qk9NiNxSrtD0N\nNAhqSQbNR2LUOomKdyHp6KghPVLSVvsxGlu7VK8etgGx8Lte6tQAOcTjMw9yedJHlViaGis7\nPH79RMShJ5DOZfFsVJclXPCcDH91nUxPP0+8q0wxO6NnL59lfRI9pJrmVXl8R2IOZRCkWmSA\ncvZZSzlcUX7WUOnWIGISVW+ETsZMxk68WzpivlgsvKXIGrWQ/2S5pqapZ9urrig9mz67FFI/\nLLJ4ejmiKyPpfn0gz6jhIW0id+KcvTEQFznIrlNB9AZZPPzEcNoMgO0bfEz9wFUsSrkqMrFp\n0FGptYP9RRVIdlURVEfO0waJN+hkRDI7ndlX99GWzXlLdHWgaaUO7bVwcFk4F9cl1gUdOxeS\ns8Kp/v2DmCQy1XSeO3KWrekskI5mmq/3/JwahyE1esYSkvDByk12o2GLxVz5pa3z4Oj47TJM\nxDQxTgehD1i6VjOWvVuJ4yGIqTJEQ2RCs0mHStcEHsOjCSlPKE17nyU5Z3r0pDCSPKKWT9Uv\npcR684W+o4b48mePI+kocuzMuF4qeV1kG9whL1ZzJGZUliJ5vnOeO6OZyPFgnNWvMJyrWqa9\nk3ZYlCrB+RALmI8iiVaJ5UVPzOMJ6fZAks9WOJJgZyfoH/gkkd7G7Ih00jLUmaarm5plPQ7R\n29pF2IiHv0nbjvagHlVNq5GIsHhi2ahUziwzznjpsBJLHC1xGX2notBGjQyTI5ai6Lxo3cjZ\ntxMlHBbOxXWJdMHHy1YkKjEnwbfzN0NiknL5WiqUs876sNCoBgpRkHHw3mxfsNbtjh1hHZLy\nqc3bI6GxU1dyreePFZUu+pJh8a7Zp8i7k0HaoOiU5IUP0I7fyvPoCPcnJFW1rJzoETfB+ZyV\nsOR096ahszKt5q40DtuZCkJHxWMnY/ISZI7Eyk+/TsmqtpbnUqqtIbMJEyV9fa5zUWG84DTb\nJU/c+O08jwSxSyEZK4g4NTJu1oNXGFGANMBPsTC0x1zVPFhzHI4zL7rymPKjHq45EsMU1bfe\nsOSIqiudzA/rSudNGp8NtWlHqqaukwMfH+EOhVTSScXSMWq69izzVdDd9XlEsOUOORT+6Fkq\nnMzytveryTn2+g7H505gi1O3OX+s7NUW08QuwS6K1RQ8C9PCGM6XhNQe7JjO+0YFn3Vuuv6a\nrmIBbzQ9koolJrDR5qT87VJIO2KJa5VRX3N99ld9jpAenUn5ewAhFa+v3KPVRmcyel0RuAT3\nvmZfH6VMUGh1Ppdp3GxpytD3L6Q8anv8kwoisIpG7ITWa96MZlpcXUZteYiLzq+ZVkMTfLd3\nWdBF3oftHbl6n6bJI/W4cDkVwov2XGVvXfUwv9eHQgN65VFuPMftxAzsV0jDNtQ/t18j1r1y\nFFcy5IXvQlJKfApjbo4V21fZ441WVBXPg3bMZmL23RTHLXHUbmmC8/YuC7igmbVE0PG/au6v\nJbEH3Z6SU4Xu+mmNQx4eGQ7bFH3D7KhVA2IwI9ER/449d5DymNG5kCaaBytSNgd1c1pzoSuU\n2XUy7fWa2quQhgpP4hcBpuxIwzTJhY+Jyy1D2aXqMvIWFA0vJe2KjNEwzI6qomANKoXEXSp7\n2qtzzMhYMU05D8M543preKB9/MGUUsfCoEIS4ddM56VHTaPZXeJd0KLPIrid4Q2dcuBtUsdt\n9AKyLyo8DzXOSAMeNPXN/zBWTmNsRnnlBkRP5thzVzJsZc/wah8bapgdUiM1DND5pDGKLduK\nTKptJHXJDlJNxfh0XhpVtJnfJdyF2inI36rp2KrSr0l0g6B5HhUSs+K0MnrQUXSysxSSZVg0\n7VQbFhqzlFvQ7JAsJMue9qqO0foUNUocWhZyujvZkywDzo6WD1L/bupk0Wgb3kg9diokuV0M\nh+wdg61NVotukEy2Luls+6Tg3Npz+oji6F3nR6lfWQ1WXZI2QrDU/NBC2E95VXHKqCBYnjcz\nJ0kHnEPJCxrvqV+nGjOh5soPlyaXWeJuvAyYPJSQzOSwHnZOhqlSNclFVFRSXwJdReZF/fa+\nLN+0vTScjEaqMPpu5gDkNpH4M3MwxR2JZsvYkYb5cdM0LGpcAn76RSqYfy91bFEhEeaJJtnj\n64rHXoVkbRp6RewbdtbkiTZ9Gmkrw4lOaO7j1p4XPd8uRaC8KvjmxR5JzfAQ2d+Dw65TgYoM\nqr525q0cDNeFt9BZmdYJSdwe4Fq3AuqoQ5EyGgppLYJOJOhk/F1aR3ub3gmf+xASm8/80Mlp\nEonMmRCzJYux/8tSkx2J3MCsgOUfpL2eSRaPnlUVrdxFhDHVrdMlpV0WMm8e43HJIatQKqzT\nmePnhqc5Kx3/D0+Xk32aWboOyHOlqB9GSLdDtxPGH3x981fFJO2K6bByeT2U1BxYhuWx3J/N\nr7Dd5c75kayluTVLRA5euOkb8SyQxaS0ayivNAfd8KCGmw8UClIMvxyunqc8SWRacgo61nH4\nb0dtDOHL7FvDzja94fjciZD4Ckf/q4/RtWvozf9mdm9nchLFBQCLIT/vD5DpNayyg/ScKiE+\nKNqYT6ocVDe4J8VteB0My1XGyJfZXwWcq9aut4I8fbhgOpKTIQzSTCpYD0Q0otYMV2J1fCQh\nZamI6VYLLqksMYW6drR1seqTVPPVTji3U24VJotOl4CyOhzJHkRdqlhYmqyh00iJIEU96ZBI\n1DJ7xjDkefukhz15LB/ct7HSCWmp2hHLMF0XyfrwUELiVSX2HXmMr54kxd6KOzS7pZsqJ1cv\nmzE+o322eb3Yws1xJ+ZXNBGqyQGpuhRjt5o4azQPkrxc6BQkamtIZaW0y8iJnGaUpYanvT+g\nZl4+ytHQYiIV8DhCUtohqUsU2VqlmJihhwdBdHRqkoDFkC0yvY3tSLIAbCH1RU0bd/1wZF2y\nfLDgi0OXaaSJUORAjCEtCoucq8qUbz8PfCAibk+FfKYHFT2QkOiDrhFRN06KvRPDjtPx9DI1\ndXYlcQ0pIVnCpT2cPSDHyRs7U2rWgrrMVEPnKSsmTtbcikJiQQlVuXunGoicCjZtwplaQL3R\n7k5IzdPnjt06YWaeT15SMdgborRuTwMx6AipY85p1dSkIJlt/XogAyo0YmlYU0h1a2RuoPqY\nHYtWhgSWve1OSAtPH7lqqonB2GXyOqY6FaaheDFlLKHa71ioFW09f/OaLUf9PLV2bHe2RyEt\nOX111snO4CimcaYSvXA0Txt7ZJWHSemqtD9hoOuyotT3J6SJ09dQeBVN+1Z204nhFRwbVzQt\nl7bN8VT2CJPRUnpcT+o7FNLkjTmiOG5n66enRcDVBpsSEFhJW9f7YivCfPYopI7Vs3kuyb/7\n9q6pysAYpdDoiyUVzojfUkhD78poG81bp8nuq4ddt4uWvY1nRLWt9+I1CJbY7oTE5zQrhL4S\nZ+mmDT3Z1aY0z/j1wYkvkVc8smUOxYi9IqSGmwfKDgtspFduIpNtrVJj3l2HbHJGZsJKe00n\ns0HTvFewNyHx8uwTyw+mfO76hJS2Ye9moxxQb7Xr1INtj5UjkzWx11n1IwtZ2OeP5Zh5S2a3\nsI+Qtv0z8tCJWB1j7Iw7HiEjaUomhj4U/JNuXoCF49Pqe2dCYiohl1D5by6y60FfSOb0mW5J\nKSRjRnNL3mEoP3aemrLLvXBd5voeaSmzVupF2hI95eDGo6FnVBZoq8S/nmwsKdxiN+qfz5Ud\noHPcX8HG2J2QOrXpsKRZxZ7I932kuWE9LK1Fw0Sbs2NPNvlmC2tJ3WqTZCpdIbkbykjd8RVm\nTKd0w5QhE3/Ch7RjDZWFO0xaGj6U5Q5V5Ds5/mk3L0DneCG5Y+xLSGqayV/FHcleapKGn8xu\nWQPmn/e0J5u1NIuLF15J2OZI+GjI37pmeCyGFTJ2YjAbT3SMbuVZjcx9Y1hvlCljjeHjHEYg\nhMTGqRvYxz1zlexLSGxVpbOdT+Z573v4X27khZ3t53Nk5gzdddyf6qJbkmKkJlmZkyZODpyB\n9HXJa9ha0HPUlhmWYXvwfRNhh3nTjZL2RdOhZoA/qsXM6mQsGLJBZx6nufKWKX2UnC+dDOsS\n5kLXZMprpLVD2DVNjPHjevXLFrJl3sNRjnKpYjbiEvVbuTRytyIPNBHysLlHsBEYQyMDEDub\n/Ud/LKncW4tLEvF03BHJlJcmPpVeHnVyvLQUzeQG/qnALnEuzKk2ZpBtVLIGiC3Vt2NPSCHQ\nPDORMSe8mQqLRa51LJRkr42lnDAhdcyyLXZuJqePRMVHJm2Q+Gle8jwwmbNYZI5EE8uh0Ed2\nYo2Dh2GtZLLTEJWZ4uJ07E1IRDv0TzOdKlX9vAtbia3iHbHeDaXZ92bdbn8Mx6gyhrnWgQ5n\nO13uoqZ4Z5UENtRhrrlho6GxXVFfwkWuIFKCfDhipMwmi84IWU+d3SY752fdjsOkyDGLjsZo\n/FjMSC8t3DORXeJd8HQOz8l5uUnzrNLetywza7kocor5fiFVpX/iQUeqosixdyIm2llOojOK\nPqr+iTH5VO3ZOwuBx0tywIZOhiMdOasXjUSYseuTrjHKj8if7KnCShI5Jmd1GdqQB5O9Cqnj\n6ex0ddGHvqlMlc40sdb3F9b5A52ivjXdGIQWjLoQcXb9pLJRsjbG37rYaHWwtLlLrxUvy0FW\nF08m26O4kPi4Le3r7ZaFk5skskqo/ImepOwTGw+LWTbOU1e0aLNjIeU5IKurmEy1jOdFbjgo\nn2phaeHwLYdONhUSj0jOfB+OiJMMbxASH4v+mwdgSVc6NaqQu1SrB8khL2+5liurqkSpNlwh\niaWS2EkMpyPvz5euYZZoYzIYbZBlwWS3QnLXzI4+SiHlWrAMiq5MKXYRGAH1OVf9ir0Kw5OL\nglokRFHxhUAaz8dVyYiMWjmQFUycpGSPhZQnWwrI+Lz5UJKmg/QSZ2w1gxNqhateTJ0OxfM2\nNPFPBXZZxAUrHV1lanWUOTQMdrJP7sla2UnlFab7+b284dlxqTiT8pz7K/OFkqEJFS741sBr\nMrGDXlr4CAZbfgIsF+6bZLwfeRiykvheZu9IYtAqEjPOaxP/VGCXpVyQoetqUyN3C0hbq/Pp\nnbdTXp6IBo8qAn6AlobplZa+65I+mKNQeirFb0Scxruq4ic2Cj2zdlhD4dNYcX2bozrauZDY\namfMVaF92VqdT//8aOBNSI8qArs0SmU6VhZVawp10zYCppEKD4aN8XmUojFkJMz6Nkfj3LmQ\n2AVQ7YzKrHAbLT7dBb+mawOymzJjlkZVmTom62q8HttYm4v61reWo2mbEITrcpUui7nI822u\nNoUu6jKjxaXhv3So6pxpvy2i65NJGlCd5pWX07tgNFS2dYwsiI3GVumylIs0XIHwVbj0orBv\nYF0k17nMDoy+JXM1rtplQEZEOrfVxaROFSFVHR45tRTU5Xz3exfS7YEU0lVWbI2R+eqbcRus\nJfXh64Q/eIfsc1Z85KKkavwsIhlaUs/UuIjbyktjZUI/dWVZUCs9tZKgpMsk57ExFxMiaO8S\n6oJuCPRIIl+1ND91yV9tioonpVS6UiO9dBzcXDlsKz768qacIGaXmbcGahUFO57UzbFRp06K\neUhOpIVB5MDYqLTzUmBV0Gj00I1clE3vT0iGCshskoLoG4h6yq06YaPvVrhS07KjR5k5P+xb\nUKWn9UIii4kRGhlTEt1EnnJSkrobYzh1R9PpZKlIC4Mga2FiSTPDVoF65/SAaDQ5CdqPadpK\niXfCZ1shDQm+/ZH4g0q0eZYLSVYcnYkkrnnoKSFdas5YjWkkQr7cnhqeNGTW5xADD/r2hyMk\n2ljvY+U5oKq3FSQVObixLOYeXQ5ZRsJmif7Nk6POWQPiFZKHRTvwI3SIegTO8QJbCynnhZYK\n1UrWCZlLWSZG0pgN3UFUjzZKzRlhs0JgbcVAeCjCjjzGwy58bcgyYzZOIzWTC1VEpBYTHq38\nWw2iy8kxhTQkxwk0j4oXvjmgfLA3JCNlq5ZYK/QIvGT5bCqkPke9QkgSrHKQRVvetMXewF9v\n8YWOuaHn+Z5Bp5ae4zGSwEhdWtPPj8kRleEhEV/MSK4n1YvOwRDL0JuGyp3yVcMuRR2ISgRT\nWJcXh8Ek+5sWvi0C5kBERk7w0T+MkLphwRETLnOiCoPOFKuLjs5Qtk+KeujZ2/GqVwerakIe\nV/vPYIiIjAXaJyGxHrkrO6bCc47zeFiQ5DmLYuio6s4qZhI0Sw9/7oVsZa2TodEHQ1NKBCqB\nnpBURAY7FFJfkXTShhE6E2GTbXZimm6a4vcuRKKlMZ1kq02nXBEJdeyorN9ODXnIRC257tVw\nxRG2+5Gn7Iwcghet/kv0L6SV7Uhy2MNSkxI7wtPHBe3PkHmiEzlTNm4dnOMFpgvp6/14ieV4\n+prugo+vk7PslESfUXrcURLpkp8oIWnJ6iTnQHU8ZHUmtUEXadKb+iUGOmldDZN5Zs+Ymunf\nOqFya7GWaDoaNv5O9qZGdHtVwjJtvADYV53YQjSMbEivISQ6AdZ2RUfPY1GsKKTfF5L91+ku\nEhn+9bEjD1Q17BiZOLVC5dyz3FJzibc39j5vLrrhNC8IbyBinGzucg9SPjIQUkUkF0b598Pq\nRyeGPIyVPBVzQOVgZTU718JR7R01mGukmJo+EyLLSVmUMyRD4OlX3p223apCOqXDv+/Ls5/P\nQzrNcJFniE84Pc0qQHaTPaglkkAJNTVUnnmeGGWizCc6Wgil/nwUvHaydRGEGTrpKTInRjys\n9WSoImEVWVV5kX/JWRDqll1yhDzi3Fe7o/2s0GUIcmDctt92VSEd0vfw/Dsd5rigpeGthFpI\nMhmyQOVTUXbCUlFoIhAxCbxuhXsdLzvVdxVWiGplfHIgsuB4Olka+6CLReRlla0hMslGe1ax\nNCKjQz7MnI850KGPtcgu7kdIoh7CXOgR+qXt9XDbeUaYfLx2UmR6sW4dGTvT22ByYZowNgf6\naDgSSkqmOT9a0qq83Buj6/vSrKnVSxv0QxPtjVaVhWDOoWxTNhHU5ULgjiRaq1SOzXtVGXfe\nRkFN2FMuG/HNYDhR5T93MN37f8gNl0SgLiZVnHL4VrjFmCqWe+2Mhm54tDPtZ1KEeUkukQAA\nCgZJREFUY3QuzJ1yUW67opD+XiN9/lyezX2NxNtWFqXsUOphnnbrZkSx9ElrsH1j3Y/96VWk\nbE+juB6kPeq3at5LxpStjhij3XNkIy5rQjPaO9G0TEd5plvimtzlyiuZ4pffSS74AtM/FgY4\naRGymqjarJCRaSCRS6i63rTuhWDkM3l9p0+mm6KI5XyQK4H1sqLqSFvSiplzxGaYKzhrTLVy\nJmIjp8icDNu12Kjr5mpNIXVfp8v7SIfj+8T3kfI4yeRfT4x1EEeyFbOTsskqMNeJmIBCapKP\n30nEy3roXHSy4kRp5JPickdcEpojV3HyKzfu1hojG4IxTJ0mM3dWfkZyZ+RAjMqK1x2I48s7\n4TNdSHNd0Kmmy4q7Qw/NuI2umJvk9ZLFIERVTHTKtVS7CPAh83nvaAaMg31n6jf36FsM4RBL\nchT0KF/KBx+GJV6wbqWygeY82fLndnh+1EyZPZl54o1PTMemlnR0Ir85cyawwLZC6sRGnc+U\nZoebyJMszjJTYnpob3GIVYwXekdnJpFRtO1jfQ+xaMqDvVNWffnkkEhWVKpMVcJ0eugYcniD\nHdJrcGJuI8k2xMLSe5NKoRgpU4JYTEg/Gd0wXXzSrLdUBsf+JHpsJiQ6G+rI9UF24I/989S/\nlarOqmQOnXJSzVUsV4wTOjXNl4NqIZEwOpEBPsc8UUNPmitVjKxQeQC8jIhWyMhYuFR9g5BI\nci0hSQlkL8QotUkdD9ZEUNQZt8+adzI6EnMekC4LFr4/iR4hQipVj+uCrzQdmxjTqCEVUs2O\nkHLllXoNU9z1GS4JiSgpUVsjqcg6oD2po6xfdpB4IUsHz1X2z1VnLCz5sw4kZ0Zl0qoj+TGX\nfmOwovRJRdNc6/FSIZEH4YyFyxZDoXMSBstLjsiaKmcKCywkJJ59p5OqAXOh0R0sP5bMWMGb\nwelqNafWDF3UlxWd6iZcEyH1c51yifO09FFqO1Y+uQtj5PSWIat3meg8PGp1sNZpJyJPQgJG\n0nKt53DUo3DGEqRGnh1RISX6X6IxI3B7PCVChDTNhcgoHf31tNPBPsZmP5/trIkmhvi8kHIu\nCUmsE/yvQiJoWxVGHoc46IzWO+k36NhhZ/FRjpMJC3o0TzcvQkl9V75pZGs8vNxDJYiN2gjZ\nG4c7w/sSktAPPXx9MDs4RnjSxkyRhrqbZ0wZoJNoLP92L15bLIxOmtEWc8/CSeKrFEIh08qS\nSEr/rD5P2R8xNDyQXYdYE+EZs2wkkhzX5oa0SUMybHdA/kjbu8S5cEYyNjuOGatPhSm7IOsD\naAk1CSaG25Icz0aDMSNaZqHZ32BNZsLwQo/MHrgZlHVugrnJgYR8sc/pEZew0NyHkJdE92RF\n/9BgalrF+LwYItZGlLzV7K0opKgv9gFwf6wopLgv9gFwb6wopMW+RgHA5qwoJPGicAkXAGwE\ndiQAAlj3NdISX+wD4A5Y8/Z3wBf7ALhP1n0fae4X+wC4U/b2yQYA7hIICYAAICQAAoCQAAgA\nQgIgAAgJgAAgJAACgJAACABCAiAACAmAAO5USADsjAlVHi+cXfi+kwg2DwARBAUAIT13AIgA\nQnqECDYPABFASI8QweYBIAII6REi2DwARAAhPUIEmweACCCkR4hg8wAQAYT0CBFsHgAigJAe\nIYLNA0AEENIjRLB5AIgAQnqECDYPABFASI8QweYBIIIHEBIADwOEBEAAEBIAAUBIAAQAIQEQ\nAIQEQAAQEgABQEgABAAhARAAhARAABASAAFASAAEACEBEACEBEAAEBIAAUBIAASwmZBOh3Q4\n/a7u9uNlcEsiWDeYr1vSNwrg+y2lt58NI/i13a4VwUdf86FhbCWk18uP/r+s7fZ0cXv45RGs\nG8zv4Zr0jQL43DoFP4drBD/bRPDd/1sTtu+pYWwkpK90+O6+D+lrXbff6e33vCS9sQhWDuZ4\nncmtAjj8+fo9ptNmEbydff8tadtMwp+La83bvieHsZGQTunz77//0vu6bo/X4Z5TSSJYN5h/\nt39+Z6MA/l3K+DcdNosgbTkJH+n15t/2PTmMjYR0TOeN/TsdN/F+TiWJYNVgfvqZ3CiAt/Td\nP90ogtuV7VnK60fwt4zchGT7nhzGRkIiq9L6/KZXFsGqwbymn6ujjQJ4Sd374XKFu1UE77dL\nu/ctIviWnoTvyWE8o5A+ztv3ZlX0r9tUSCkdLy/1t4ug+zjfbTh8bBUBhBTEz+HYbVVFl2uG\njYV0vtnwts1+cOX9cmPsvYOQ5rOhkH4PryKCFYN5Od923lhI59dIP+cbvBtF8HG+tPuT8geE\nNJ/DdkJ6fZERrBfM2+We0NXRJgGwStkogpd0foH2e5byJhHcXNi+J4exkZCuN0d+1r9r9/Py\n+iMjWC8Y+u/PbxIAewdgowjSxhGwu3bS9+QwNhLS+2Vp/rzcv1mTz/SqI1gvGCqkTQLoff2c\n87BRBNdF//JO1kaTcHmwfU8OYyMhbfTJhp9BRxt+suE2kxsF8Pfq6Pf8CuXfZhGc0vmjbKfN\nPltxE9JjfLLh70L5zOt4w1De8oZAI1g5mNtMbhTAu+l2zQhet42gf/1j+54axlZCun4CeG2v\n5MqKRrByMLeZ3CqAz1fD7aoRmG5Xi6AXku17ahhbCQmAhwJCAiAACAmAACAkAAKAkAAIAEIC\nIAAICYAAICQAAoCQAAgAQgIgAAgJgAAgJAACgJAACABCAiAACAmAACAkAAKAkAAIAEICIAAI\nCYAAICQAAoCQAAgAQgIgAAgJgAAgJAACgJAACABCAiAACAmAACAkAAKAkAAIAEICIAAICYAA\nICQAAoCQAAgAQtoD9r9W3/5v2IPFwFzsAQjp7sFc7AEI6e7BXOwBCOnuwVzsgYtkUvo5psP7\n5cDpkE43IX28pMPH3+Nr+vr771d62y7MZwZC2gM3IR3SH2clvZ6fHC9Hj+en6bXrftLh78/D\n4XfbUJ8VCGkP3IT0+tt9pJeu+5cO39334Xz083zw9zV9/m1Nfxp7T/+2jvVJgZD2wE1IX7en\nx8uzz+vT8w70m47deZ/6uDyCDYCQ9sBNSP3T212G69Mb3fni7u9l1IZRPjUQ0h6oE1J3Sqft\nYnxyIKQ9UBJSboUdaUMgpD0ghHQ831vovvLTK8e/10ivG0X49EBIe0AI6TPftbvcwOsuNxn+\n/V3YvaePjUN9ViCkPSCEdH3z6O3y9PKWUjr8dL+Hy/tIuLjbBghpD0ghde/skw3p7U89b7dP\nNuDibhMgJAACgJAACABCAiAACAmAACAkAAKAkAAIAEICIAAICYAAICQAAoCQAAgAQgIgAAgJ\ngAAgJAACgJAACABCAiAACAmAACAkAAKAkAAIAEICIAAICYAAICQAAoCQAAgAQgIgAAgJgAAg\nJAACgJAACABCAiAACAmAACAkAAKAkAAIAEICIID/JKATprWM9RAAAAAASUVORK5CYII=", "text/plain": [ "plot without title" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot(-log(p.value, 10), type = \"p\",lwd = 2, ylim = c(0, 3))" ] }, { "cell_type": "code", "execution_count": 85, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 0.0301309681396463
    2. \n", "\t
  2. 0.0175916492958266
    3. \n", "\t
  3. 0.0301309681396463
    4. \n", "\t
  4. 0.0374504988094523
    5. \n", "\t
  5. 0.0175916492958266
    6. \n", "\t
  6. 0.0323186543751581
    7. \n", "\t
  7. 0.00988457403114586
    8. \n", "\t
  8. 0.0175916492958266
    9. \n", "\t
  9. 0.0351733047983836
    10. \n", "\t
  10. 0.0207723954791442
    11. \n", "\t
  11. 0.0301309681396463
    12. \n", "\t
  12. 0.0323186543751583
    13. \n", "\t
  13. 0.0376688357105577
    14. \n", "\t
  14. 0.0374504988094523
    15. \n", "\t
  15. 0.0301309681396463
    16. \n", "\t
  16. 0.0118078833271804
    17. \n", "\t
  17. 0.0118078833271804
    18. \n", "\t
  18. 0.032127101177511
    19. \n", "\t
  19. 0.0118078833271804
    20. \n", "\t
  20. 0.0374504988094523
    21. \n", "\t
  21. 0.0323186543751583
    22. \n", "\t
  22. 0.0175916492958266
    23. \n", "\t
  23. 0.0118078833271804
    24. \n", "\t
  24. 0.0376688357105577
    25. \n", "\t
  25. 0.032127101177511
    26. \n", "\t
  26. 0.0301309681396461
    27. \n", "\t
  27. 0.0351733047983836
    28. \n", "\t
  28. 0.0323186543751583
    29. \n", "\t
  29. 0.0374504988094523
    30. \n", "\t
  30. 0.0351733047983836
    31. \n", "\t
  31. 0.0207723954791442
    32. \n", "\t
  32. 0.0207723954791442
    33. \n", "\t
  33. 0.0207723954791442
    34. \n", "\t
  34. 0.0207723954791442
    35. \n", "\t
  35. 0.0301309681396463
    36. \n", "\t
  36. 0.0207723954791442
    37. \n", "\t
  37. 0.0351733047983836
    38. \n", "\t
  38. 0.032127101177511
    39. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.0301309681396463\n", "\\item 0.0175916492958266\n", "\\item 0.0301309681396463\n", "\\item 0.0374504988094523\n", "\\item 0.0175916492958266\n", "\\item 0.0323186543751581\n", "\\item 0.00988457403114586\n", "\\item 0.0175916492958266\n", "\\item 0.0351733047983836\n", "\\item 0.0207723954791442\n", "\\item 0.0301309681396463\n", "\\item 0.0323186543751583\n", "\\item 0.0376688357105577\n", "\\item 0.0374504988094523\n", "\\item 0.0301309681396463\n", "\\item 0.0118078833271804\n", "\\item 0.0118078833271804\n", "\\item 0.032127101177511\n", "\\item 0.0118078833271804\n", "\\item 0.0374504988094523\n", "\\item 0.0323186543751583\n", "\\item 0.0175916492958266\n", "\\item 0.0118078833271804\n", "\\item 0.0376688357105577\n", "\\item 0.032127101177511\n", "\\item 0.0301309681396461\n", "\\item 0.0351733047983836\n", "\\item 0.0323186543751583\n", "\\item 0.0374504988094523\n", "\\item 0.0351733047983836\n", "\\item 0.0207723954791442\n", "\\item 0.0207723954791442\n", "\\item 0.0207723954791442\n", "\\item 0.0207723954791442\n", "\\item 0.0301309681396463\n", "\\item 0.0207723954791442\n", "\\item 0.0351733047983836\n", "\\item 0.032127101177511\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.0301309681396463\n", "2. 0.0175916492958266\n", "3. 0.0301309681396463\n", "4. 0.0374504988094523\n", "5. 0.0175916492958266\n", "6. 0.0323186543751581\n", "7. 0.00988457403114586\n", "8. 0.0175916492958266\n", "9. 0.0351733047983836\n", "10. 0.0207723954791442\n", "11. 0.0301309681396463\n", "12. 0.0323186543751583\n", "13. 0.0376688357105577\n", "14. 0.0374504988094523\n", "15. 0.0301309681396463\n", "16. 0.0118078833271804\n", "17. 0.0118078833271804\n", "18. 0.032127101177511\n", "19. 0.0118078833271804\n", "20. 0.0374504988094523\n", "21. 0.0323186543751583\n", "22. 0.0175916492958266\n", "23. 0.0118078833271804\n", "24. 0.0376688357105577\n", "25. 0.032127101177511\n", "26. 0.0301309681396461\n", "27. 0.0351733047983836\n", "28. 0.0323186543751583\n", "29. 0.0374504988094523\n", "30. 0.0351733047983836\n", "31. 0.0207723954791442\n", "32. 0.0207723954791442\n", "33. 0.0207723954791442\n", "34. 0.0207723954791442\n", "35. 0.0301309681396463\n", "36. 0.0207723954791442\n", "37. 0.0351733047983836\n", "38. 0.032127101177511\n", "\n", "\n" ], "text/plain": [ " [1] 0.030130968 0.017591649 0.030130968 0.037450499 0.017591649 0.032318654\n", " [7] 0.009884574 0.017591649 0.035173305 0.020772395 0.030130968 0.032318654\n", "[13] 0.037668836 0.037450499 0.030130968 0.011807883 0.011807883 0.032127101\n", "[19] 0.011807883 0.037450499 0.032318654 0.017591649 0.011807883 0.037668836\n", "[25] 0.032127101 0.030130968 0.035173305 0.032318654 0.037450499 0.035173305\n", "[31] 0.020772395 0.020772395 0.020772395 0.020772395 0.030130968 0.020772395\n", "[37] 0.035173305 0.032127101" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ " V1 \n", " Min. :0.009885 \n", " 1st Qu.:0.213697 \n", " Median :0.428511 \n", " Mean :0.482848 \n", " 3rd Qu.:0.704740 \n", " Max. :0.996872 " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "38" ], "text/latex": [ "38" ], "text/markdown": [ "38" ], "text/plain": [ "[1] 38" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "k_snp = which(p.value < 0.05)\n", "p.value[k_snp]\n", "summary(p.value)\n", "length(k_snp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "12. Stworzyć macierz H na podsatwie metody Single-Step zaproponowanej przez VanRadena $\\mathbf{H}^{-1} = \\mathbf{A}^{-1} + \\left(\\begin{array}{cc} 0 & 0 \\\\ 0 & \\mathbf{G}^{-1}-\\mathbf{A}_{22}^{-1} \\end{array} \\right)$" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
1.833333e+00 5.000000e-01 2.775558e-15-6.666667e-01 5.551115e-16-1.00000 -5.662137e-15-1.880440e-15
5.000000e-01 2.000000e+00 5.000000e-01-1.332268e-15-1.000000e+00-1.00000 6.661338e-16-2.844947e-16
1.054712e-15 5.000000e-01 2.000000e+00-3.885781e-16-1.000000e+00 0.50000 -5.551115e-16-1.000000e+00
-6.666667e-01-1.443290e-15-8.881784e-16 2.775710e+01 2.706656e+0127.11523 2.709734e+01 2.706377e+01
-1.110223e-16-1.000000e+00-1.000000e+00 2.706656e+01 2.775375e+0127.08998 2.710357e+01 2.708614e+01
-1.000000e+00-1.000000e+00 5.000000e-01 2.711523e+01 2.708998e+0127.76475 2.706306e+01 2.706698e+01
0.000000e+00-4.440892e-16-6.106227e-16 2.709734e+01 2.710357e+0127.06306 2.770903e+01 2.712700e+01
-1.325329e-15-3.785167e-15-1.000000e+00 2.706377e+01 2.708614e+0127.06698 2.712700e+01 2.775610e+01
\n" ], "text/latex": [ "\\begin{tabular}{llllllll}\n", "\t 1.833333e+00 & 5.000000e-01 & 2.775558e-15 & -6.666667e-01 & 5.551115e-16 & -1.00000 & -5.662137e-15 & -1.880440e-15\\\\\n", "\t 5.000000e-01 & 2.000000e+00 & 5.000000e-01 & -1.332268e-15 & -1.000000e+00 & -1.00000 & 6.661338e-16 & -2.844947e-16\\\\\n", "\t 1.054712e-15 & 5.000000e-01 & 2.000000e+00 & -3.885781e-16 & -1.000000e+00 & 0.50000 & -5.551115e-16 & -1.000000e+00\\\\\n", "\t -6.666667e-01 & -1.443290e-15 & -8.881784e-16 & 2.775710e+01 & 2.706656e+01 & 27.11523 & 2.709734e+01 & 2.706377e+01\\\\\n", "\t -1.110223e-16 & -1.000000e+00 & -1.000000e+00 & 2.706656e+01 & 2.775375e+01 & 27.08998 & 2.710357e+01 & 2.708614e+01\\\\\n", "\t -1.000000e+00 & -1.000000e+00 & 5.000000e-01 & 2.711523e+01 & 2.708998e+01 & 27.76475 & 2.706306e+01 & 2.706698e+01\\\\\n", "\t 0.000000e+00 & -4.440892e-16 & -6.106227e-16 & 2.709734e+01 & 2.710357e+01 & 27.06306 & 2.770903e+01 & 2.712700e+01\\\\\n", "\t -1.325329e-15 & -3.785167e-15 & -1.000000e+00 & 2.706377e+01 & 2.708614e+01 & 27.06698 & 2.712700e+01 & 2.775610e+01\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 1.833333e+00 | 5.000000e-01 | 2.775558e-15 | -6.666667e-01 | 5.551115e-16 | -1.00000 | -5.662137e-15 | -1.880440e-15 |\n", "| 5.000000e-01 | 2.000000e+00 | 5.000000e-01 | -1.332268e-15 | -1.000000e+00 | -1.00000 | 6.661338e-16 | -2.844947e-16 |\n", "| 1.054712e-15 | 5.000000e-01 | 2.000000e+00 | -3.885781e-16 | -1.000000e+00 | 0.50000 | -5.551115e-16 | -1.000000e+00 |\n", "| -6.666667e-01 | -1.443290e-15 | -8.881784e-16 | 2.775710e+01 | 2.706656e+01 | 27.11523 | 2.709734e+01 | 2.706377e+01 |\n", "| -1.110223e-16 | -1.000000e+00 | -1.000000e+00 | 2.706656e+01 | 2.775375e+01 | 27.08998 | 2.710357e+01 | 2.708614e+01 |\n", "| -1.000000e+00 | -1.000000e+00 | 5.000000e-01 | 2.711523e+01 | 2.708998e+01 | 27.76475 | 2.706306e+01 | 2.706698e+01 |\n", "| 0.000000e+00 | -4.440892e-16 | -6.106227e-16 | 2.709734e+01 | 2.710357e+01 | 27.06306 | 2.770903e+01 | 2.712700e+01 |\n", "| -1.325329e-15 | -3.785167e-15 | -1.000000e+00 | 2.706377e+01 | 2.708614e+01 | 27.06698 | 2.712700e+01 | 2.775610e+01 |\n", "\n" ], "text/plain": [ " [,1] [,2] [,3] [,4] [,5] \n", "[1,] 1.833333e+00 5.000000e-01 2.775558e-15 -6.666667e-01 5.551115e-16\n", "[2,] 5.000000e-01 2.000000e+00 5.000000e-01 -1.332268e-15 -1.000000e+00\n", "[3,] 1.054712e-15 5.000000e-01 2.000000e+00 -3.885781e-16 -1.000000e+00\n", "[4,] -6.666667e-01 -1.443290e-15 -8.881784e-16 2.775710e+01 2.706656e+01\n", "[5,] -1.110223e-16 -1.000000e+00 -1.000000e+00 2.706656e+01 2.775375e+01\n", "[6,] -1.000000e+00 -1.000000e+00 5.000000e-01 2.711523e+01 2.708998e+01\n", "[7,] 0.000000e+00 -4.440892e-16 -6.106227e-16 2.709734e+01 2.710357e+01\n", "[8,] -1.325329e-15 -3.785167e-15 -1.000000e+00 2.706377e+01 2.708614e+01\n", " [,6] [,7] [,8] \n", "[1,] -1.00000 -5.662137e-15 -1.880440e-15\n", "[2,] -1.00000 6.661338e-16 -2.844947e-16\n", "[3,] 0.50000 -5.551115e-16 -1.000000e+00\n", "[4,] 27.11523 2.709734e+01 2.706377e+01\n", "[5,] 27.08998 2.710357e+01 2.708614e+01\n", "[6,] 27.76475 2.706306e+01 2.706698e+01\n", "[7,] 27.06306 2.770903e+01 2.712700e+01\n", "[8,] 27.06698 2.712700e+01 2.775610e+01" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "invP = ginv(G) - ginv(A)[4:8, 4:8]\n", "invH = ginv(A) + rbind(cbind(matrix(0, 3, 3), matrix(0, 3, 5)), cbind(matrix(0, 5, 3), invP))\n", "invH" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [], "source": [ "mme4 = function(y, X, Z, invH, sigma_a, sigma_e) {\n", " alpha = sigma_e / sigma_a\n", " C = rbind(cbind(t(X)%*%X, t(X)%*%Z),\n", " cbind(t(Z)%*%X, t(Z)%*%Z+invH*c(alpha)))\n", " rhs = rbind(t(X)%*%y, t(Z)%*%y)\n", " invC = ginv(C)\n", " estimators = invC%*%rhs\n", " list(C = C, est = estimators)\n", "}" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
\n", "\t
$C
\n", "\t\t
\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
3 0 0.000000e+00 0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00 0.00000 1.000000e+00 1.000000e+00
0 2 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.000000e+00 1.00000 0.000000e+00 0.000000e+00
0 0 3.666667e+00 1.000000e+00 5.551115e-15-1.333333e+00 1.110223e-15-2.00000 -1.132427e-14-3.760880e-15
0 0 1.000000e+00 4.000000e+00 1.000000e+00-2.664535e-15-2.000000e+00-2.00000 1.332268e-15-5.689893e-16
0 0 2.109424e-15 1.000000e+00 4.000000e+00-7.771561e-16-2.000000e+00 1.00000 -1.110223e-15-2.000000e+00
1 0 -1.333333e+00-2.886580e-15-1.776357e-15 5.651419e+01 5.413312e+0154.23046 5.419468e+01 5.412755e+01
0 1 -2.220446e-16-2.000000e+00-2.000000e+00 5.413312e+01 5.650749e+0154.17996 5.420714e+01 5.417229e+01
0 1 -2.000000e+00-2.000000e+00 1.000000e+00 5.423046e+01 5.417996e+0156.52950 5.412612e+01 5.413396e+01
1 0 0.000000e+00-8.881784e-16-1.221245e-15 5.419468e+01 5.420714e+0154.12612 5.641806e+01 5.425400e+01
1 0 -2.650657e-15-7.570333e-15-2.000000e+00 5.412755e+01 5.417229e+0154.13396 5.425400e+01 5.651221e+01
\n", "
\n", "\t
$est
\n", "\t\t
\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
4.16814146
3.64849892
0.16057057
-0.31766683
0.11610105
0.27185488
-0.45130721
-0.04569063
-0.27810985
0.50183059
\n", "
\n", "
\n" ], "text/latex": [ "\\begin{description}\n", "\\item[\\$C] \\begin{tabular}{llllllllll}\n", "\t 3 & 0 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 1.000000e+00 & 0.000000e+00 & 0.00000 & 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.00000 & 0.000000e+00 & 0.000000e+00\\\\\n", "\t 0 & 0 & 3.666667e+00 & 1.000000e+00 & 5.551115e-15 & -1.333333e+00 & 1.110223e-15 & -2.00000 & -1.132427e-14 & -3.760880e-15\\\\\n", "\t 0 & 0 & 1.000000e+00 & 4.000000e+00 & 1.000000e+00 & -2.664535e-15 & -2.000000e+00 & -2.00000 & 1.332268e-15 & -5.689893e-16\\\\\n", "\t 0 & 0 & 2.109424e-15 & 1.000000e+00 & 4.000000e+00 & -7.771561e-16 & -2.000000e+00 & 1.00000 & -1.110223e-15 & -2.000000e+00\\\\\n", "\t 1 & 0 & -1.333333e+00 & -2.886580e-15 & -1.776357e-15 & 5.651419e+01 & 5.413312e+01 & 54.23046 & 5.419468e+01 & 5.412755e+01\\\\\n", "\t 0 & 1 & -2.220446e-16 & -2.000000e+00 & -2.000000e+00 & 5.413312e+01 & 5.650749e+01 & 54.17996 & 5.420714e+01 & 5.417229e+01\\\\\n", "\t 0 & 1 & -2.000000e+00 & -2.000000e+00 & 1.000000e+00 & 5.423046e+01 & 5.417996e+01 & 56.52950 & 5.412612e+01 & 5.413396e+01\\\\\n", "\t 1 & 0 & 0.000000e+00 & -8.881784e-16 & -1.221245e-15 & 5.419468e+01 & 5.420714e+01 & 54.12612 & 5.641806e+01 & 5.425400e+01\\\\\n", "\t 1 & 0 & -2.650657e-15 & -7.570333e-15 & -2.000000e+00 & 5.412755e+01 & 5.417229e+01 & 54.13396 & 5.425400e+01 & 5.651221e+01\\\\\n", "\\end{tabular}\n", "\n", "\\item[\\$est] \\begin{tabular}{l}\n", "\t 4.16814146\\\\\n", "\t 3.64849892\\\\\n", "\t 0.16057057\\\\\n", "\t -0.31766683\\\\\n", "\t 0.11610105\\\\\n", "\t 0.27185488\\\\\n", "\t -0.45130721\\\\\n", "\t -0.04569063\\\\\n", "\t -0.27810985\\\\\n", "\t 0.50183059\\\\\n", "\\end{tabular}\n", "\n", "\\end{description}\n" ], "text/markdown": [ "$C\n", ": \n", "| 3 | 0 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 1.000000e+00 | 0.000000e+00 | 0.00000 | 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.00000 | 0.000000e+00 | 0.000000e+00 |\n", "| 0 | 0 | 3.666667e+00 | 1.000000e+00 | 5.551115e-15 | -1.333333e+00 | 1.110223e-15 | -2.00000 | -1.132427e-14 | -3.760880e-15 |\n", "| 0 | 0 | 1.000000e+00 | 4.000000e+00 | 1.000000e+00 | -2.664535e-15 | -2.000000e+00 | -2.00000 | 1.332268e-15 | -5.689893e-16 |\n", "| 0 | 0 | 2.109424e-15 | 1.000000e+00 | 4.000000e+00 | -7.771561e-16 | -2.000000e+00 | 1.00000 | -1.110223e-15 | -2.000000e+00 |\n", "| 1 | 0 | -1.333333e+00 | -2.886580e-15 | -1.776357e-15 | 5.651419e+01 | 5.413312e+01 | 54.23046 | 5.419468e+01 | 5.412755e+01 |\n", "| 0 | 1 | -2.220446e-16 | -2.000000e+00 | -2.000000e+00 | 5.413312e+01 | 5.650749e+01 | 54.17996 | 5.420714e+01 | 5.417229e+01 |\n", "| 0 | 1 | -2.000000e+00 | -2.000000e+00 | 1.000000e+00 | 5.423046e+01 | 5.417996e+01 | 56.52950 | 5.412612e+01 | 5.413396e+01 |\n", "| 1 | 0 | 0.000000e+00 | -8.881784e-16 | -1.221245e-15 | 5.419468e+01 | 5.420714e+01 | 54.12612 | 5.641806e+01 | 5.425400e+01 |\n", "| 1 | 0 | -2.650657e-15 | -7.570333e-15 | -2.000000e+00 | 5.412755e+01 | 5.417229e+01 | 54.13396 | 5.425400e+01 | 5.651221e+01 |\n", "\n", "\n", "$est\n", ": \n", "| 4.16814146 |\n", "| 3.64849892 |\n", "| 0.16057057 |\n", "| -0.31766683 |\n", "| 0.11610105 |\n", "| 0.27185488 |\n", "| -0.45130721 |\n", "| -0.04569063 |\n", "| -0.27810985 |\n", "| 0.50183059 |\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "$C\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.551115e-15 -1.333333e+00\n", " [4,] 0 0 1.000000e+00 4.000000e+00 1.000000e+00 -2.664535e-15\n", " [5,] 0 0 2.109424e-15 1.000000e+00 4.000000e+00 -7.771561e-16\n", " [6,] 1 0 -1.333333e+00 -2.886580e-15 -1.776357e-15 5.651419e+01\n", " [7,] 0 1 -2.220446e-16 -2.000000e+00 -2.000000e+00 5.413312e+01\n", " [8,] 0 1 -2.000000e+00 -2.000000e+00 1.000000e+00 5.423046e+01\n", " [9,] 1 0 0.000000e+00 -8.881784e-16 -1.221245e-15 5.419468e+01\n", "[10,] 1 0 -2.650657e-15 -7.570333e-15 -2.000000e+00 5.412755e+01\n", " [,7] [,8] [,9] [,10]\n", " [1,] 0.000000e+00 0.00000 1.000000e+00 1.000000e+00\n", " [2,] 1.000000e+00 1.00000 0.000000e+00 0.000000e+00\n", " [3,] 1.110223e-15 -2.00000 -1.132427e-14 -3.760880e-15\n", " [4,] -2.000000e+00 -2.00000 1.332268e-15 -5.689893e-16\n", " [5,] -2.000000e+00 1.00000 -1.110223e-15 -2.000000e+00\n", " [6,] 5.413312e+01 54.23046 5.419468e+01 5.412755e+01\n", " [7,] 5.650749e+01 54.17996 5.420714e+01 5.417229e+01\n", " [8,] 5.417996e+01 56.52950 5.412612e+01 5.413396e+01\n", " [9,] 5.420714e+01 54.12612 5.641806e+01 5.425400e+01\n", "[10,] 5.417229e+01 54.13396 5.425400e+01 5.651221e+01\n", "\n", "$est\n", " [,1]\n", " [1,] 4.16814146\n", " [2,] 3.64849892\n", " [3,] 0.16057057\n", " [4,] -0.31766683\n", " [5,] 0.11610105\n", " [6,] 0.27185488\n", " [7,] -0.45130721\n", " [8,] -0.04569063\n", " [9,] -0.27810985\n", "[10,] 0.50183059\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mme4(y, X, Z, invH, 20, 40)" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 0.353771302024288
    2. \n", "\t
  2. 0.529097527052498
    3. \n", "\t
  3. 1.58406895095227
    4. \n", "\t
  4. -0.54690915195531
    5. \n", "\t
  5. -1.89555567544748
    6. \n", "\t
  6. 3.33864133147547
    7. \n", "\t
  7. -0.112859231722328
    8. \n", "\t
  8. 1.87389719543237
    9. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.353771302024288\n", "\\item 0.529097527052498\n", "\\item 1.58406895095227\n", "\\item -0.54690915195531\n", "\\item -1.89555567544748\n", "\\item 3.33864133147547\n", "\\item -0.112859231722328\n", "\\item 1.87389719543237\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.353771302024288\n", "2. 0.529097527052498\n", "3. 1.58406895095227\n", "4. -0.54690915195531\n", "5. -1.89555567544748\n", "6. 3.33864133147547\n", "7. -0.112859231722328\n", "8. 1.87389719543237\n", "\n", "\n" ], "text/plain": [ "[1] 0.3537713 0.5290975 1.5840690 -0.5469092 -1.8955557 3.3386413 -0.1128592\n", "[8] 1.8738972" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 0.594786770216258
    2. \n", "\t
  2. 0.727390903883529
    3. \n", "\t
  3. 1.25859801007004
    4. \n", "\t
  4. 0
    5. \n", "\t
  5. 0
    6. \n", "\t
  6. 1.82719493526976
    7. \n", "\t
  7. 0
    8. \n", "\t
  8. 1.36890364724197
    9. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.594786770216258\n", "\\item 0.727390903883529\n", "\\item 1.25859801007004\n", "\\item 0\n", "\\item 0\n", "\\item 1.82719493526976\n", "\\item 0\n", "\\item 1.36890364724197\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.594786770216258\n", "2. 0.727390903883529\n", "3. 1.25859801007004\n", "4. 0\n", "5. 0\n", "6. 1.82719493526976\n", "7. 0\n", "8. 1.36890364724197\n", "\n", "\n" ], "text/plain": [ "[1] 0.5947868 0.7273909 1.2585980 0.0000000 0.0000000 1.8271949 0.0000000\n", "[8] 1.3689036" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "C = as.matrix(mme4(y, X, Z, invH, 20, 40)$C)\n", "invC = ginv(C)\n", "\n", "invC22 = invC[3:10, 3:10]\n", "(r2 = diag(1 - invC22*2))\n", "r2[r2 < 0] = 0\n", "(r = sqrt(r2))\n" ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [], "source": [ "EM_new = function(y, X, Z, invH, sigma_a, sigma_e) {\n", " n = nrow(X)\n", " p = ncol(X) \n", " q = nrow(invH) \n", " \n", " t = 1 #iteration number 1\n", " tmp = 0.1 #test for convergance\n", " \n", " while (tmp > 0.00001) {\n", " mme_new = mme4(y, X, Z, invH, sigma_a, sigma_e)\n", " C_new = ginv(mme_new$C)\n", " Ck = C_new[(p+1):(p+q), (p+1):(p+q)]\n", " mme2 = mme_new$est\n", " \n", " a = as.matrix(mme2[(p+1):(p+q)])\n", " sigma_a_new = (t(a)%*%invH%*%a + sum(diag(invH%*%Ck))*c(sigma_e))/q\n", " if (sigma_a_new < 0) {\n", " sigma_a_new = 0.01\n", " }\n", " \n", " res = as.matrix(y-X%*%as.matrix(mme2[1:p]) - Z%*%as.matrix(mme2[(p+1):(p+q)]))\n", " X.tmp1 = cbind(X,Z) %*% C_new\n", " X.tmp2 = t(cbind(X,Z))\n", " sigma_e_new = (t(res)%*%res + sum(diag(X.tmp1%*%X.tmp2))*c(sigma_e))/n\n", " if (sigma_e_new < 0) {\n", " sigma_e_new = 0.01\n", " } \n", " tmp = max(abs(sigma_a - sigma_a_new), abs(sigma_e - sigma_e_new))\n", " sigma_a = sigma_a_new\n", " sigma_e = sigma_e_new\n", " \n", " t = t + 1\n", " }\n", " list(t = t, sigma_a = sigma_a, sigma_e = sigma_e)\n", "}" ] }, { "cell_type": "code", "execution_count": 69, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\t
$t
\n", "\t\t
13479
\n", "\t
$sigma_a
\n", "\t\t
\n", "\n", "\t\n", "\n", "
0.004213116
\n", "
\n", "\t
$sigma_e
\n", "\t\t
\n", "\n", "\t\n", "\n", "
0.5440923
\n", "
\n", "
\n" ], "text/latex": [ "\\begin{description}\n", "\\item[\\$t] 13479\n", "\\item[\\$sigma\\_a] \\begin{tabular}{l}\n", "\t 0.004213116\\\\\n", "\\end{tabular}\n", "\n", "\\item[\\$sigma\\_e] \\begin{tabular}{l}\n", "\t 0.5440923\\\\\n", "\\end{tabular}\n", "\n", "\\end{description}\n" ], "text/markdown": [ "$t\n", ": 13479\n", "$sigma_a\n", ": \n", "| 0.004213116 |\n", "\n", "\n", "$sigma_e\n", ": \n", "| 0.5440923 |\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "$t\n", "[1] 13479\n", "\n", "$sigma_a\n", " [,1]\n", "[1,] 0.004213116\n", "\n", "$sigma_e\n", " [,1]\n", "[1,] 0.5440923\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "EM_new(y, X, Z, invH, 1.04, 2.03)" ] } ], "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 }