{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "#### Model mieszany + estymacja parametrów wariancji + dodanie markerów genetycznych" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rozważmy model: " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathbf{y} = \\mathbf{X\\beta+Za+\\epsilon}$, gdzie" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathbf{y}$ - wektor z wartościami fenotypowmi ($n\\times 1$ w naszym przypadku n = ilość zwierząt) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathbf{\\beta}$ - wektor efektów stałych ($p\\times 1$)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathbf{a}$ - wektor efektów losowych ($q\\times 1$)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathbf{\\epsilon}$ - wektor błędów ($n\\times 1$)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathbf{X}$ - macierz wystąpien dla efektów stałych ($n\\times p$)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathbf{Z}$ - macierz wystąpien dla efektów losowych ($n\\times q$)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Zakładamy, że macierze $\\mathbf{X}$ i $\\mathbf{Z}$ są niezależne. Dodatkowo zakładamy, że:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$E(\\mathbf{y}) = \\mathbf{X\\beta}$, $E(\\mathbf{a})=E(\\mathbf{\\epsilon})=0$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$V=Var(\\mathbf{y})=Var(\\mathbf{Za})+ Var(\\mathbf{\\epsilon})=\\mathbf{Z}Var(\\mathbf{a})\\mathbf{Z}^{T}+\\mathbf{I}Var(\\mathbf{\\epsilon})\\mathbf{I}^{T}=\\mathbf{ZGZ}^{T}+\\mathbf{R}$, gdzie" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathbf{G}=\\mathbf{A}\\cdot\\sigma^{2}_{a}$, $\\mathbf{R}=\\mathbf{I}\\cdot\\sigma^{2}_{\\epsilon}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$ \\left[ \\begin{array}{cc}\n", " \\mathbf{X}^{T}\\mathbf{X} & \\mathbf{X}^{T}\\mathbf{Z} \\\\\n", " \\mathbf{Z}^{T}\\mathbf{X} & \\mathbf{Z}^{T}\\mathbf{Z}+\\mathbf{A}^{-1}\\alpha \\end{array}\\right]\\cdot\n", " \\left[ \\begin{array}{c}\n", " \\widehat{\\mathbf{\\beta}} \\\\\n", " \\widehat{\\mathbf{a}} \\end{array}\\right]=\n", " \\left[ \\begin{array}{c}\n", " \\mathbf{X}^{T}\\mathbf{y} \\\\\n", " \\mathbf{Z}^{T}\\mathbf{y} \\end{array}\\right]$, gdzie" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\alpha=\\frac{\\sigma^{2}_{\\epsilon}}{\\sigma^{2}_{a}}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$ \\left[ \\begin{array}{c}\n", " \\widehat{\\mathbf{\\beta}} \\\\\n", " \\widehat{\\mathbf{a}} \\end{array}\\right]=\n", " \\left[ \\begin{array}{cc}\n", " \\mathbf{X}^{T}\\mathbf{X} & \\mathbf{X}^{T}\\mathbf{Z} \\\\\n", " \\mathbf{Z}^{T}\\mathbf{X} & \\mathbf{Z}^{T}\\mathbf{Z}+\\mathbf{A}^{-1}\\alpha \\end{array}\\right]^{-1}\\cdot\n", " \\left[ \\begin{array}{c}\n", " \\mathbf{X}^{T}\\mathbf{y} \\\\\n", " \\mathbf{Z}^{T}\\mathbf{y} \\end{array}\\right]$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$ C = \\left[ \\begin{array}{cc}\n", " \\mathbf{X}^{T}\\mathbf{X} & \\mathbf{X}^{T}\\mathbf{Z} \\\\\n", " \\mathbf{Z}^{T}\\mathbf{X} & \\mathbf{Z}^{T}\\mathbf{Z}+\\mathbf{A}^{-1}\\alpha \\end{array}\\right] = \n", " \\left[ \\begin{array}{cc}\n", " \\mathbf{C}_{11} & \\mathbf{C}_{12} \\\\ \\mathbf{C}_{21} & \\mathbf{C}_{22} \\end{array}\\right]$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$ C^{-1} = \\left[ \\begin{array}{cc}\n", " \\mathbf{C}^{11} & \\mathbf{C}^{12} \\\\ \\mathbf{C}^{21} & \\mathbf{C}^{22} \\end{array}\\right]$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dokładność oszacowania wynosi: $r^2 = diag(1-\\mathbf{C}^{22}\\cdot\\alpha$)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Liczymy macierz A" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "\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": 4, "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" ] }, { "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(Amatrix(ped)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. Definiujemy kolejne źródła informacji" ] }, { "cell_type": "code", "execution_count": 5, "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": 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": 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.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.440093e-15-1.333333e+00 1.554312e-15-2.000000e+00-9.880985e-15-4.996004e-15
0 0 1.000000e+00 4.000000e+00 1.000000e+00-4.440892e-16-2.000000e+00-2.000000e+00-1.110223e-15-6.938894e-17
0 0 2.220446e-15 1.000000e+00 4.000000e+00 5.551115e-16-2.000000e+00 1.000000e+00-5.551115e-16-2.000000e+00
1 0 -1.333333e+00-2.220446e-16-3.108624e-15 4.666667e+00 1.000000e+00-4.440892e-16-2.000000e+00 1.727785e-15
0 1 2.220446e-16-2.000000e+00-2.000000e+00 1.000000e+00 6.000000e+00-1.776357e-15-2.000000e+00-2.636780e-15
0 1 -2.000000e+00-2.000000e+00 1.000000e+00-2.109424e-15-3.330669e-15 6.000000e+00 2.664535e-15-2.000000e+00
1 0 7.771561e-16-3.330669e-15-8.881784e-16-2.000000e+00-2.000000e+00 1.332268e-15 5.000000e+00 1.346145e-15
1 0 -3.497203e-15-8.590351e-15-2.000000e+00-5.280498e-15 2.303713e-15-2.000000e+00 6.231127e-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.440093e-15 & -1.333333e+00 & 1.554312e-15 & -2.000000e+00 & -9.880985e-15 & -4.996004e-15\\\\\n", "\t 0 & 0 & 1.000000e+00 & 4.000000e+00 & 1.000000e+00 & -4.440892e-16 & -2.000000e+00 & -2.000000e+00 & -1.110223e-15 & -6.938894e-17\\\\\n", "\t 0 & 0 & 2.220446e-15 & 1.000000e+00 & 4.000000e+00 & 5.551115e-16 & -2.000000e+00 & 1.000000e+00 & -5.551115e-16 & -2.000000e+00\\\\\n", "\t 1 & 0 & -1.333333e+00 & -2.220446e-16 & -3.108624e-15 & 4.666667e+00 & 1.000000e+00 & -4.440892e-16 & -2.000000e+00 & 1.727785e-15\\\\\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 & -2.636780e-15\\\\\n", "\t 0 & 1 & -2.000000e+00 & -2.000000e+00 & 1.000000e+00 & -2.109424e-15 & -3.330669e-15 & 6.000000e+00 & 2.664535e-15 & -2.000000e+00\\\\\n", "\t 1 & 0 & 7.771561e-16 & -3.330669e-15 & -8.881784e-16 & -2.000000e+00 & -2.000000e+00 & 1.332268e-15 & 5.000000e+00 & 1.346145e-15\\\\\n", "\t 1 & 0 & -3.497203e-15 & -8.590351e-15 & -2.000000e+00 & -5.280498e-15 & 2.303713e-15 & -2.000000e+00 & 6.231127e-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.440093e-15 | -1.333333e+00 | 1.554312e-15 | -2.000000e+00 | -9.880985e-15 | -4.996004e-15 |\n", "| 0 | 0 | 1.000000e+00 | 4.000000e+00 | 1.000000e+00 | -4.440892e-16 | -2.000000e+00 | -2.000000e+00 | -1.110223e-15 | -6.938894e-17 |\n", "| 0 | 0 | 2.220446e-15 | 1.000000e+00 | 4.000000e+00 | 5.551115e-16 | -2.000000e+00 | 1.000000e+00 | -5.551115e-16 | -2.000000e+00 |\n", "| 1 | 0 | -1.333333e+00 | -2.220446e-16 | -3.108624e-15 | 4.666667e+00 | 1.000000e+00 | -4.440892e-16 | -2.000000e+00 | 1.727785e-15 |\n", "| 0 | 1 | 2.220446e-16 | -2.000000e+00 | -2.000000e+00 | 1.000000e+00 | 6.000000e+00 | -1.776357e-15 | -2.000000e+00 | -2.636780e-15 |\n", "| 0 | 1 | -2.000000e+00 | -2.000000e+00 | 1.000000e+00 | -2.109424e-15 | -3.330669e-15 | 6.000000e+00 | 2.664535e-15 | -2.000000e+00 |\n", "| 1 | 0 | 7.771561e-16 | -3.330669e-15 | -8.881784e-16 | -2.000000e+00 | -2.000000e+00 | 1.332268e-15 | 5.000000e+00 | 1.346145e-15 |\n", "| 1 | 0 | -3.497203e-15 | -8.590351e-15 | -2.000000e+00 | -5.280498e-15 | 2.303713e-15 | -2.000000e+00 | 6.231127e-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.440093e-15 -1.333333e+00\n", " [4,] 0 0 1.000000e+00 4.000000e+00 1.000000e+00 -4.440892e-16\n", " [5,] 0 0 2.220446e-15 1.000000e+00 4.000000e+00 5.551115e-16\n", " [6,] 1 0 -1.333333e+00 -2.220446e-16 -3.108624e-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.109424e-15\n", " [9,] 1 0 7.771561e-16 -3.330669e-15 -8.881784e-16 -2.000000e+00\n", "[10,] 1 0 -3.497203e-15 -8.590351e-15 -2.000000e+00 -5.280498e-15\n", " [,7] [,8] [,9] [,10]\n", " [1,] 0.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00\n", " [2,] 1.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00\n", " [3,] 1.554312e-15 -2.000000e+00 -9.880985e-15 -4.996004e-15\n", " [4,] -2.000000e+00 -2.000000e+00 -1.110223e-15 -6.938894e-17\n", " [5,] -2.000000e+00 1.000000e+00 -5.551115e-16 -2.000000e+00\n", " [6,] 1.000000e+00 -4.440892e-16 -2.000000e+00 1.727785e-15\n", " [7,] 6.000000e+00 -1.776357e-15 -2.000000e+00 -2.636780e-15\n", " [8,] -3.330669e-15 6.000000e+00 2.664535e-15 -2.000000e+00\n", " [9,] -2.000000e+00 1.332268e-15 5.000000e+00 1.346145e-15\n", "[10,] 2.303713e-15 -2.000000e+00 6.231127e-15 5.000000e+00\n", "\n", "$est\n", " [,1]\n", " [1,] 4.358502330\n", " [2,] 3.404430006\n", " [3,] 0.098444576\n", " [4,] -0.018770099\n", " [5,] -0.041084203\n", " [6,] -0.008663123\n", " [7,] -0.185732099\n", " [8,] 0.176872088\n", " [9,] -0.249458555\n", "[10,] 0.182614688\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "library(MASS)\n", "\n", "mme = function(y, X, Z, A, sigma_a, sigma_e) {\n", " alpha = sigma_e / sigma_a\n", " invA = ginv(A)\n", " C = rbind(cbind(t(X)%*%X, t(X)%*%Z),\n", " cbind(t(Z)%*%X, t(Z)%*%Z+invA*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", "}\n", "\n", "mme(y, X, Z, A, 20, 40)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "-0.00572209096278567" ], "text/latex": [ "-0.00572209096278567" ], "text/markdown": [ "-0.00572209096278567" ], "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": 19, "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
  3. 99.1697811662677
  4. \n", "\t
  5. 97.7499790740002
  6. \n", "\t
  7. 99.8128679962949
  8. \n", "\t
  9. 88.5463208053009
  10. \n", "\t
  11. 111.618104804649
  12. \n", "\t
  13. 84.4915330653251
  14. \n", "\t
  15. 111.983495048983
  16. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 106.627918039179\n", "\\item 99.1697811662677\n", "\\item 97.7499790740002\n", "\\item 99.8128679962949\n", "\\item 88.5463208053009\n", "\\item 111.618104804649\n", "\\item 84.4915330653251\n", "\\item 111.983495048983\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 106.627918039179\n", "2. 99.1697811662677\n", "3. 97.7499790740002\n", "4. 99.8128679962949\n", "5. 88.5463208053009\n", "6. 111.618104804649\n", "7. 84.4915330653251\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": "markdown", "metadata": {}, "source": [ "4. Dokładność wartości hodowlanych" ] }, { "cell_type": "code", "execution_count": 8, "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", "\n", "
0.59556097 0.15730213 -0.164119413-0.083624565-0.13059083 -0.26455749 -0.14827804 -0.16632621 -0.2842464 -0.2378790
0.15730213 0.80245865 -0.132863260-0.241250738-0.11196840 -0.08730803 -0.29891465 -0.30600266 -0.1859495 -0.1986488
-0.16411941 -0.13286326 0.471094211 0.006928037 0.03264668 0.21954371 0.04495225 0.22077427 0.1386223 0.1341923
-0.08362457 -0.24125074 0.006928037 0.492095721-0.01030797 0.02039033 0.23734577 0.24515571 0.1198194 0.1106640
-0.13059083 -0.11196840 0.032646682-0.010307967 0.45645878 0.04812709 0.20132326 0.02261354 0.1258983 0.2177471
-0.26455749 -0.08730803 0.219543709 0.020390333 0.04812709 0.42768015 0.04704420 0.12757186 0.2428012 0.1231911
-0.14827804 -0.29891465 0.044952254 0.237345770 0.20132326 0.04704420 0.42810675 0.16972255 0.2197160 0.1780739
-0.16632621 -0.30600266 0.220774267 0.245155707 0.02261354 0.12757186 0.16972255 0.44228277 0.1521830 0.2192238
-0.28424641 -0.18594950 0.138622268 0.119819354 0.12589831 0.24280124 0.21971599 0.15218301 0.4418562 0.1680818
-0.23787901 -0.19864885 0.134192262 0.110664009 0.21774710 0.12319108 0.17807393 0.21922376 0.1680818 0.4223641
\n" ], "text/latex": [ "\\begin{tabular}{llllllllll}\n", "\t 0.59556097 & 0.15730213 & -0.164119413 & -0.083624565 & -0.13059083 & -0.26455749 & -0.14827804 & -0.16632621 & -0.2842464 & -0.2378790 \\\\\n", "\t 0.15730213 & 0.80245865 & -0.132863260 & -0.241250738 & -0.11196840 & -0.08730803 & -0.29891465 & -0.30600266 & -0.1859495 & -0.1986488 \\\\\n", "\t -0.16411941 & -0.13286326 & 0.471094211 & 0.006928037 & 0.03264668 & 0.21954371 & 0.04495225 & 0.22077427 & 0.1386223 & 0.1341923 \\\\\n", "\t -0.08362457 & -0.24125074 & 0.006928037 & 0.492095721 & -0.01030797 & 0.02039033 & 0.23734577 & 0.24515571 & 0.1198194 & 0.1106640 \\\\\n", "\t -0.13059083 & -0.11196840 & 0.032646682 & -0.010307967 & 0.45645878 & 0.04812709 & 0.20132326 & 0.02261354 & 0.1258983 & 0.2177471 \\\\\n", "\t -0.26455749 & -0.08730803 & 0.219543709 & 0.020390333 & 0.04812709 & 0.42768015 & 0.04704420 & 0.12757186 & 0.2428012 & 0.1231911 \\\\\n", "\t -0.14827804 & -0.29891465 & 0.044952254 & 0.237345770 & 0.20132326 & 0.04704420 & 0.42810675 & 0.16972255 & 0.2197160 & 0.1780739 \\\\\n", "\t -0.16632621 & -0.30600266 & 0.220774267 & 0.245155707 & 0.02261354 & 0.12757186 & 0.16972255 & 0.44228277 & 0.1521830 & 0.2192238 \\\\\n", "\t -0.28424641 & -0.18594950 & 0.138622268 & 0.119819354 & 0.12589831 & 0.24280124 & 0.21971599 & 0.15218301 & 0.4418562 & 0.1680818 \\\\\n", "\t -0.23787901 & -0.19864885 & 0.134192262 & 0.110664009 & 0.21774710 & 0.12319108 & 0.17807393 & 0.21922376 & 0.1680818 & 0.4223641 \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 0.59556097 | 0.15730213 | -0.164119413 | -0.083624565 | -0.13059083 | -0.26455749 | -0.14827804 | -0.16632621 | -0.2842464 | -0.2378790 |\n", "| 0.15730213 | 0.80245865 | -0.132863260 | -0.241250738 | -0.11196840 | -0.08730803 | -0.29891465 | -0.30600266 | -0.1859495 | -0.1986488 |\n", "| -0.16411941 | -0.13286326 | 0.471094211 | 0.006928037 | 0.03264668 | 0.21954371 | 0.04495225 | 0.22077427 | 0.1386223 | 0.1341923 |\n", "| -0.08362457 | -0.24125074 | 0.006928037 | 0.492095721 | -0.01030797 | 0.02039033 | 0.23734577 | 0.24515571 | 0.1198194 | 0.1106640 |\n", "| -0.13059083 | -0.11196840 | 0.032646682 | -0.010307967 | 0.45645878 | 0.04812709 | 0.20132326 | 0.02261354 | 0.1258983 | 0.2177471 |\n", "| -0.26455749 | -0.08730803 | 0.219543709 | 0.020390333 | 0.04812709 | 0.42768015 | 0.04704420 | 0.12757186 | 0.2428012 | 0.1231911 |\n", "| -0.14827804 | -0.29891465 | 0.044952254 | 0.237345770 | 0.20132326 | 0.04704420 | 0.42810675 | 0.16972255 | 0.2197160 | 0.1780739 |\n", "| -0.16632621 | -0.30600266 | 0.220774267 | 0.245155707 | 0.02261354 | 0.12757186 | 0.16972255 | 0.44228277 | 0.1521830 | 0.2192238 |\n", "| -0.28424641 | -0.18594950 | 0.138622268 | 0.119819354 | 0.12589831 | 0.24280124 | 0.21971599 | 0.15218301 | 0.4418562 | 0.1680818 |\n", "| -0.23787901 | -0.19864885 | 0.134192262 | 0.110664009 | 0.21774710 | 0.12319108 | 0.17807393 | 0.21922376 | 0.1680818 | 0.4223641 |\n", "\n" ], "text/plain": [ " [,1] [,2] [,3] [,4] [,5] [,6] \n", " [1,] 0.59556097 0.15730213 -0.164119413 -0.083624565 -0.13059083 -0.26455749\n", " [2,] 0.15730213 0.80245865 -0.132863260 -0.241250738 -0.11196840 -0.08730803\n", " [3,] -0.16411941 -0.13286326 0.471094211 0.006928037 0.03264668 0.21954371\n", " [4,] -0.08362457 -0.24125074 0.006928037 0.492095721 -0.01030797 0.02039033\n", " [5,] -0.13059083 -0.11196840 0.032646682 -0.010307967 0.45645878 0.04812709\n", " [6,] -0.26455749 -0.08730803 0.219543709 0.020390333 0.04812709 0.42768015\n", " [7,] -0.14827804 -0.29891465 0.044952254 0.237345770 0.20132326 0.04704420\n", " [8,] -0.16632621 -0.30600266 0.220774267 0.245155707 0.02261354 0.12757186\n", " [9,] -0.28424641 -0.18594950 0.138622268 0.119819354 0.12589831 0.24280124\n", "[10,] -0.23787901 -0.19864885 0.134192262 0.110664009 0.21774710 0.12319108\n", " [,7] [,8] [,9] [,10] \n", " [1,] -0.14827804 -0.16632621 -0.2842464 -0.2378790\n", " [2,] -0.29891465 -0.30600266 -0.1859495 -0.1986488\n", " [3,] 0.04495225 0.22077427 0.1386223 0.1341923\n", " [4,] 0.23734577 0.24515571 0.1198194 0.1106640\n", " [5,] 0.20132326 0.02261354 0.1258983 0.2177471\n", " [6,] 0.04704420 0.12757186 0.2428012 0.1231911\n", " [7,] 0.42810675 0.16972255 0.2197160 0.1780739\n", " [8,] 0.16972255 0.44228277 0.1521830 0.2192238\n", " [9,] 0.21971599 0.15218301 0.4418562 0.1680818\n", "[10,] 0.17807393 0.21922376 0.1680818 0.4223641" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 0.0578115770821044
  2. \n", "\t
  3. 0.0158085581151144
  4. \n", "\t
  5. 0.0870824309247233
  6. \n", "\t
  7. 0.144639692852924
  8. \n", "\t
  9. 0.143786506530158
  10. \n", "\t
  11. 0.11543446872744
  12. \n", "\t
  13. 0.116287655050207
  14. \n", "\t
  15. 0.155271707028943
  16. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.0578115770821044\n", "\\item 0.0158085581151144\n", "\\item 0.0870824309247233\n", "\\item 0.144639692852924\n", "\\item 0.143786506530158\n", "\\item 0.11543446872744\n", "\\item 0.116287655050207\n", "\\item 0.155271707028943\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.0578115770821044\n", "2. 0.0158085581151144\n", "3. 0.0870824309247233\n", "4. 0.144639692852924\n", "5. 0.143786506530158\n", "6. 0.11543446872744\n", "7. 0.116287655050207\n", "8. 0.155271707028943\n", "\n", "\n" ], "text/plain": [ "[1] 0.05781158 0.01580856 0.08708243 0.14463969 0.14378651 0.11543447 0.11628766\n", "[8] 0.15527171" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 0.240440381554564
  2. \n", "\t
  3. 0.125732088645319
  4. \n", "\t
  5. 0.295097324496044
  6. \n", "\t
  7. 0.38031525456248
  8. \n", "\t
  9. 0.379191912532635
  10. \n", "\t
  11. 0.33975648445238
  12. \n", "\t
  13. 0.341009757998517
  14. \n", "\t
  15. 0.394045310883079
  16. \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": [ "5. Istotność efektów stałych" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Test Walda:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$W = \\frac{\\widehat{\\beta}}{se(\\widehat{\\beta})}\\sim\\mathcal{N}(0,1)$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$H_{0}: \\beta=0$ vs. $H_{1}: \\beta\\neq 0$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$var(\\widehat{\\beta}) = \\left(\\mathbf{X}^{T}\\mathbf{V}^{-1}\\mathbf{X}\\right)^{-1}$" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "scrolled": false }, "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
  3. 5.66554023294585
  4. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 4.88082357807458\n", "\\item 5.66554023294585\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 4.88082357807458\n", "2. 5.66554023294585\n", "\n", "\n" ], "text/plain": [ "[1] 4.880824 5.665540" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 0.892985017822407
  2. \n", "\t
  3. 0.600901214346598
  4. \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.371865195985295
  2. \n", "\t
  3. 0.547905784351094
  4. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.371865195985295\n", "\\item 0.547905784351094\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.371865195985295\n", "2. 0.547905784351094\n", "\n", "\n" ], "text/plain": [ "[1] 0.3718652 0.5479058" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "G = A*20\n", "R = diag(5)*40\n", "V = Z%*%G%*%t(Z) + R\n", "V\n", "\n", "(varB = ginv(t(X)%*%ginv(V)%*%X))\n", "(seB = sqrt(diag(varB)))\n", "\n", "(testWalda = mme(y, X, Z, A, 20, 40)$est[1:2] / seB)\n", "(p_value = 2*pnorm(abs(testWalda), lower.tail = FALSE))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "6. Estymacja parametrów wariancji" ] }, { "cell_type": "code", "execution_count": 10, "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": 11, "metadata": {}, "outputs": [], "source": [ "sigma_a = 100.01 #starting value for random effect\n", "sigma_e = 9.01 #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": [ "#### $\\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": 14, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
\n", "\t
$t
\n", "\t\t
657
\n", "\t
$sigma_a
\n", "\t\t
\n", "\n", "\t\n", "\n", "
0.6694952
\n", "
\n", "\t
$sigma_e
\n", "\t\t
\n", "\n", "\t\n", "\n", "
0.005739873
\n", "
\n", "
\n" ], "text/latex": [ "\\begin{description}\n", "\\item[\\$t] 657\n", "\\item[\\$sigma\\_a] \\begin{tabular}{l}\n", "\t 0.6694952\\\\\n", "\\end{tabular}\n", "\n", "\\item[\\$sigma\\_e] \\begin{tabular}{l}\n", "\t 0.005739873\\\\\n", "\\end{tabular}\n", "\n", "\\end{description}\n" ], "text/markdown": [ "$t\n", ": 657\n", "$sigma_a\n", ": \n", "| 0.6694952 |\n", "\n", "\n", "$sigma_e\n", ": \n", "| 0.005739873 |\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "$t\n", "[1] 657\n", "\n", "$sigma_a\n", " [,1]\n", "[1,] 0.6694952\n", "\n", "$sigma_e\n", " [,1]\n", "[1,] 0.005739873\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "\t\n", "\n", "
0.675235
\n" ], "text/latex": [ "\\begin{tabular}{l}\n", "\t 0.675235\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 0.675235 |\n", "\n" ], "text/plain": [ " [,1] \n", "[1,] 0.675235" ] }, "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": [ "7. Odziedziczalność" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ "0.970839380554511" ], "text/latex": [ "0.970839380554511" ], "text/markdown": [ "0.970839380554511" ], "text/plain": [ "[1] 0.9708394" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "(h2 = 0.6555448 / (0.6694952 + 0.005739873))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "scrolled": true }, "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", "\n", "
4.44171842
3.48850299
0.28002080
-0.13721162
-0.08518631
0.05357608
-0.58605426
0.40904828
-0.93033132
0.55159998
\n" ], "text/latex": [ "\\begin{tabular}{l}\n", "\t 4.44171842\\\\\n", "\t 3.48850299\\\\\n", "\t 0.28002080\\\\\n", "\t -0.13721162\\\\\n", "\t -0.08518631\\\\\n", "\t 0.05357608\\\\\n", "\t -0.58605426\\\\\n", "\t 0.40904828\\\\\n", "\t -0.93033132\\\\\n", "\t 0.55159998\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 4.44171842 |\n", "| 3.48850299 |\n", "| 0.28002080 |\n", "| -0.13721162 |\n", "| -0.08518631 |\n", "| 0.05357608 |\n", "| -0.58605426 |\n", "| 0.40904828 |\n", "| -0.93033132 |\n", "| 0.55159998 |\n", "\n" ], "text/plain": [ " [,1] \n", " [1,] 4.44171842\n", " [2,] 3.48850299\n", " [3,] 0.28002080\n", " [4,] -0.13721162\n", " [5,] -0.08518631\n", " [6,] 0.05357608\n", " [7,] -0.58605426\n", " [8,] 0.40904828\n", " [9,] -0.93033132\n", "[10,] 0.55159998" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mme(y, X, Z, A, wyniki$sigma_a, wyniki$sigma_e)$est" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "-0.0555672951330806" ], "text/latex": [ "-0.0555672951330806" ], "text/markdown": [ "-0.0555672951330806" ], "text/plain": [ "[1] -0.0555673" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "0.501885441017205" ], "text/latex": [ "0.501885441017205" ], "text/markdown": [ "0.501885441017205" ], "text/plain": [ "[1] 0.5018854" ] }, "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.37325
97.74998 99.40985
99.81287102.17467
88.54632 89.43012
111.61810109.25740
84.49153 82.57044
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.37325\\\\\n", "\t 97.74998 & 99.40985\\\\\n", "\t 99.81287 & 102.17467\\\\\n", "\t 88.54632 & 89.43012\\\\\n", "\t 111.61810 & 109.25740\\\\\n", "\t 84.49153 & 82.57044\\\\\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.37325 |\n", "| 97.74998 | 99.40985 |\n", "| 99.81287 | 102.17467 |\n", "| 88.54632 | 89.43012 |\n", "| 111.61810 | 109.25740 |\n", "| 84.49153 | 82.57044 |\n", "| 111.98350 | 112.09773 |\n", "\n" ], "text/plain": [ " est est2 \n", "[1,] 106.62792 106.68655\n", "[2,] 99.16978 98.37325\n", "[3,] 97.74998 99.40985\n", "[4,] 99.81287 102.17467\n", "[5,] 88.54632 89.43012\n", "[6,] 111.61810 109.25740\n", "[7,] 84.49153 82.57044\n", "[8,] 111.98350 112.09773" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "est2 = mme(y, X, Z, A, wyniki$sigma_a, wyniki$sigma_e)$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": "markdown", "metadata": {}, "source": [ "8. Test Walda z nowymi $\\sigma^{2}_{a}$ i $\\sigma^{2}_{\\epsilon}$" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
0.67523500.00000000.16737380.33474760.0836869
0.00000000.67523500.16737380.33474760.2510607
0.16737380.16737380.67523500.16737380.3347476
0.33474760.33474760.16737380.67523500.1673738
0.08368690.25106070.33474760.16737380.6752350
\n" ], "text/latex": [ "\\begin{tabular}{lllll}\n", "\t 0.6752350 & 0.0000000 & 0.1673738 & 0.3347476 & 0.0836869\\\\\n", "\t 0.0000000 & 0.6752350 & 0.1673738 & 0.3347476 & 0.2510607\\\\\n", "\t 0.1673738 & 0.1673738 & 0.6752350 & 0.1673738 & 0.3347476\\\\\n", "\t 0.3347476 & 0.3347476 & 0.1673738 & 0.6752350 & 0.1673738\\\\\n", "\t 0.0836869 & 0.2510607 & 0.3347476 & 0.1673738 & 0.6752350\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 0.6752350 | 0.0000000 | 0.1673738 | 0.3347476 | 0.0836869 |\n", "| 0.0000000 | 0.6752350 | 0.1673738 | 0.3347476 | 0.2510607 |\n", "| 0.1673738 | 0.1673738 | 0.6752350 | 0.1673738 | 0.3347476 |\n", "| 0.3347476 | 0.3347476 | 0.1673738 | 0.6752350 | 0.1673738 |\n", "| 0.0836869 | 0.2510607 | 0.3347476 | 0.1673738 | 0.6752350 |\n", "\n" ], "text/plain": [ " [,1] [,2] [,3] [,4] [,5] \n", "[1,] 0.6752350 0.0000000 0.1673738 0.3347476 0.0836869\n", "[2,] 0.0000000 0.6752350 0.1673738 0.3347476 0.2510607\n", "[3,] 0.1673738 0.1673738 0.6752350 0.1673738 0.3347476\n", "[4,] 0.3347476 0.3347476 0.1673738 0.6752350 0.1673738\n", "[5,] 0.0836869 0.2510607 0.3347476 0.1673738 0.6752350" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\n", "
0.34510400.2065156
0.20651560.3626311
\n" ], "text/latex": [ "\\begin{tabular}{ll}\n", "\t 0.3451040 & 0.2065156\\\\\n", "\t 0.2065156 & 0.3626311\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 0.3451040 | 0.2065156 |\n", "| 0.2065156 | 0.3626311 |\n", "\n" ], "text/plain": [ " [,1] [,2] \n", "[1,] 0.3451040 0.2065156\n", "[2,] 0.2065156 0.3626311" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 0.587455495770437
  2. \n", "\t
  3. 0.602188575392395
  4. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.587455495770437\n", "\\item 0.602188575392395\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.587455495770437\n", "2. 0.602188575392395\n", "\n", "\n" ], "text/plain": [ "[1] 0.5874555 0.6021886" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 7.5609445323216
  2. \n", "\t
  3. 5.79304080974346
  4. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 7.5609445323216\n", "\\item 5.79304080974346\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 7.5609445323216\n", "2. 5.79304080974346\n", "\n", "\n" ], "text/plain": [ "[1] 7.560945 5.793041" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 4.00152869669226e-14
  2. \n", "\t
  3. 6.91233023736263e-09
  4. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 4.00152869669226e-14\n", "\\item 6.91233023736263e-09\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 4.00152869669226e-14\n", "2. 6.91233023736263e-09\n", "\n", "\n" ], "text/plain": [ "[1] 4.001529e-14 6.912330e-09" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "G = A*c(wyniki$sigma_a)\n", "R = diag(5)*c(wyniki$sigma_e)\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, wyniki$sigma_a, wyniki$sigma_e)$est[1:2] / seB)\n", "(p_value = 2*pnorm(abs(testWalda), lower.tail = FALSE))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "9. Dokładność oceny" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 0.189935890368132
  2. \n", "\t
  3. 0.0776482185901066
  4. \n", "\t
  5. 0.291998383422087
  6. \n", "\t
  7. 0.482418666354055
  8. \n", "\t
  9. 0.459808525687841
  10. \n", "\t
  11. 0.457012381382101
  12. \n", "\t
  13. 0.479622522048315
  14. \n", "\t
  15. 0.483237179540816
  16. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.189935890368132\n", "\\item 0.0776482185901066\n", "\\item 0.291998383422087\n", "\\item 0.482418666354055\n", "\\item 0.459808525687841\n", "\\item 0.457012381382101\n", "\\item 0.479622522048315\n", "\\item 0.483237179540816\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.189935890368132\n", "2. 0.0776482185901066\n", "3. 0.291998383422087\n", "4. 0.482418666354055\n", "5. 0.459808525687841\n", "6. 0.457012381382101\n", "7. 0.479622522048315\n", "8. 0.483237179540816\n", "\n", "\n" ], "text/plain": [ "[1] 0.18993589 0.07764822 0.29199838 0.48241867 0.45980853 0.45701238 0.47962252\n", "[8] 0.48323718" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 0.435816349358457
  2. \n", "\t
  3. 0.278654299428713
  4. \n", "\t
  5. 0.540368747636359
  6. \n", "\t
  7. 0.694563651765665
  8. \n", "\t
  9. 0.678091826884708
  10. \n", "\t
  11. 0.676026908770724
  12. \n", "\t
  13. 0.692547848201346
  14. \n", "\t
  15. 0.695152630391928
  16. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.435816349358457\n", "\\item 0.278654299428713\n", "\\item 0.540368747636359\n", "\\item 0.694563651765665\n", "\\item 0.678091826884708\n", "\\item 0.676026908770724\n", "\\item 0.692547848201346\n", "\\item 0.695152630391928\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.435816349358457\n", "2. 0.278654299428713\n", "3. 0.540368747636359\n", "4. 0.694563651765665\n", "5. 0.678091826884708\n", "6. 0.676026908770724\n", "7. 0.692547848201346\n", "8. 0.695152630391928\n", "\n", "\n" ], "text/plain": [ "[1] 0.4358163 0.2786543 0.5403687 0.6945637 0.6780918 0.6760269 0.6925478\n", "[8] 0.6951526" ] }, "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.4358163
0.12573210.2786543
0.29509730.5403687
0.38031530.6945637
0.37919190.6780918
0.33975650.6760269
0.34100980.6925478
0.39404530.6951526
\n" ], "text/latex": [ "\\begin{tabular}{ll}\n", " r & r\\_2\\\\\n", "\\hline\n", "\t 0.2404404 & 0.4358163\\\\\n", "\t 0.1257321 & 0.2786543\\\\\n", "\t 0.2950973 & 0.5403687\\\\\n", "\t 0.3803153 & 0.6945637\\\\\n", "\t 0.3791919 & 0.6780918\\\\\n", "\t 0.3397565 & 0.6760269\\\\\n", "\t 0.3410098 & 0.6925478\\\\\n", "\t 0.3940453 & 0.6951526\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| r | r_2 |\n", "|---|---|\n", "| 0.2404404 | 0.4358163 |\n", "| 0.1257321 | 0.2786543 |\n", "| 0.2950973 | 0.5403687 |\n", "| 0.3803153 | 0.6945637 |\n", "| 0.3791919 | 0.6780918 |\n", "| 0.3397565 | 0.6760269 |\n", "| 0.3410098 | 0.6925478 |\n", "| 0.3940453 | 0.6951526 |\n", "\n" ], "text/plain": [ " r r_2 \n", "[1,] 0.2404404 0.4358163\n", "[2,] 0.1257321 0.2786543\n", "[3,] 0.2950973 0.5403687\n", "[4,] 0.3803153 0.6945637\n", "[5,] 0.3791919 0.6780918\n", "[6,] 0.3397565 0.6760269\n", "[7,] 0.3410098 0.6925478\n", "[8,] 0.3940453 0.6951526" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "C = as.matrix(mme(y, X, Z, A, wyniki$sigma_a, wyniki$sigma_e)$C)\n", "invC = ginv(C)\n", "\n", "invC22 = invC[3:10, 3:10]\n", "alpha = wyniki$sigma_e / wyniki$sigma_a\n", "(r_22 = diag(1 - invC22*c(alpha)))\n", "(r_2 = sqrt(r_22))\n", "\n", "cbind(r, r_2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "10. Korelacje oceny" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "0.986408003048259" ], "text/latex": [ "0.986408003048259" ], "text/markdown": [ "0.986408003048259" ], "text/plain": [ "[1] 0.986408" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "0.977062095530124" ], "text/latex": [ "0.977062095530124" ], "text/markdown": [ "0.977062095530124" ], "text/plain": [ "[1] 0.9770621" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cor(est, est2)\n", "cor(r, r_2)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "scrolled": true }, "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": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\n", "Call:\n", "lm(formula = est ~ est2)\n", "\n", "Residuals:\n", " Min 1Q Median 3Q Max \n", "-2.33224 -1.18757 0.04123 1.00186 2.48653 \n", "\n", "Coefficients:\n", " Estimate Std. Error t value Pr(>|t|) \n", "(Intercept) 1.35920 6.73740 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": [ "11. Macierz z genotypami" ] }, { "cell_type": "code", "execution_count": 211, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 'BB'
  2. \n", "\t
  3. 'AB'
  4. \n", "\t
  5. 'BB'
  6. \n", "\t
  7. 'AA'
  8. \n", "\t
  9. 'AA'
  10. \n", "\t
  11. 'AA'
  12. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 'BB'\n", "\\item 'AB'\n", "\\item 'BB'\n", "\\item 'AA'\n", "\\item 'AA'\n", "\\item 'AA'\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 'BB'\n", "2. 'AB'\n", "3. 'BB'\n", "4. 'AA'\n", "5. 'AA'\n", "6. 'AA'\n", "\n", "\n" ], "text/plain": [ "[1] \"BB\" \"AB\" \"BB\" \"AA\" \"AA\" \"AA\"" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Genotypes = sample(c(\"AA\", \"AB\", \"BB\"), 5000, replace = TRUE)\n", "head(Genotypes)" ] }, { "cell_type": "code", "execution_count": 212, "metadata": {}, "outputs": [], "source": [ "Genotypes = as.data.frame(matrix(Genotypes, 5, 1000))" ] }, { "cell_type": "code", "execution_count": 213, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
V1V2V3V4V5V6V7V8V9V10...V991V992V993V994V995V996V997V998V999V1000
BB AA BB AA BB AA AA AB AB BB ...BB AB BB AA AA AB AB AB AB AB
AB AA BB BB AB BB AA AB AA BB ...AA BB BB AA AA AA BB BB BB AA
BB AA AA AA BB AA BB AA BB AA ...AB AA AA BB BB AB AA AB AB AB
AA AB BB BB BB BB AB BB AB BB ...AA AA BB BB AB AB AA AB BB AA
AA AB AA BB AA AA AA BB AA BB ...BB AB AA BB BB AB AA BB BB AA
\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 BB & AA & BB & AA & BB & AA & AA & AB & AB & BB & ... & BB & AB & BB & AA & AA & AB & AB & AB & AB & AB \\\\\n", "\t AB & AA & BB & BB & AB & BB & AA & AB & AA & BB & ... & AA & BB & BB & AA & AA & AA & BB & BB & BB & AA \\\\\n", "\t BB & AA & AA & AA & BB & AA & BB & AA & BB & AA & ... & AB & AA & AA & BB & BB & AB & AA & AB & AB & AB \\\\\n", "\t AA & AB & BB & BB & BB & BB & AB & BB & AB & BB & ... & AA & AA & BB & BB & AB & AB & AA & AB & BB & AA \\\\\n", "\t AA & AB & AA & BB & AA & AA & AA & BB & AA & BB & ... & BB & AB & AA & BB & BB & AB & AA & BB & BB & AA \\\\\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", "| BB | AA | BB | AA | BB | AA | AA | AB | AB | BB | ... | BB | AB | BB | AA | AA | AB | AB | AB | AB | AB |\n", "| AB | AA | BB | BB | AB | BB | AA | AB | AA | BB | ... | AA | BB | BB | AA | AA | AA | BB | BB | BB | AA |\n", "| BB | AA | AA | AA | BB | AA | BB | AA | BB | AA | ... | AB | AA | AA | BB | BB | AB | AA | AB | AB | AB |\n", "| AA | AB | BB | BB | BB | BB | AB | BB | AB | BB | ... | AA | AA | BB | BB | AB | AB | AA | AB | BB | AA |\n", "| AA | AB | AA | BB | AA | AA | AA | BB | AA | BB | ... | BB | AB | AA | BB | BB | AB | AA | BB | BB | AA |\n", "\n" ], "text/plain": [ " V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 ... V991 V992 V993 V994 V995 V996 V997 V998\n", "1 BB AA BB AA BB AA AA AB AB BB ... BB AB BB AA AA AB AB AB \n", "2 AB AA BB BB AB BB AA AB AA BB ... AA BB BB AA AA AA BB BB \n", "3 BB AA AA AA BB AA BB AA BB AA ... AB AA AA BB BB AB AA AB \n", "4 AA AB BB BB BB BB AB BB AB BB ... AA AA BB BB AB AB AA AB \n", "5 AA AB AA BB AA AA AA BB AA BB ... BB AB AA BB BB AB AA BB \n", " V999 V1000\n", "1 AB AB \n", "2 BB AA \n", "3 AB AB \n", "4 BB AA \n", "5 BB AA " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "head(Genotypes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "11. Macierz G" ] }, { "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": 214, "metadata": {}, "outputs": [], "source": [ "tmp = data.frame(matrix(0, 3, 1))\n", "colnames(tmp) = c(\"Var1\")\n", "tmp[, 1] = c(\"AA\", \"AB\", \"BB\")\n", "tmp = merge(tmp, as.data.frame(table(Genotypes[, 1])), by = \"Var1\", all.x = TRUE)\n", "tmp[is.na(tmp)] = 0" ] }, { "cell_type": "code", "execution_count": 215, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\n", "
Var1Freq
AA2
AB1
BB2
\n" ], "text/latex": [ "\\begin{tabular}{r|ll}\n", " Var1 & Freq\\\\\n", "\\hline\n", "\t AA & 2 \\\\\n", "\t AB & 1 \\\\\n", "\t BB & 2 \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| Var1 | Freq |\n", "|---|---|\n", "| AA | 2 |\n", "| AB | 1 |\n", "| BB | 2 |\n", "\n" ], "text/plain": [ " Var1 Freq\n", "1 AA 2 \n", "2 AB 1 \n", "3 BB 2 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "tmp" ] }, { "cell_type": "code", "execution_count": 216, "metadata": {}, "outputs": [], "source": [ "p = matrix(0, 3, 1000)\n", "for (i in 1:1000) {\n", " tmp = data.frame(matrix(0, 3, 1))\n", " colnames(tmp) = c(\"Var1\")\n", " tmp[, 1] = c(\"AA\", \"AB\", \"BB\")\n", " tmp = merge(tmp, as.data.frame(table(Genotypes[, i])), by = \"Var1\", all.x = TRUE)\n", " tmp[is.na(tmp)] = 0\n", " p[, i] = tmp$Freq\n", "}" ] }, { "cell_type": "code", "execution_count": 217, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\n", "
2 3 2 2 1 3 3 1 2 1 ...2 2 2 2 2 1 3 0 0 3
1 2 0 0 1 0 1 2 2 0 ...1 2 0 0 1 4 1 3 2 2
2 0 3 3 3 2 1 2 1 4 ...2 1 3 3 2 0 1 2 3 0
\n" ], "text/latex": [ "\\begin{tabular}{llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll}\n", "\t 2 & 3 & 2 & 2 & 1 & 3 & 3 & 1 & 2 & 1 & ... & 2 & 2 & 2 & 2 & 2 & 1 & 3 & 0 & 0 & 3 \\\\\n", "\t 1 & 2 & 0 & 0 & 1 & 0 & 1 & 2 & 2 & 0 & ... & 1 & 2 & 0 & 0 & 1 & 4 & 1 & 3 & 2 & 2 \\\\\n", "\t 2 & 0 & 3 & 3 & 3 & 2 & 1 & 2 & 1 & 4 & ... & 2 & 1 & 3 & 3 & 2 & 0 & 1 & 2 & 3 & 0 \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 2 | 3 | 2 | 2 | 1 | 3 | 3 | 1 | 2 | 1 | ... | 2 | 2 | 2 | 2 | 2 | 1 | 3 | 0 | 0 | 3 |\n", "| 1 | 2 | 0 | 0 | 1 | 0 | 1 | 2 | 2 | 0 | ... | 1 | 2 | 0 | 0 | 1 | 4 | 1 | 3 | 2 | 2 |\n", "| 2 | 0 | 3 | 3 | 3 | 2 | 1 | 2 | 1 | 4 | ... | 2 | 1 | 3 | 3 | 2 | 0 | 1 | 2 | 3 | 0 |\n", "\n" ], "text/plain": [ " [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]\n", "[1,] 2 3 2 2 1 3 3 1 2 1 ... 2 2 2 \n", "[2,] 1 2 0 0 1 0 1 2 2 0 ... 1 2 0 \n", "[3,] 2 0 3 3 3 2 1 2 1 4 ... 2 1 3 \n", " [,15] [,16] [,17] [,18] [,19] [,20] [,21]\n", "[1,] 2 2 1 3 0 0 3 \n", "[2,] 0 1 4 1 3 2 2 \n", "[3,] 3 2 0 1 2 3 0 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p" ] }, { "cell_type": "code", "execution_count": 218, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 3
  2. \n", "\t
  3. 1000
  4. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 3\n", "\\item 1000\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 3\n", "2. 1000\n", "\n", "\n" ], "text/plain": [ "[1] 3 1000" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p = as.data.frame(p)\n", "dim(p)" ] }, { "cell_type": "code", "execution_count": 220, "metadata": {}, "outputs": [], "source": [ "colMax <- function(data) sapply(data, max, na.rm = TRUE)\n", "k_max = which(colMax(p) == 5)" ] }, { "cell_type": "code", "execution_count": 221, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\t
V36
\n", "\t\t
36
\n", "\t
V58
\n", "\t\t
58
\n", "\t
V108
\n", "\t\t
108
\n", "\t
V131
\n", "\t\t
131
\n", "\t
V186
\n", "\t\t
186
\n", "\t
V278
\n", "\t\t
278
\n", "\t
V510
\n", "\t\t
510
\n", "\t
V573
\n", "\t\t
573
\n", "\t
V751
\n", "\t\t
751
\n", "\t
V850
\n", "\t\t
850
\n", "
\n" ], "text/latex": [ "\\begin{description*}\n", "\\item[V36] 36\n", "\\item[V58] 58\n", "\\item[V108] 108\n", "\\item[V131] 131\n", "\\item[V186] 186\n", "\\item[V278] 278\n", "\\item[V510] 510\n", "\\item[V573] 573\n", "\\item[V751] 751\n", "\\item[V850] 850\n", "\\end{description*}\n" ], "text/markdown": [ "V36\n", ": 36V58\n", ": 58V108\n", ": 108V131\n", ": 131V186\n", ": 186V278\n", ": 278V510\n", ": 510V573\n", ": 573V751\n", ": 751V850\n", ": 850\n", "\n" ], "text/plain": [ " V36 V58 V108 V131 V186 V278 V510 V573 V751 V850 \n", " 36 58 108 131 186 278 510 573 751 850 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "k_max" ] }, { "cell_type": "code", "execution_count": 222, "metadata": {}, "outputs": [], "source": [ "Genotypes = Genotypes[,-k_max]" ] }, { "cell_type": "code", "execution_count": 223, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 5
  2. \n", "\t
  3. 990
  4. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 5\n", "\\item 990\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 5\n", "2. 990\n", "\n", "\n" ], "text/plain": [ "[1] 5 990" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "dim(Genotypes)" ] }, { "cell_type": "code", "execution_count": 224, "metadata": {}, "outputs": [], "source": [ "p2 = numeric(ncol(Genotypes))\n", "for (i in 1:ncol(Genotypes)) {\n", " p2[i] = (2*p[1, i] + p[2, i]) / 10\n", "}" ] }, { "cell_type": "code", "execution_count": 225, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " Min. 1st Qu. Median Mean 3rd Qu. Max. \n", " 0.0000 0.4000 0.5000 0.4937 0.6000 1.0000 " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "990" ], "text/latex": [ "990" ], "text/markdown": [ "990" ], "text/plain": [ "[1] 990" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "summary(p2)\n", "length(p2)" ] }, { "cell_type": "code", "execution_count": 226, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 58
  2. \n", "\t
  3. 510
  4. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 58\n", "\\item 510\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 58\n", "2. 510\n", "\n", "\n" ], "text/plain": [ "[1] 58 510" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "(k = which(p2 == 0))\n", "p2 = p2[-k]" ] }, { "cell_type": "code", "execution_count": 227, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 107
  2. \n", "\t
  3. 130
  4. \n", "\t
  5. 571
  6. \n", "\t
  7. 749
  8. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 107\n", "\\item 130\n", "\\item 571\n", "\\item 749\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 107\n", "2. 130\n", "3. 571\n", "4. 749\n", "\n", "\n" ], "text/plain": [ "[1] 107 130 571 749" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "(k2 = which(p2 == 1))\n", "p2 = p2[-k2]" ] }, { "cell_type": "code", "execution_count": 228, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 5
  2. \n", "\t
  3. 984
  4. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 5\n", "\\item 984\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 5\n", "2. 984\n", "\n", "\n" ], "text/plain": [ "[1] 5 984" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Genotypes2 = Genotypes[, -c(k, k2)]\n", "dim(Genotypes2)" ] }, { "cell_type": "code", "execution_count": 229, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
V1V2V3V4V5V6V7V8V9V10...V991V992V993V994V995V996V997V998V999V1000
BB AA BB AA BB AA AA AB AB BB ...BB AB BB AA AA AB AB AB AB AB
AB AA BB BB AB BB AA AB AA BB ...AA BB BB AA AA AA BB BB BB AA
BB AA AA AA BB AA BB AA BB AA ...AB AA AA BB BB AB AA AB AB AB
AA AB BB BB BB BB AB BB AB BB ...AA AA BB BB AB AB AA AB BB AA
AA AB AA BB AA AA AA BB AA BB ...BB AB AA BB BB AB AA BB BB AA
\n" ], "text/latex": [ "\\begin{tabular}{r|llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll}\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 BB & AA & BB & AA & BB & AA & AA & AB & AB & BB & ... & BB & AB & BB & AA & AA & AB & AB & AB & AB & AB \\\\\n", "\t AB & AA & BB & BB & AB & BB & AA & AB & AA & BB & ... & AA & BB & BB & AA & AA & AA & BB & BB & BB & AA \\\\\n", "\t BB & AA & AA & AA & BB & AA & BB & AA & BB & AA & ... & AB & AA & AA & BB & BB & AB & AA & AB & AB & AB \\\\\n", "\t AA & AB & BB & BB & BB & BB & AB & BB & AB & BB & ... & AA & AA & BB & BB & AB & AB & AA & AB & BB & AA \\\\\n", "\t AA & AB & AA & BB & AA & AA & AA & BB & AA & BB & ... & BB & AB & AA & BB & BB & AB & AA & BB & BB & AA \\\\\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", "| BB | AA | BB | AA | BB | AA | AA | AB | AB | BB | ... | BB | AB | BB | AA | AA | AB | AB | AB | AB | AB |\n", "| AB | AA | BB | BB | AB | BB | AA | AB | AA | BB | ... | AA | BB | BB | AA | AA | AA | BB | BB | BB | AA |\n", "| BB | AA | AA | AA | BB | AA | BB | AA | BB | AA | ... | AB | AA | AA | BB | BB | AB | AA | AB | AB | AB |\n", "| AA | AB | BB | BB | BB | BB | AB | BB | AB | BB | ... | AA | AA | BB | BB | AB | AB | AA | AB | BB | AA |\n", "| AA | AB | AA | BB | AA | AA | AA | BB | AA | BB | ... | BB | AB | AA | BB | BB | AB | AA | BB | BB | AA |\n", "\n" ], "text/plain": [ " V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 ... V991 V992 V993 V994 V995 V996 V997 V998\n", "1 BB AA BB AA BB AA AA AB AB BB ... BB AB BB AA AA AB AB AB \n", "2 AB AA BB BB AB BB AA AB AA BB ... AA BB BB AA AA AA BB BB \n", "3 BB AA AA AA BB AA BB AA BB AA ... AB AA AA BB BB AB AA AB \n", "4 AA AB BB BB BB BB AB BB AB BB ... AA AA BB BB AB AB AA AB \n", "5 AA AB AA BB AA AA AA BB AA BB ... BB AB AA BB BB AB AA BB \n", " V999 V1000\n", "1 AB AB \n", "2 BB AA \n", "3 AB AB \n", "4 BB AA \n", "5 BB AA " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Genotypes2" ] }, { "cell_type": "code", "execution_count": 238, "metadata": {}, "outputs": [], "source": [ "p2 = ifelse(p2 > 0.5, 1-p2, p2)" ] }, { "cell_type": "code", "execution_count": 239, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " Min. 1st Qu. Median Mean 3rd Qu. Max. \n", " 0.1000 0.3000 0.4000 0.3535 0.4000 0.5000 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "summary(p2)" ] }, { "cell_type": "code", "execution_count": 240, "metadata": {}, "outputs": [], "source": [ "Genotypes2_2 = matrix(0, 5, ncol(Genotypes2))\n", "for (i in 1:ncol(Genotypes2)) {\n", " Genotypes2_2[Genotypes2[,i] == \"AB\", i] = 1\n", " Genotypes2_2[Genotypes2[,i] == \"BB\", i] = 2\n", "}" ] }, { "cell_type": "code", "execution_count": 241, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
2 0 2 0 2 0 0 1 1 2 ...2 1 2 0 0 1 1 1 1 1
1 0 2 2 1 2 0 1 0 2 ...0 2 2 0 0 0 2 2 2 0
2 0 0 0 2 0 2 0 2 0 ...1 0 0 2 2 1 0 1 1 1
0 1 2 2 2 2 1 2 1 2 ...0 0 2 2 1 1 0 1 2 0
0 1 0 2 0 0 0 2 0 2 ...2 1 0 2 2 1 0 2 2 0
\n" ], "text/latex": [ "\\begin{tabular}{llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll}\n", "\t 2 & 0 & 2 & 0 & 2 & 0 & 0 & 1 & 1 & 2 & ... & 2 & 1 & 2 & 0 & 0 & 1 & 1 & 1 & 1 & 1 \\\\\n", "\t 1 & 0 & 2 & 2 & 1 & 2 & 0 & 1 & 0 & 2 & ... & 0 & 2 & 2 & 0 & 0 & 0 & 2 & 2 & 2 & 0 \\\\\n", "\t 2 & 0 & 0 & 0 & 2 & 0 & 2 & 0 & 2 & 0 & ... & 1 & 0 & 0 & 2 & 2 & 1 & 0 & 1 & 1 & 1 \\\\\n", "\t 0 & 1 & 2 & 2 & 2 & 2 & 1 & 2 & 1 & 2 & ... & 0 & 0 & 2 & 2 & 1 & 1 & 0 & 1 & 2 & 0 \\\\\n", "\t 0 & 1 & 0 & 2 & 0 & 0 & 0 & 2 & 0 & 2 & ... & 2 & 1 & 0 & 2 & 2 & 1 & 0 & 2 & 2 & 0 \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 2 | 0 | 2 | 0 | 2 | 0 | 0 | 1 | 1 | 2 | ... | 2 | 1 | 2 | 0 | 0 | 1 | 1 | 1 | 1 | 1 |\n", "| 1 | 0 | 2 | 2 | 1 | 2 | 0 | 1 | 0 | 2 | ... | 0 | 2 | 2 | 0 | 0 | 0 | 2 | 2 | 2 | 0 |\n", "| 2 | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 2 | 0 | ... | 1 | 0 | 0 | 2 | 2 | 1 | 0 | 1 | 1 | 1 |\n", "| 0 | 1 | 2 | 2 | 2 | 2 | 1 | 2 | 1 | 2 | ... | 0 | 0 | 2 | 2 | 1 | 1 | 0 | 1 | 2 | 0 |\n", "| 0 | 1 | 0 | 2 | 0 | 0 | 0 | 2 | 0 | 2 | ... | 2 | 1 | 0 | 2 | 2 | 1 | 0 | 2 | 2 | 0 |\n", "\n" ], "text/plain": [ " [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]\n", "[1,] 2 0 2 0 2 0 0 1 1 2 ... 2 1 2 \n", "[2,] 1 0 2 2 1 2 0 1 0 2 ... 0 2 2 \n", "[3,] 2 0 0 0 2 0 2 0 2 0 ... 1 0 0 \n", "[4,] 0 1 2 2 2 2 1 2 1 2 ... 0 0 2 \n", "[5,] 0 1 0 2 0 0 0 2 0 2 ... 2 1 0 \n", " [,15] [,16] [,17] [,18] [,19] [,20] [,21]\n", "[1,] 0 0 1 1 1 1 1 \n", "[2,] 0 0 0 2 2 2 0 \n", "[3,] 2 2 1 0 1 1 1 \n", "[4,] 2 1 1 0 1 2 0 \n", "[5,] 2 2 1 0 2 2 0 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Genotypes2_2" ] }, { "cell_type": "code", "execution_count": 242, "metadata": {}, "outputs": [], "source": [ "M = matrix(0, 5, length(p2))\n", "for (i in 1:5) {\n", " for (j in 1:length(p2)) {\n", " if (Genotypes2_2[i, j] == 0) {\n", " M[i, j] = - 2*p2[j]\n", " } else if (Genotypes2_2[i, j] == 1) {\n", " M[i, j] = 1 - 2*p2[j]\n", " } else {\n", " M[i, j] = 2 - 2*p2[j]\n", " }\n", " }\n", "}" ] }, { "cell_type": "code", "execution_count": 243, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
1 -0.4 1.2-0.8 1.4-0.8-0.6 0.2 0.2 1.6... 1 0.4 1.4-0.6-0.4 0.8 0.80 0 0
0 -0.4 1.2 1.2 0.4 1.2-0.6 0.2-0.8 1.6... -1 1.4 1.4-0.6-0.4-0.2 1.81 1 -1
1 -0.4-0.8-0.8 1.4-0.8 1.4-0.8 1.2-0.4... 0 -0.6-0.6 1.4 1.6 0.8-0.20 0 0
-1 0.6 1.2 1.2 1.4 1.2 0.4 1.2 0.2 1.6... -1 -0.6 1.4 1.4 0.6 0.8-0.20 1 -1
-1 0.6-0.8 1.2-0.6-0.8-0.6 1.2-0.8 1.6... 1 0.4-0.6 1.4 1.6 0.8-0.21 1 -1
\n" ], "text/latex": [ "\\begin{tabular}{llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll}\n", "\t 1 & -0.4 & 1.2 & -0.8 & 1.4 & -0.8 & -0.6 & 0.2 & 0.2 & 1.6 & ... & 1 & 0.4 & 1.4 & -0.6 & -0.4 & 0.8 & 0.8 & 0 & 0 & 0 \\\\\n", "\t 0 & -0.4 & 1.2 & 1.2 & 0.4 & 1.2 & -0.6 & 0.2 & -0.8 & 1.6 & ... & -1 & 1.4 & 1.4 & -0.6 & -0.4 & -0.2 & 1.8 & 1 & 1 & -1 \\\\\n", "\t 1 & -0.4 & -0.8 & -0.8 & 1.4 & -0.8 & 1.4 & -0.8 & 1.2 & -0.4 & ... & 0 & -0.6 & -0.6 & 1.4 & 1.6 & 0.8 & -0.2 & 0 & 0 & 0 \\\\\n", "\t -1 & 0.6 & 1.2 & 1.2 & 1.4 & 1.2 & 0.4 & 1.2 & 0.2 & 1.6 & ... & -1 & -0.6 & 1.4 & 1.4 & 0.6 & 0.8 & -0.2 & 0 & 1 & -1 \\\\\n", "\t -1 & 0.6 & -0.8 & 1.2 & -0.6 & -0.8 & -0.6 & 1.2 & -0.8 & 1.6 & ... & 1 & 0.4 & -0.6 & 1.4 & 1.6 & 0.8 & -0.2 & 1 & 1 & -1 \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 1 | -0.4 | 1.2 | -0.8 | 1.4 | -0.8 | -0.6 | 0.2 | 0.2 | 1.6 | ... | 1 | 0.4 | 1.4 | -0.6 | -0.4 | 0.8 | 0.8 | 0 | 0 | 0 |\n", "| 0 | -0.4 | 1.2 | 1.2 | 0.4 | 1.2 | -0.6 | 0.2 | -0.8 | 1.6 | ... | -1 | 1.4 | 1.4 | -0.6 | -0.4 | -0.2 | 1.8 | 1 | 1 | -1 |\n", "| 1 | -0.4 | -0.8 | -0.8 | 1.4 | -0.8 | 1.4 | -0.8 | 1.2 | -0.4 | ... | 0 | -0.6 | -0.6 | 1.4 | 1.6 | 0.8 | -0.2 | 0 | 0 | 0 |\n", "| -1 | 0.6 | 1.2 | 1.2 | 1.4 | 1.2 | 0.4 | 1.2 | 0.2 | 1.6 | ... | -1 | -0.6 | 1.4 | 1.4 | 0.6 | 0.8 | -0.2 | 0 | 1 | -1 |\n", "| -1 | 0.6 | -0.8 | 1.2 | -0.6 | -0.8 | -0.6 | 1.2 | -0.8 | 1.6 | ... | 1 | 0.4 | -0.6 | 1.4 | 1.6 | 0.8 | -0.2 | 1 | 1 | -1 |\n", "\n" ], "text/plain": [ " [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]\n", "[1,] 1 -0.4 1.2 -0.8 1.4 -0.8 -0.6 0.2 0.2 1.6 ... 1 0.4 1.4 \n", "[2,] 0 -0.4 1.2 1.2 0.4 1.2 -0.6 0.2 -0.8 1.6 ... -1 1.4 1.4 \n", "[3,] 1 -0.4 -0.8 -0.8 1.4 -0.8 1.4 -0.8 1.2 -0.4 ... 0 -0.6 -0.6 \n", "[4,] -1 0.6 1.2 1.2 1.4 1.2 0.4 1.2 0.2 1.6 ... -1 -0.6 1.4 \n", "[5,] -1 0.6 -0.8 1.2 -0.6 -0.8 -0.6 1.2 -0.8 1.6 ... 1 0.4 -0.6 \n", " [,15] [,16] [,17] [,18] [,19] [,20] [,21]\n", "[1,] -0.6 -0.4 0.8 0.8 0 0 0 \n", "[2,] -0.6 -0.4 -0.2 1.8 1 1 -1 \n", "[3,] 1.4 1.6 0.8 -0.2 0 0 0 \n", "[4,] 1.4 0.6 0.8 -0.2 0 1 -1 \n", "[5,] 1.4 1.6 0.8 -0.2 1 1 -1 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "M" ] }, { "cell_type": "code", "execution_count": 246, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Initial data: \n", "\tNumber of Individuals: 5 \n", "\tNumber of Markers: 984 \n", "\n", "Missing data check: \n", "\tTotal SNPs: 984 \n", "\t 0 SNPs dropped due to missing data threshold of 1 \n", "\tTotal of: 984 SNPs \n", "MAF check: \n", "\tNo SNPs with MAF below 0 \n", "Monomorphic check: \n", "\tNo monomorphic SNPs \n", "Summary check: \n", "\tInitial: 984 SNPs \n", "\tFinal: 984 SNPs ( 0 SNPs removed) \n", " \n", "Completed! Time = 0.04 seconds \n" ] }, { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
1.2288695-0.2286356-0.3235886-0.3399598-0.3366855
-0.2286356 1.2536601-0.2702652-0.3591375-0.3956219
-0.3235886-0.2702652 1.2227887-0.2927172-0.3362178
-0.3399598-0.3591375-0.2927172 1.2017400-0.2099256
-0.3366855-0.3956219-0.3362178-0.2099256 1.2784508
\n" ], "text/latex": [ "\\begin{tabular}{lllll}\n", "\t 1.2288695 & -0.2286356 & -0.3235886 & -0.3399598 & -0.3366855\\\\\n", "\t -0.2286356 & 1.2536601 & -0.2702652 & -0.3591375 & -0.3956219\\\\\n", "\t -0.3235886 & -0.2702652 & 1.2227887 & -0.2927172 & -0.3362178\\\\\n", "\t -0.3399598 & -0.3591375 & -0.2927172 & 1.2017400 & -0.2099256\\\\\n", "\t -0.3366855 & -0.3956219 & -0.3362178 & -0.2099256 & 1.2784508\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 1.2288695 | -0.2286356 | -0.3235886 | -0.3399598 | -0.3366855 |\n", "| -0.2286356 | 1.2536601 | -0.2702652 | -0.3591375 | -0.3956219 |\n", "| -0.3235886 | -0.2702652 | 1.2227887 | -0.2927172 | -0.3362178 |\n", "| -0.3399598 | -0.3591375 | -0.2927172 | 1.2017400 | -0.2099256 |\n", "| -0.3366855 | -0.3956219 | -0.3362178 | -0.2099256 | 1.2784508 |\n", "\n" ], "text/plain": [ " [,1] [,2] [,3] [,4] [,5] \n", "[1,] 1.2288695 -0.2286356 -0.3235886 -0.3399598 -0.3366855\n", "[2,] -0.2286356 1.2536601 -0.2702652 -0.3591375 -0.3956219\n", "[3,] -0.3235886 -0.2702652 1.2227887 -0.2927172 -0.3362178\n", "[4,] -0.3399598 -0.3591375 -0.2927172 1.2017400 -0.2099256\n", "[5,] -0.3366855 -0.3956219 -0.3362178 -0.2099256 1.2784508" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Gmatrix(SNPmatrix = Genotypes2_2,\n", " method = \"VanRaden\",\n", " missingValue = -9,\n", " maf = 0.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "12. Wrzucenie macierzy G do modelu" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathbf{y} = \\mathbf{X\\beta+Z_{1}a+Z_{2}g+\\epsilon}$" ] }, { "cell_type": "code", "execution_count": 32, "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.341304960
3.401214252
-0.003452949
-0.055133673
0.054219329
-0.072283409
0.054390036
0.087512029
-0.246672434
0.245158275
-0.425100763
0.335020177
\n" ], "text/latex": [ "\\begin{tabular}{l}\n", "\t 4.341304960\\\\\n", "\t 3.401214252\\\\\n", "\t -0.003452949\\\\\n", "\t -0.055133673\\\\\n", "\t 0.054219329\\\\\n", "\t -0.072283409\\\\\n", "\t 0.054390036\\\\\n", "\t 0.087512029\\\\\n", "\t -0.246672434\\\\\n", "\t 0.245158275\\\\\n", "\t -0.425100763\\\\\n", "\t 0.335020177\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 4.341304960 |\n", "| 3.401214252 |\n", "| -0.003452949 |\n", "| -0.055133673 |\n", "| 0.054219329 |\n", "| -0.072283409 |\n", "| 0.054390036 |\n", "| 0.087512029 |\n", "| -0.246672434 |\n", "| 0.245158275 |\n", "| -0.425100763 |\n", "| 0.335020177 |\n", "\n" ], "text/plain": [ " [,1] \n", " [1,] 4.341304960\n", " [2,] 3.401214252\n", " [3,] -0.003452949\n", " [4,] -0.055133673\n", " [5,] 0.054219329\n", " [6,] -0.072283409\n", " [7,] 0.054390036\n", " [8,] 0.087512029\n", " [9,] -0.246672434\n", "[10,] 0.245158275\n", "[11,] -0.425100763\n", "[12,] 0.335020177" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mme2 = function(y, X, Z1, Z2, A, G, sigma_a, sigma_g, sigma_e) {\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 = mme2(y, X, diag(5), diag(5), A[4:8, 4:8], G, 10, 30, 40)$est)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "13. Efekty markeró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": 33, "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": 34, "metadata": {}, "outputs": [], "source": [ "snp = snp%*%t(M)%*%ginv(A)[4:8, 4:8]" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [], "source": [ "snp = snp %*% est[8:12]" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
-0.0013010329
-0.0006039244
0.0011944542
-0.0003928861
-0.0002694417
-0.0002621532
\n" ], "text/latex": [ "\\begin{tabular}{l}\n", "\t -0.0013010329\\\\\n", "\t -0.0006039244\\\\\n", "\t 0.0011944542\\\\\n", "\t -0.0003928861\\\\\n", "\t -0.0002694417\\\\\n", "\t -0.0002621532\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| -0.0013010329 |\n", "| -0.0006039244 |\n", "| 0.0011944542 |\n", "| -0.0003928861 |\n", "| -0.0002694417 |\n", "| -0.0002621532 |\n", "\n" ], "text/plain": [ " [,1] \n", "[1,] -0.0013010329\n", "[2,] -0.0006039244\n", "[3,] 0.0011944542\n", "[4,] -0.0003928861\n", "[5,] -0.0002694417\n", "[6,] -0.0002621532" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "head(snp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "14. Czy efekty SNP są statystycznie istotne?" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "0.000776102778820242" ], "text/latex": [ "0.000776102778820242" ], "text/markdown": [ "0.000776102778820242" ], "text/plain": [ "[1] 0.0007761028" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "-5.11112732482814e-20" ], "text/latex": [ "-5.11112732482814e-20" ], "text/markdown": [ "-5.11112732482814e-20" ], "text/plain": [ "[1] -5.111127e-20" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "(s = sd(snp))\n", "(m= mean(snp))" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "1" ], "text/latex": [ "1" ], "text/markdown": [ "1" ], "text/plain": [ "[1] 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "snp2 = snp / s\n", "sd(snp2)" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "W = snp2\n", "p.value = 2*pnorm(abs(W), lower.tail = FALSE)" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "40" ], "text/latex": [ "40" ], "text/markdown": [ "40" ], "text/plain": [ "[1] 40" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "length(which(p.value < 0.05))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "15. Estymacja efektów markerów - metoda GBLUP" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathbf{y} = \\mathbf{X\\beta+Z_{1}a+Z_{2}g+\\epsilon}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$Z_{2} = \\{-1, 0, 1\\}$ odpowiednio dla genotypów $\\{AA, AB, BB\\}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Efekt $\\mathbf{g}\\sim\\mathcal{N}(0, \\mathbf{I}\\cdot\\frac{\\sigma^{2}_{a}}{N_{SNP}})$" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [], "source": [ "Z2 = matrix(0, 5, length(p2))\n", "for (i in 1:5) {\n", " for (j in 1:length(p2)) {\n", " if (Genotypes2[i, j] == \"AA\") {\n", " Z2[i, j] = -1\n", " } else if (Genotypes2[i, j] == \"AB\") {\n", " Z2[i, j] = 0\n", " } else {\n", " Z2[i, j] = 1\n", " }\n", " }\n", "}" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
-1 -1 -1 1 1 1 0 1 1 0 ... 0 0 -1 -1 1 1 0 1 0 1
-1 1 1 1 -1 1 0 1 1 1 ... 1 0 0 1 -1 1 1 0 1 0
0 0 1 0 -1 0 -1 1 0 1 ...-1 -1 1 -1 -1 1 1 1 1 1
-1 -1 1 -1 -1 0 -1 0 0 1 ... 0 1 1 -1 0 0 -1 -1 0 0
1 1 -1 0 -1 1 -1 0 -1 0 ... 1 1 0 0 1 0 1 0 -1 -1
\n" ], "text/latex": [ "\\begin{tabular}{llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll}\n", "\t -1 & -1 & -1 & 1 & 1 & 1 & 0 & 1 & 1 & 0 & ... & 0 & 0 & -1 & -1 & 1 & 1 & 0 & 1 & 0 & 1 \\\\\n", "\t -1 & 1 & 1 & 1 & -1 & 1 & 0 & 1 & 1 & 1 & ... & 1 & 0 & 0 & 1 & -1 & 1 & 1 & 0 & 1 & 0 \\\\\n", "\t 0 & 0 & 1 & 0 & -1 & 0 & -1 & 1 & 0 & 1 & ... & -1 & -1 & 1 & -1 & -1 & 1 & 1 & 1 & 1 & 1 \\\\\n", "\t -1 & -1 & 1 & -1 & -1 & 0 & -1 & 0 & 0 & 1 & ... & 0 & 1 & 1 & -1 & 0 & 0 & -1 & -1 & 0 & 0 \\\\\n", "\t 1 & 1 & -1 & 0 & -1 & 1 & -1 & 0 & -1 & 0 & ... & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & -1 & -1 \\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| -1 | -1 | -1 | 1 | 1 | 1 | 0 | 1 | 1 | 0 | ... | 0 | 0 | -1 | -1 | 1 | 1 | 0 | 1 | 0 | 1 |\n", "| -1 | 1 | 1 | 1 | -1 | 1 | 0 | 1 | 1 | 1 | ... | 1 | 0 | 0 | 1 | -1 | 1 | 1 | 0 | 1 | 0 |\n", "| 0 | 0 | 1 | 0 | -1 | 0 | -1 | 1 | 0 | 1 | ... | -1 | -1 | 1 | -1 | -1 | 1 | 1 | 1 | 1 | 1 |\n", "| -1 | -1 | 1 | -1 | -1 | 0 | -1 | 0 | 0 | 1 | ... | 0 | 1 | 1 | -1 | 0 | 0 | -1 | -1 | 0 | 0 |\n", "| 1 | 1 | -1 | 0 | -1 | 1 | -1 | 0 | -1 | 0 | ... | 1 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | -1 | -1 |\n", "\n" ], "text/plain": [ " [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]\n", "[1,] -1 -1 -1 1 1 1 0 1 1 0 ... 0 0 -1 \n", "[2,] -1 1 1 1 -1 1 0 1 1 1 ... 1 0 0 \n", "[3,] 0 0 1 0 -1 0 -1 1 0 1 ... -1 -1 1 \n", "[4,] -1 -1 1 -1 -1 0 -1 0 0 1 ... 0 1 1 \n", "[5,] 1 1 -1 0 -1 1 -1 0 -1 0 ... 1 1 0 \n", " [,15] [,16] [,17] [,18] [,19] [,20] [,21]\n", "[1,] -1 1 1 0 1 0 1 \n", "[2,] 1 -1 1 1 0 1 0 \n", "[3,] -1 -1 1 1 1 1 1 \n", "[4,] -1 0 0 -1 -1 0 0 \n", "[5,] 0 1 0 1 0 -1 -1 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "head(Z2)" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [], "source": [ "sigma_a = 0.6555448 \n", "sigma_e = 0.01799635\n", "sigma_g = sigma_a / length(p2)\n", "G = diag(length(p2))" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [], "source": [ "res2 = mme2(y, X, diag(5), Z2, A[4:8, 4:8], G, sigma_a, sigma_g, sigma_e)" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 999
  2. \n", "\t
  3. 999
  4. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 999\n", "\\item 999\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 999\n", "2. 999\n", "\n", "\n" ], "text/plain": [ "[1] 999 999" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "dim(res2$C)" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 0.00107632007119809
  2. \n", "\t
  3. 0.00056523214301077
  4. \n", "\t
  5. -0.00117104593778478
  6. \n", "\t
  7. 0.000505113920098791
  8. \n", "\t
  9. 0.000350269830614274
  10. \n", "\t
  11. 0.000329979004801596
  12. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.00107632007119809\n", "\\item 0.00056523214301077\n", "\\item -0.00117104593778478\n", "\\item 0.000505113920098791\n", "\\item 0.000350269830614274\n", "\\item 0.000329979004801596\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.00107632007119809\n", "2. 0.00056523214301077\n", "3. -0.00117104593778478\n", "4. 0.000505113920098791\n", "5. 0.000350269830614274\n", "6. 0.000329979004801596\n", "\n", "\n" ], "text/plain": [ "[1] 0.0010763201 0.0005652321 -0.0011710459 0.0005051139 0.0003502698\n", "[6] 0.0003299790" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "est_snp2 = res2$est[8:999]\n", "head(est_snp2)" ] }, { "cell_type": "code", "execution_count": 57, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "0.000691430096493191" ], "text/latex": [ "0.000691430096493191" ], "text/markdown": [ "0.000691430096493191" ], "text/plain": [ "[1] 0.0006914301" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "-8.63480392554709e-06" ], "text/latex": [ "-8.63480392554709e-06" ], "text/markdown": [ "-8.63480392554709e-06" ], "text/plain": [ "[1] -8.634804e-06" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "(s = sd(est_snp2))\n", "(m= mean(est_snp2))" ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "1" ], "text/latex": [ "1" ], "text/markdown": [ "1" ], "text/plain": [ "[1] 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "est_snp2 = est_snp2 / s\n", "sd(est_snp2)" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [], "source": [ "W2 = est_snp2\n", "p.value2 = 2*pnorm(abs(W2), lower.tail = FALSE)" ] }, { "cell_type": "code", "execution_count": 61, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "34" ], "text/latex": [ "34" ], "text/markdown": [ "34" ], "text/plain": [ "[1] 34" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "length(which(p.value2 < 0.05))" ] }, { "cell_type": "code", "execution_count": 62, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 38
  2. \n", "\t
  3. 42
  4. \n", "\t
  5. 69
  6. \n", "\t
  7. 75
  8. \n", "\t
  9. 93
  10. \n", "\t
  11. 107
  12. \n", "\t
  13. 136
  14. \n", "\t
  15. 138
  16. \n", "\t
  17. 145
  18. \n", "\t
  19. 204
  20. \n", "\t
  21. 223
  22. \n", "\t
  23. 230
  24. \n", "\t
  25. 257
  26. \n", "\t
  27. 268
  28. \n", "\t
  29. 302
  30. \n", "\t
  31. 312
  32. \n", "\t
  33. 346
  34. \n", "\t
  35. 349
  36. \n", "\t
  37. 376
  38. \n", "\t
  39. 434
  40. \n", "\t
  41. 436
  42. \n", "\t
  43. 452
  44. \n", "\t
  45. 462
  46. \n", "\t
  47. 463
  48. \n", "\t
  49. 507
  50. \n", "\t
  51. 543
  52. \n", "\t
  53. 561
  54. \n", "\t
  55. 576
  56. \n", "\t
  57. 604
  58. \n", "\t
  59. 626
  60. \n", "\t
  61. 695
  62. \n", "\t
  63. 715
  64. \n", "\t
  65. 733
  66. \n", "\t
  67. 739
  68. \n", "\t
  69. 768
  70. \n", "\t
  71. 787
  72. \n", "\t
  73. 791
  74. \n", "\t
  75. 799
  76. \n", "\t
  77. 811
  78. \n", "\t
  79. 972
  80. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 38\n", "\\item 42\n", "\\item 69\n", "\\item 75\n", "\\item 93\n", "\\item 107\n", "\\item 136\n", "\\item 138\n", "\\item 145\n", "\\item 204\n", "\\item 223\n", "\\item 230\n", "\\item 257\n", "\\item 268\n", "\\item 302\n", "\\item 312\n", "\\item 346\n", "\\item 349\n", "\\item 376\n", "\\item 434\n", "\\item 436\n", "\\item 452\n", "\\item 462\n", "\\item 463\n", "\\item 507\n", "\\item 543\n", "\\item 561\n", "\\item 576\n", "\\item 604\n", "\\item 626\n", "\\item 695\n", "\\item 715\n", "\\item 733\n", "\\item 739\n", "\\item 768\n", "\\item 787\n", "\\item 791\n", "\\item 799\n", "\\item 811\n", "\\item 972\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 38\n", "2. 42\n", "3. 69\n", "4. 75\n", "5. 93\n", "6. 107\n", "7. 136\n", "8. 138\n", "9. 145\n", "10. 204\n", "11. 223\n", "12. 230\n", "13. 257\n", "14. 268\n", "15. 302\n", "16. 312\n", "17. 346\n", "18. 349\n", "19. 376\n", "20. 434\n", "21. 436\n", "22. 452\n", "23. 462\n", "24. 463\n", "25. 507\n", "26. 543\n", "27. 561\n", "28. 576\n", "29. 604\n", "30. 626\n", "31. 695\n", "32. 715\n", "33. 733\n", "34. 739\n", "35. 768\n", "36. 787\n", "37. 791\n", "38. 799\n", "39. 811\n", "40. 972\n", "\n", "\n" ], "text/plain": [ " [1] 38 42 69 75 93 107 136 138 145 204 223 230 257 268 302 312 346 349 376\n", "[20] 434 436 452 462 463 507 543 561 576 604 626 695 715 733 739 768 787 791 799\n", "[39] 811 972" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 38
  2. \n", "\t
  3. 42
  4. \n", "\t
  5. 69
  6. \n", "\t
  7. 75
  8. \n", "\t
  9. 93
  10. \n", "\t
  11. 107
  12. \n", "\t
  13. 136
  14. \n", "\t
  15. 138
  16. \n", "\t
  17. 145
  18. \n", "\t
  19. 204
  20. \n", "\t
  21. 223
  22. \n", "\t
  23. 230
  24. \n", "\t
  25. 257
  26. \n", "\t
  27. 268
  28. \n", "\t
  29. 312
  30. \n", "\t
  31. 346
  32. \n", "\t
  33. 349
  34. \n", "\t
  35. 376
  36. \n", "\t
  37. 436
  38. \n", "\t
  39. 452
  40. \n", "\t
  41. 462
  42. \n", "\t
  43. 463
  44. \n", "\t
  45. 507
  46. \n", "\t
  47. 561
  48. \n", "\t
  49. 576
  50. \n", "\t
  51. 604
  52. \n", "\t
  53. 626
  54. \n", "\t
  55. 695
  56. \n", "\t
  57. 739
  58. \n", "\t
  59. 768
  60. \n", "\t
  61. 787
  62. \n", "\t
  63. 791
  64. \n", "\t
  65. 811
  66. \n", "\t
  67. 972
  68. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 38\n", "\\item 42\n", "\\item 69\n", "\\item 75\n", "\\item 93\n", "\\item 107\n", "\\item 136\n", "\\item 138\n", "\\item 145\n", "\\item 204\n", "\\item 223\n", "\\item 230\n", "\\item 257\n", "\\item 268\n", "\\item 312\n", "\\item 346\n", "\\item 349\n", "\\item 376\n", "\\item 436\n", "\\item 452\n", "\\item 462\n", "\\item 463\n", "\\item 507\n", "\\item 561\n", "\\item 576\n", "\\item 604\n", "\\item 626\n", "\\item 695\n", "\\item 739\n", "\\item 768\n", "\\item 787\n", "\\item 791\n", "\\item 811\n", "\\item 972\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 38\n", "2. 42\n", "3. 69\n", "4. 75\n", "5. 93\n", "6. 107\n", "7. 136\n", "8. 138\n", "9. 145\n", "10. 204\n", "11. 223\n", "12. 230\n", "13. 257\n", "14. 268\n", "15. 312\n", "16. 346\n", "17. 349\n", "18. 376\n", "19. 436\n", "20. 452\n", "21. 462\n", "22. 463\n", "23. 507\n", "24. 561\n", "25. 576\n", "26. 604\n", "27. 626\n", "28. 695\n", "29. 739\n", "30. 768\n", "31. 787\n", "32. 791\n", "33. 811\n", "34. 972\n", "\n", "\n" ], "text/plain": [ " [1] 38 42 69 75 93 107 136 138 145 204 223 230 257 268 312 346 349 376 436\n", "[20] 452 462 463 507 561 576 604 626 695 739 768 787 791 811 972" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "which(p.value < 0.05)\n", "which(p.value2 < 0.05)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "16. Model wielocechowy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Mamy $\\textbf{y}_{1} = \\textbf{X}_{1}\\beta_{1}+\\textbf{Z}_{1}\\textbf{a}_{1}+\\epsilon_{1}$ oraz $\\textbf{y}_{2} = \\textbf{X}_{2}\\beta_{2}+\\textbf{Z}_{2}\\textbf{a}_{2}+\\epsilon_{2}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$ \\left[ \\begin{array}{c}\n", " \\mathbf{y}_{1} \\\\\n", " \\mathbf{y}_{2}\\end{array}\\right] = \n", " \\left[ \\begin{array}{cc}\n", " \\mathbf{X}_{1} & \\mathbf{0} \\\\\n", " \\mathbf{0} & \\mathbf{X}_{2} \\end{array}\\right]\\cdot\n", " \\left[ \\begin{array}{c}\n", " \\mathbf{\\beta}_{1} \\\\\n", " \\mathbf{\\beta}_{2} \\end{array}\\right] +\n", " \\left[ \\begin{array}{cc}\n", " \\mathbf{Z}_{1} & \\mathbf{0} \\\\\n", " \\mathbf{0} & \\mathbf{Z}_{2} \\end{array}\\right]\\cdot\n", " \\left[ \\begin{array}{c}\n", " \\mathbf{a}_{1} \\\\\n", " \\mathbf{a}_{2} \\end{array}\\right] +\n", " \\left[ \\begin{array}{c}\n", " \\mathbf{\\epsilon}_{1} \\\\\n", " \\mathbf{\\epsilon}_{2}\\end{array}\\right]$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Korelacja pomiedzy efektami" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$var\\left[ \\begin{array}{c}\n", " \\mathbf{a}_{1} \\\\\n", " \\mathbf{a}_{2} \\\\ \n", " \\mathbf{\\epsilon}_{1} \\\\\n", " \\mathbf{\\epsilon}_{2}\n", " \\end{array}\\right] = \n", " \\left[ \\begin{array}{cccc}\n", " \\mathbf{g}_{11}\\mathbf{A} & \\mathbf{g}_{12}\\mathbf{A} & \\mathbf{0} & \\mathbf{0} \\\\\n", " \\mathbf{g}_{21}\\mathbf{A} & \\mathbf{g}_{22}\\mathbf{A} & \\mathbf{0} & \\mathbf{0} \\\\ \n", " \\mathbf{0} & \\mathbf{0} & \\mathbf{r}_{11}\\mathbf{I} & \\mathbf{r}_{12}\\mathbf{I} \\\\\n", " \\mathbf{0} & \\mathbf{0} & \\mathbf{r}_{21}\\mathbf{I} & \\mathbf{r}_{22}\\mathbf{I}\n", " \\end{array}\\right] $" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Układ równań mieszanych" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$ \\left[ \\begin{array}{cc}\n", " \\mathbf{X}^{T}\\mathbf{R}^{-1}\\mathbf{X} & \\mathbf{X}^{T}\\mathbf{R}^{-1}\\mathbf{Z} \\\\\n", " \\mathbf{Z}^{T}\\mathbf{R}^{-1}\\mathbf{X} & \\mathbf{Z}^{T}\\mathbf{R}^{-1}\\mathbf{Z}+\\mathbf{A}^{-1}\\otimes \\mathbf{G}^{-1} \\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{R}^{-1}\\mathbf{y} \\\\\n", " \\mathbf{Z}^{T}\\mathbf{R}^{-1}\\mathbf{y} \\end{array}\\right]$, gdzie" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$ \\mathbf{y} = \\left[ \\begin{array}{c}\n", " \\mathbf{y}_{1} \\\\\n", " \\mathbf{y}_{2}\\end{array}\\right]$, $ \\mathbf{X} = \\left[ \\begin{array}{cc}\n", " \\mathbf{X}_{1} & \\mathbf{0} \\\\\n", " \\mathbf{0} & \\mathbf{X}_{2} \\end{array}\\right]$, \n", "$ \\mathbf{\\beta} = \\left[ \\begin{array}{c}\n", " \\mathbf{\\beta}_{1} \\\\\n", " \\mathbf{\\beta}_{2} \\end{array}\\right]$, \n", "$ \\mathbf{Z} = \\left[ \\begin{array}{cc}\n", " \\mathbf{Z}_{1} & \\mathbf{0} \\\\\n", " \\mathbf{0} & \\mathbf{Z}_{2} \\end{array}\\right]$,\n", "$ \\mathbf{a} = \\left[ \\begin{array}{c}\n", " \\mathbf{a}_{1} \\\\\n", " \\mathbf{a}_{2} \\end{array}\\right]$,\n", "$ \\mathbf{\\epsilon} = \\left[ \\begin{array}{c}\n", " \\mathbf{\\epsilon}_{1} \\\\\n", " \\mathbf{\\epsilon}_{2}\\end{array}\\right]$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Macierz A" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": true }, "outputs": [], "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": 4, "metadata": { "scrolled": true }, "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" ] }, { "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(Amatrix(ped)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pozostałe źródła informacji" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\n", "
2018
1840
\n" ], "text/latex": [ "\\begin{tabular}{ll}\n", "\t 20 & 18\\\\\n", "\t 18 & 40\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 20 | 18 |\n", "| 18 | 40 |\n", "\n" ], "text/plain": [ " [,1] [,2]\n", "[1,] 20 18 \n", "[2,] 18 40 " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\n", "
4011
1130
\n" ], "text/latex": [ "\\begin{tabular}{ll}\n", "\t 40 & 11\\\\\n", "\t 11 & 30\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 40 | 11 |\n", "| 11 | 30 |\n", "\n" ], "text/plain": [ " [,1] [,2]\n", "[1,] 40 11 \n", "[2,] 11 30 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "y1 = as.matrix(c(4.5, 2.9, 3.9, 3.5, 5.0))\n", "y2 = as.matrix(c(6.8, 5.0, 6.8, 6.0, 7.5))\n", "\n", "(G = matrix(c(20, 18, 18, 40), 2, 2))\n", "(R = matrix(c(40, 11, 11, 30), 2, 2))" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/html": [ "0.333333333333333" ], "text/latex": [ "0.333333333333333" ], "text/markdown": [ "0.333333333333333" ], "text/plain": [ "[1] 0.3333333" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "0.571428571428571" ], "text/latex": [ "0.571428571428571" ], "text/markdown": [ "0.571428571428571" ], "text/plain": [ "[1] 0.5714286" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "20/60\n", "40/70" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "sex = c(1, 0, 0, 1, 1)\n", "X1 = matrix(0, 5, 2)\n", "X1[, 1] = sex\n", "X1[, 2] = 1-sex\n", "X2 = X1" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "I = diag(5)\n", "Z1 = matrix(0, 5, 8)\n", "Z1[1:5, 4:8] = I\n", "Z2 = Z1" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "library(MASS)\n", "\n", "mme3 = function(y, X, Z, A, G, R) {\n", " invA = ginv(A)\n", " invG = ginv(G)\n", " R = kronecker(R, diag(5))\n", " invR = ginv(R)\n", " C = rbind(cbind(t(X)%*%invR%*%X, t(X)%*%invR%*%Z),\n", " cbind(t(Z)%*%invR%*%X, t(Z)%*%invR%*%Z+kronecker(invG, invA)))\n", " rhs = rbind(t(X)%*%invR%*%y, t(Z)%*%invR%*%y)\n", " invC = ginv(C)\n", " estimators = invC%*%rhs\n", " list(C = C, est = estimators)\n", "}" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "scrolled": true }, "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", "\n", "
0001000000000000
0000100000000000
0000010000000000
0000001000000000
0000000100000000
0000000000010000
0000000000001000
0000000000000100
0000000000000010
0000000000000001
\n" ], "text/latex": [ "\\begin{tabular}{llllllllllllllll}\n", "\t 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\\n", "\t 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\\n", "\t 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\\n", "\t 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\\n", "\t 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\\n", "\t 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\\\\n", "\t 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\\\\n", "\t 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\\\\n", "\t 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\\\\n", "\t 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |\n", "| 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |\n", "| 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |\n", "| 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |\n", "| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |\n", "| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 |\n", "| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |\n", "| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 |\n", "| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 |\n", "| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 |\n", "\n" ], "text/plain": [ " [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]\n", " [1,] 0 0 0 1 0 0 0 0 0 0 0 0 0 \n", " [2,] 0 0 0 0 1 0 0 0 0 0 0 0 0 \n", " [3,] 0 0 0 0 0 1 0 0 0 0 0 0 0 \n", " [4,] 0 0 0 0 0 0 1 0 0 0 0 0 0 \n", " [5,] 0 0 0 0 0 0 0 1 0 0 0 0 0 \n", " [6,] 0 0 0 0 0 0 0 0 0 0 0 1 0 \n", " [7,] 0 0 0 0 0 0 0 0 0 0 0 0 1 \n", " [8,] 0 0 0 0 0 0 0 0 0 0 0 0 0 \n", " [9,] 0 0 0 0 0 0 0 0 0 0 0 0 0 \n", "[10,] 0 0 0 0 0 0 0 0 0 0 0 0 0 \n", " [,14] [,15] [,16]\n", " [1,] 0 0 0 \n", " [2,] 0 0 0 \n", " [3,] 0 0 0 \n", " [4,] 0 0 0 \n", " [5,] 0 0 0 \n", " [6,] 0 0 0 \n", " [7,] 0 0 0 \n", " [8,] 1 0 0 \n", " [9,] 0 1 0 \n", "[10,] 0 0 1 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "y = as.matrix(c(y1, y2))\n", "\n", "X = matrix(0, 10, 4)\n", "X[1:5, 1:2] = X1\n", "X[6:10, 3:4] = X2\n", "\n", "Z = matrix(0, 10, 16)\n", "Z[1:5, 1:8] = Z1\n", "Z[6:10, 9:16] = Z2\n", "Z" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": true }, "outputs": [], "source": [ "results = mme3(y, X, Z, A, G, R)" ] }, { "cell_type": "code", "execution_count": 12, "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", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
0.08341057 0.00000000 -0.03058387 0.00000000 0.000000e+00 0.000000e+00 0.000000e+00 2.780352e-02 0.000000e+00 0.000000e+00 2.780352e-02 2.780352e-02 0.000000e+00 0.000000e+00 0.000000e+00-1.019462e-02 0.000000e+00 0.000000e+00-1.019462e-02-1.019462e-02
0.00000000 0.05560704 0.00000000 -0.02038925 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 2.780352e-02 2.780352e-02 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00-1.019462e-02-1.019462e-02 0.000000e+00 0.000000e+00
-0.03058387 0.00000000 0.11121409 0.00000000 0.000000e+00 0.000000e+00 0.000000e+00-1.019462e-02 0.000000e+00 0.000000e+00-1.019462e-02-1.019462e-02 0.000000e+00 0.000000e+00 0.000000e+00 3.707136e-02 0.000000e+00 0.000000e+00 3.707136e-02 3.707136e-02
0.00000000 -0.02038925 0.00000000 0.07414272 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00-1.019462e-02-1.019462e-02 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 3.707136e-02 3.707136e-02 0.000000e+00 0.000000e+00
0.00000000 0.00000000 0.00000000 0.00000000 1.540616e-01 4.201681e-02 2.332401e-16-5.602241e-02 3.731842e-17-8.403361e-02-4.478211e-16-1.559793e-16-6.932773e-02-1.890756e-02-1.049581e-16 2.521008e-02-1.679329e-17 3.781513e-02 2.015195e-16 7.019070e-17
0.00000000 0.00000000 0.00000000 0.00000000 4.201681e-02 1.680672e-01 4.201681e-02-6.530724e-17-8.403361e-02-8.403361e-02 7.463684e-17-9.621155e-18-1.890756e-02-7.563025e-02-1.890756e-02 2.938826e-17 3.781513e-02 3.781513e-02-3.358658e-17 4.329520e-18
0.00000000 0.00000000 0.00000000 0.00000000 9.329605e-17 4.201681e-02 1.680672e-01 1.865921e-17-8.403361e-02 4.201681e-02-3.265362e-17-8.403361e-02-4.198322e-17-1.890756e-02-7.563025e-02-8.396645e-18 3.781513e-02-1.890756e-02 1.469413e-17 3.781513e-02
0.02780352 0.00000000 -0.01019462 0.00000000 -5.602241e-02-1.399441e-16-1.072905e-16 1.818651e-01 4.201681e-02 5.597763e-17-8.403361e-02 2.128316e-17 2.521008e-02 6.297484e-17 4.828071e-17-7.952236e-02-1.890756e-02-2.518993e-17 3.781513e-02-9.577423e-18
0.00000000 0.02780352 0.00000000 -0.01019462 9.329605e-18-8.403361e-02-8.403361e-02 4.201681e-02 2.378876e-01-7.463684e-17-8.403361e-02-8.163405e-17-4.198322e-18 3.781513e-02 3.781513e-02-1.890756e-02-1.047324e-01 3.358658e-17 3.781513e-02 3.673532e-17
0.00000000 0.02780352 0.00000000 -0.01019462 -8.403361e-02-8.403361e-02 4.201681e-02-9.796086e-17-6.530724e-17 2.378876e-01 9.796086e-17-8.403361e-02 3.781513e-02 3.781513e-02-1.890756e-02 4.408238e-17 2.938826e-17-1.047324e-01-4.408238e-17 3.781513e-02
0.02780352 0.00000000 -0.01019462 0.00000000 4.664803e-18-3.731842e-17-9.329605e-18-8.403361e-02-8.403361e-02-4.664803e-18 1.958707e-01 1.137046e-17-2.099161e-18 1.679329e-17 4.198322e-18 3.781513e-02 3.781513e-02 2.099161e-18-8.582488e-02-5.116705e-18
0.02780352 0.00000000 -0.01019462 0.00000000 -1.020426e-16-3.221629e-16-8.403361e-02-2.183711e-16 9.737775e-17-8.403361e-02 2.040851e-16 1.958707e-01 4.591915e-17 1.449733e-16 3.781513e-02 9.826698e-17-4.381999e-17 3.781513e-02-9.183830e-17-8.582488e-02
0.00000000 0.00000000 0.00000000 0.00000000 -6.932773e-02-1.890756e-02-1.049581e-16 2.521008e-02-1.679329e-17 3.781513e-02 2.015195e-16 7.019070e-17 7.703081e-02 2.100840e-02 1.166201e-16-2.801120e-02 1.865921e-17-4.201681e-02-2.239105e-16-7.798967e-17
0.00000000 0.00000000 0.00000000 0.00000000 -1.890756e-02-7.563025e-02-1.890756e-02 2.938826e-17 3.781513e-02 3.781513e-02-3.358658e-17 4.329520e-18 2.100840e-02 8.403361e-02 2.100840e-02-3.265362e-17-4.201681e-02-4.201681e-02 3.731842e-17-4.810578e-18
0.00000000 0.00000000 0.00000000 0.00000000 -4.198322e-17-1.890756e-02-7.563025e-02-8.396645e-18 3.781513e-02-1.890756e-02 1.469413e-17 3.781513e-02 4.664803e-17 2.100840e-02 8.403361e-02 9.329605e-18-4.201681e-02 2.100840e-02-1.632681e-17-4.201681e-02
-0.01019462 0.00000000 0.03707136 0.00000000 2.521008e-02 6.297484e-17 4.828071e-17-7.952236e-02-1.890756e-02-2.518993e-17 3.781513e-02-9.577423e-18-2.801120e-02-6.997204e-17-5.364523e-17 1.141022e-01 2.100840e-02 2.798882e-17-4.201681e-02 1.064158e-17
0.00000000 -0.01019462 0.00000000 0.03707136 -4.198322e-18 3.781513e-02 3.781513e-02-1.890756e-02-1.047324e-01 3.358658e-17 3.781513e-02 3.673532e-17 4.664803e-18-4.201681e-02-4.201681e-02 2.100840e-02 1.421134e-01-3.731842e-17-4.201681e-02-4.081702e-17
0.00000000 -0.01019462 0.00000000 0.03707136 3.781513e-02 3.781513e-02-1.890756e-02 4.408238e-17 2.938826e-17-1.047324e-01-4.408238e-17 3.781513e-02-4.201681e-02-4.201681e-02 2.100840e-02-4.898043e-17-3.265362e-17 1.421134e-01 4.898043e-17-4.201681e-02
-0.01019462 0.00000000 0.03707136 0.00000000 -2.099161e-18 1.679329e-17 4.198322e-18 3.781513e-02 3.781513e-02 2.099161e-18-8.582488e-02-5.116705e-18 2.332401e-18-1.865921e-17-4.664803e-18-4.201681e-02-4.201681e-02-2.332401e-18 1.211050e-01 5.685228e-18
-0.01019462 0.00000000 0.03707136 0.00000000 4.591915e-17 1.449733e-16 3.781513e-02 9.826698e-17-4.381999e-17 3.781513e-02-9.183830e-17-8.582488e-02-5.102128e-17-1.610815e-16-4.201681e-02-1.091855e-16 4.868888e-17-4.201681e-02 1.020426e-16 1.211050e-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", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
4.360866999
3.397261592
6.799897620
5.880295937
0.150915567
-0.015392510
-0.078391896
-0.010238959
-0.270331441
0.275808258
-0.316117562
0.243755523
0.279597973
-0.007610071
-0.170341439
-0.012670709
-0.477830262
0.517238387
-0.478983695
0.391961543
\n", "
\n", "
\n" ], "text/latex": [ "\\begin{description}\n", "\\item[\\$C] \\begin{tabular}{llllllllllllllllllll}\n", "\t 0.08341057 & 0.00000000 & -0.03058387 & 0.00000000 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 2.780352e-02 & 0.000000e+00 & 0.000000e+00 & 2.780352e-02 & 2.780352e-02 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & -1.019462e-02 & 0.000000e+00 & 0.000000e+00 & -1.019462e-02 & -1.019462e-02\\\\\n", "\t 0.00000000 & 0.05560704 & 0.00000000 & -0.02038925 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 2.780352e-02 & 2.780352e-02 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & -1.019462e-02 & -1.019462e-02 & 0.000000e+00 & 0.000000e+00\\\\\n", "\t -0.03058387 & 0.00000000 & 0.11121409 & 0.00000000 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & -1.019462e-02 & 0.000000e+00 & 0.000000e+00 & -1.019462e-02 & -1.019462e-02 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 3.707136e-02 & 0.000000e+00 & 0.000000e+00 & 3.707136e-02 & 3.707136e-02\\\\\n", "\t 0.00000000 & -0.02038925 & 0.00000000 & 0.07414272 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & -1.019462e-02 & -1.019462e-02 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 3.707136e-02 & 3.707136e-02 & 0.000000e+00 & 0.000000e+00\\\\\n", "\t 0.00000000 & 0.00000000 & 0.00000000 & 0.00000000 & 1.540616e-01 & 4.201681e-02 & 2.332401e-16 & -5.602241e-02 & 3.731842e-17 & -8.403361e-02 & -4.478211e-16 & -1.559793e-16 & -6.932773e-02 & -1.890756e-02 & -1.049581e-16 & 2.521008e-02 & -1.679329e-17 & 3.781513e-02 & 2.015195e-16 & 7.019070e-17\\\\\n", "\t 0.00000000 & 0.00000000 & 0.00000000 & 0.00000000 & 4.201681e-02 & 1.680672e-01 & 4.201681e-02 & -6.530724e-17 & -8.403361e-02 & -8.403361e-02 & 7.463684e-17 & -9.621155e-18 & -1.890756e-02 & -7.563025e-02 & -1.890756e-02 & 2.938826e-17 & 3.781513e-02 & 3.781513e-02 & -3.358658e-17 & 4.329520e-18\\\\\n", "\t 0.00000000 & 0.00000000 & 0.00000000 & 0.00000000 & 9.329605e-17 & 4.201681e-02 & 1.680672e-01 & 1.865921e-17 & -8.403361e-02 & 4.201681e-02 & -3.265362e-17 & -8.403361e-02 & -4.198322e-17 & -1.890756e-02 & -7.563025e-02 & -8.396645e-18 & 3.781513e-02 & -1.890756e-02 & 1.469413e-17 & 3.781513e-02\\\\\n", "\t 0.02780352 & 0.00000000 & -0.01019462 & 0.00000000 & -5.602241e-02 & -1.399441e-16 & -1.072905e-16 & 1.818651e-01 & 4.201681e-02 & 5.597763e-17 & -8.403361e-02 & 2.128316e-17 & 2.521008e-02 & 6.297484e-17 & 4.828071e-17 & -7.952236e-02 & -1.890756e-02 & -2.518993e-17 & 3.781513e-02 & -9.577423e-18\\\\\n", "\t 0.00000000 & 0.02780352 & 0.00000000 & -0.01019462 & 9.329605e-18 & -8.403361e-02 & -8.403361e-02 & 4.201681e-02 & 2.378876e-01 & -7.463684e-17 & -8.403361e-02 & -8.163405e-17 & -4.198322e-18 & 3.781513e-02 & 3.781513e-02 & -1.890756e-02 & -1.047324e-01 & 3.358658e-17 & 3.781513e-02 & 3.673532e-17\\\\\n", "\t 0.00000000 & 0.02780352 & 0.00000000 & -0.01019462 & -8.403361e-02 & -8.403361e-02 & 4.201681e-02 & -9.796086e-17 & -6.530724e-17 & 2.378876e-01 & 9.796086e-17 & -8.403361e-02 & 3.781513e-02 & 3.781513e-02 & -1.890756e-02 & 4.408238e-17 & 2.938826e-17 & -1.047324e-01 & -4.408238e-17 & 3.781513e-02\\\\\n", "\t 0.02780352 & 0.00000000 & -0.01019462 & 0.00000000 & 4.664803e-18 & -3.731842e-17 & -9.329605e-18 & -8.403361e-02 & -8.403361e-02 & -4.664803e-18 & 1.958707e-01 & 1.137046e-17 & -2.099161e-18 & 1.679329e-17 & 4.198322e-18 & 3.781513e-02 & 3.781513e-02 & 2.099161e-18 & -8.582488e-02 & -5.116705e-18\\\\\n", "\t 0.02780352 & 0.00000000 & -0.01019462 & 0.00000000 & -1.020426e-16 & -3.221629e-16 & -8.403361e-02 & -2.183711e-16 & 9.737775e-17 & -8.403361e-02 & 2.040851e-16 & 1.958707e-01 & 4.591915e-17 & 1.449733e-16 & 3.781513e-02 & 9.826698e-17 & -4.381999e-17 & 3.781513e-02 & -9.183830e-17 & -8.582488e-02\\\\\n", "\t 0.00000000 & 0.00000000 & 0.00000000 & 0.00000000 & -6.932773e-02 & -1.890756e-02 & -1.049581e-16 & 2.521008e-02 & -1.679329e-17 & 3.781513e-02 & 2.015195e-16 & 7.019070e-17 & 7.703081e-02 & 2.100840e-02 & 1.166201e-16 & -2.801120e-02 & 1.865921e-17 & -4.201681e-02 & -2.239105e-16 & -7.798967e-17\\\\\n", "\t 0.00000000 & 0.00000000 & 0.00000000 & 0.00000000 & -1.890756e-02 & -7.563025e-02 & -1.890756e-02 & 2.938826e-17 & 3.781513e-02 & 3.781513e-02 & -3.358658e-17 & 4.329520e-18 & 2.100840e-02 & 8.403361e-02 & 2.100840e-02 & -3.265362e-17 & -4.201681e-02 & -4.201681e-02 & 3.731842e-17 & -4.810578e-18\\\\\n", "\t 0.00000000 & 0.00000000 & 0.00000000 & 0.00000000 & -4.198322e-17 & -1.890756e-02 & -7.563025e-02 & -8.396645e-18 & 3.781513e-02 & -1.890756e-02 & 1.469413e-17 & 3.781513e-02 & 4.664803e-17 & 2.100840e-02 & 8.403361e-02 & 9.329605e-18 & -4.201681e-02 & 2.100840e-02 & -1.632681e-17 & -4.201681e-02\\\\\n", "\t -0.01019462 & 0.00000000 & 0.03707136 & 0.00000000 & 2.521008e-02 & 6.297484e-17 & 4.828071e-17 & -7.952236e-02 & -1.890756e-02 & -2.518993e-17 & 3.781513e-02 & -9.577423e-18 & -2.801120e-02 & -6.997204e-17 & -5.364523e-17 & 1.141022e-01 & 2.100840e-02 & 2.798882e-17 & -4.201681e-02 & 1.064158e-17\\\\\n", "\t 0.00000000 & -0.01019462 & 0.00000000 & 0.03707136 & -4.198322e-18 & 3.781513e-02 & 3.781513e-02 & -1.890756e-02 & -1.047324e-01 & 3.358658e-17 & 3.781513e-02 & 3.673532e-17 & 4.664803e-18 & -4.201681e-02 & -4.201681e-02 & 2.100840e-02 & 1.421134e-01 & -3.731842e-17 & -4.201681e-02 & -4.081702e-17\\\\\n", "\t 0.00000000 & -0.01019462 & 0.00000000 & 0.03707136 & 3.781513e-02 & 3.781513e-02 & -1.890756e-02 & 4.408238e-17 & 2.938826e-17 & -1.047324e-01 & -4.408238e-17 & 3.781513e-02 & -4.201681e-02 & -4.201681e-02 & 2.100840e-02 & -4.898043e-17 & -3.265362e-17 & 1.421134e-01 & 4.898043e-17 & -4.201681e-02\\\\\n", "\t -0.01019462 & 0.00000000 & 0.03707136 & 0.00000000 & -2.099161e-18 & 1.679329e-17 & 4.198322e-18 & 3.781513e-02 & 3.781513e-02 & 2.099161e-18 & -8.582488e-02 & -5.116705e-18 & 2.332401e-18 & -1.865921e-17 & -4.664803e-18 & -4.201681e-02 & -4.201681e-02 & -2.332401e-18 & 1.211050e-01 & 5.685228e-18\\\\\n", "\t -0.01019462 & 0.00000000 & 0.03707136 & 0.00000000 & 4.591915e-17 & 1.449733e-16 & 3.781513e-02 & 9.826698e-17 & -4.381999e-17 & 3.781513e-02 & -9.183830e-17 & -8.582488e-02 & -5.102128e-17 & -1.610815e-16 & -4.201681e-02 & -1.091855e-16 & 4.868888e-17 & -4.201681e-02 & 1.020426e-16 & 1.211050e-01\\\\\n", "\\end{tabular}\n", "\n", "\\item[\\$est] \\begin{tabular}{l}\n", "\t 4.360866999\\\\\n", "\t 3.397261592\\\\\n", "\t 6.799897620\\\\\n", "\t 5.880295937\\\\\n", "\t 0.150915567\\\\\n", "\t -0.015392510\\\\\n", "\t -0.078391896\\\\\n", "\t -0.010238959\\\\\n", "\t -0.270331441\\\\\n", "\t 0.275808258\\\\\n", "\t -0.316117562\\\\\n", "\t 0.243755523\\\\\n", "\t 0.279597973\\\\\n", "\t -0.007610071\\\\\n", "\t -0.170341439\\\\\n", "\t -0.012670709\\\\\n", "\t -0.477830262\\\\\n", "\t 0.517238387\\\\\n", "\t -0.478983695\\\\\n", "\t 0.391961543\\\\\n", "\\end{tabular}\n", "\n", "\\end{description}\n" ], "text/markdown": [ "$C\n", ": \n", "| 0.08341057 | 0.00000000 | -0.03058387 | 0.00000000 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.780352e-02 | 0.000000e+00 | 0.000000e+00 | 2.780352e-02 | 2.780352e-02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | -1.019462e-02 | 0.000000e+00 | 0.000000e+00 | -1.019462e-02 | -1.019462e-02 |\n", "| 0.00000000 | 0.05560704 | 0.00000000 | -0.02038925 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.780352e-02 | 2.780352e-02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | -1.019462e-02 | -1.019462e-02 | 0.000000e+00 | 0.000000e+00 |\n", "| -0.03058387 | 0.00000000 | 0.11121409 | 0.00000000 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | -1.019462e-02 | 0.000000e+00 | 0.000000e+00 | -1.019462e-02 | -1.019462e-02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 3.707136e-02 | 0.000000e+00 | 0.000000e+00 | 3.707136e-02 | 3.707136e-02 |\n", "| 0.00000000 | -0.02038925 | 0.00000000 | 0.07414272 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | -1.019462e-02 | -1.019462e-02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 3.707136e-02 | 3.707136e-02 | 0.000000e+00 | 0.000000e+00 |\n", "| 0.00000000 | 0.00000000 | 0.00000000 | 0.00000000 | 1.540616e-01 | 4.201681e-02 | 2.332401e-16 | -5.602241e-02 | 3.731842e-17 | -8.403361e-02 | -4.478211e-16 | -1.559793e-16 | -6.932773e-02 | -1.890756e-02 | -1.049581e-16 | 2.521008e-02 | -1.679329e-17 | 3.781513e-02 | 2.015195e-16 | 7.019070e-17 |\n", "| 0.00000000 | 0.00000000 | 0.00000000 | 0.00000000 | 4.201681e-02 | 1.680672e-01 | 4.201681e-02 | -6.530724e-17 | -8.403361e-02 | -8.403361e-02 | 7.463684e-17 | -9.621155e-18 | -1.890756e-02 | -7.563025e-02 | -1.890756e-02 | 2.938826e-17 | 3.781513e-02 | 3.781513e-02 | -3.358658e-17 | 4.329520e-18 |\n", "| 0.00000000 | 0.00000000 | 0.00000000 | 0.00000000 | 9.329605e-17 | 4.201681e-02 | 1.680672e-01 | 1.865921e-17 | -8.403361e-02 | 4.201681e-02 | -3.265362e-17 | -8.403361e-02 | -4.198322e-17 | -1.890756e-02 | -7.563025e-02 | -8.396645e-18 | 3.781513e-02 | -1.890756e-02 | 1.469413e-17 | 3.781513e-02 |\n", "| 0.02780352 | 0.00000000 | -0.01019462 | 0.00000000 | -5.602241e-02 | -1.399441e-16 | -1.072905e-16 | 1.818651e-01 | 4.201681e-02 | 5.597763e-17 | -8.403361e-02 | 2.128316e-17 | 2.521008e-02 | 6.297484e-17 | 4.828071e-17 | -7.952236e-02 | -1.890756e-02 | -2.518993e-17 | 3.781513e-02 | -9.577423e-18 |\n", "| 0.00000000 | 0.02780352 | 0.00000000 | -0.01019462 | 9.329605e-18 | -8.403361e-02 | -8.403361e-02 | 4.201681e-02 | 2.378876e-01 | -7.463684e-17 | -8.403361e-02 | -8.163405e-17 | -4.198322e-18 | 3.781513e-02 | 3.781513e-02 | -1.890756e-02 | -1.047324e-01 | 3.358658e-17 | 3.781513e-02 | 3.673532e-17 |\n", "| 0.00000000 | 0.02780352 | 0.00000000 | -0.01019462 | -8.403361e-02 | -8.403361e-02 | 4.201681e-02 | -9.796086e-17 | -6.530724e-17 | 2.378876e-01 | 9.796086e-17 | -8.403361e-02 | 3.781513e-02 | 3.781513e-02 | -1.890756e-02 | 4.408238e-17 | 2.938826e-17 | -1.047324e-01 | -4.408238e-17 | 3.781513e-02 |\n", "| 0.02780352 | 0.00000000 | -0.01019462 | 0.00000000 | 4.664803e-18 | -3.731842e-17 | -9.329605e-18 | -8.403361e-02 | -8.403361e-02 | -4.664803e-18 | 1.958707e-01 | 1.137046e-17 | -2.099161e-18 | 1.679329e-17 | 4.198322e-18 | 3.781513e-02 | 3.781513e-02 | 2.099161e-18 | -8.582488e-02 | -5.116705e-18 |\n", "| 0.02780352 | 0.00000000 | -0.01019462 | 0.00000000 | -1.020426e-16 | -3.221629e-16 | -8.403361e-02 | -2.183711e-16 | 9.737775e-17 | -8.403361e-02 | 2.040851e-16 | 1.958707e-01 | 4.591915e-17 | 1.449733e-16 | 3.781513e-02 | 9.826698e-17 | -4.381999e-17 | 3.781513e-02 | -9.183830e-17 | -8.582488e-02 |\n", "| 0.00000000 | 0.00000000 | 0.00000000 | 0.00000000 | -6.932773e-02 | -1.890756e-02 | -1.049581e-16 | 2.521008e-02 | -1.679329e-17 | 3.781513e-02 | 2.015195e-16 | 7.019070e-17 | 7.703081e-02 | 2.100840e-02 | 1.166201e-16 | -2.801120e-02 | 1.865921e-17 | -4.201681e-02 | -2.239105e-16 | -7.798967e-17 |\n", "| 0.00000000 | 0.00000000 | 0.00000000 | 0.00000000 | -1.890756e-02 | -7.563025e-02 | -1.890756e-02 | 2.938826e-17 | 3.781513e-02 | 3.781513e-02 | -3.358658e-17 | 4.329520e-18 | 2.100840e-02 | 8.403361e-02 | 2.100840e-02 | -3.265362e-17 | -4.201681e-02 | -4.201681e-02 | 3.731842e-17 | -4.810578e-18 |\n", "| 0.00000000 | 0.00000000 | 0.00000000 | 0.00000000 | -4.198322e-17 | -1.890756e-02 | -7.563025e-02 | -8.396645e-18 | 3.781513e-02 | -1.890756e-02 | 1.469413e-17 | 3.781513e-02 | 4.664803e-17 | 2.100840e-02 | 8.403361e-02 | 9.329605e-18 | -4.201681e-02 | 2.100840e-02 | -1.632681e-17 | -4.201681e-02 |\n", "| -0.01019462 | 0.00000000 | 0.03707136 | 0.00000000 | 2.521008e-02 | 6.297484e-17 | 4.828071e-17 | -7.952236e-02 | -1.890756e-02 | -2.518993e-17 | 3.781513e-02 | -9.577423e-18 | -2.801120e-02 | -6.997204e-17 | -5.364523e-17 | 1.141022e-01 | 2.100840e-02 | 2.798882e-17 | -4.201681e-02 | 1.064158e-17 |\n", "| 0.00000000 | -0.01019462 | 0.00000000 | 0.03707136 | -4.198322e-18 | 3.781513e-02 | 3.781513e-02 | -1.890756e-02 | -1.047324e-01 | 3.358658e-17 | 3.781513e-02 | 3.673532e-17 | 4.664803e-18 | -4.201681e-02 | -4.201681e-02 | 2.100840e-02 | 1.421134e-01 | -3.731842e-17 | -4.201681e-02 | -4.081702e-17 |\n", "| 0.00000000 | -0.01019462 | 0.00000000 | 0.03707136 | 3.781513e-02 | 3.781513e-02 | -1.890756e-02 | 4.408238e-17 | 2.938826e-17 | -1.047324e-01 | -4.408238e-17 | 3.781513e-02 | -4.201681e-02 | -4.201681e-02 | 2.100840e-02 | -4.898043e-17 | -3.265362e-17 | 1.421134e-01 | 4.898043e-17 | -4.201681e-02 |\n", "| -0.01019462 | 0.00000000 | 0.03707136 | 0.00000000 | -2.099161e-18 | 1.679329e-17 | 4.198322e-18 | 3.781513e-02 | 3.781513e-02 | 2.099161e-18 | -8.582488e-02 | -5.116705e-18 | 2.332401e-18 | -1.865921e-17 | -4.664803e-18 | -4.201681e-02 | -4.201681e-02 | -2.332401e-18 | 1.211050e-01 | 5.685228e-18 |\n", "| -0.01019462 | 0.00000000 | 0.03707136 | 0.00000000 | 4.591915e-17 | 1.449733e-16 | 3.781513e-02 | 9.826698e-17 | -4.381999e-17 | 3.781513e-02 | -9.183830e-17 | -8.582488e-02 | -5.102128e-17 | -1.610815e-16 | -4.201681e-02 | -1.091855e-16 | 4.868888e-17 | -4.201681e-02 | 1.020426e-16 | 1.211050e-01 |\n", "\n", "\n", "$est\n", ": \n", "| 4.360866999 |\n", "| 3.397261592 |\n", "| 6.799897620 |\n", "| 5.880295937 |\n", "| 0.150915567 |\n", "| -0.015392510 |\n", "| -0.078391896 |\n", "| -0.010238959 |\n", "| -0.270331441 |\n", "| 0.275808258 |\n", "| -0.316117562 |\n", "| 0.243755523 |\n", "| 0.279597973 |\n", "| -0.007610071 |\n", "| -0.170341439 |\n", "| -0.012670709 |\n", "| -0.477830262 |\n", "| 0.517238387 |\n", "| -0.478983695 |\n", "| 0.391961543 |\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "$C\n", " [,1] [,2] [,3] [,4] [,5]\n", " [1,] 0.08341057 0.00000000 -0.03058387 0.00000000 0.000000e+00\n", " [2,] 0.00000000 0.05560704 0.00000000 -0.02038925 0.000000e+00\n", " [3,] -0.03058387 0.00000000 0.11121409 0.00000000 0.000000e+00\n", " [4,] 0.00000000 -0.02038925 0.00000000 0.07414272 0.000000e+00\n", " [5,] 0.00000000 0.00000000 0.00000000 0.00000000 1.540616e-01\n", " [6,] 0.00000000 0.00000000 0.00000000 0.00000000 4.201681e-02\n", " [7,] 0.00000000 0.00000000 0.00000000 0.00000000 9.329605e-17\n", " [8,] 0.02780352 0.00000000 -0.01019462 0.00000000 -5.602241e-02\n", " [9,] 0.00000000 0.02780352 0.00000000 -0.01019462 9.329605e-18\n", "[10,] 0.00000000 0.02780352 0.00000000 -0.01019462 -8.403361e-02\n", "[11,] 0.02780352 0.00000000 -0.01019462 0.00000000 4.664803e-18\n", "[12,] 0.02780352 0.00000000 -0.01019462 0.00000000 -1.020426e-16\n", "[13,] 0.00000000 0.00000000 0.00000000 0.00000000 -6.932773e-02\n", "[14,] 0.00000000 0.00000000 0.00000000 0.00000000 -1.890756e-02\n", "[15,] 0.00000000 0.00000000 0.00000000 0.00000000 -4.198322e-17\n", "[16,] -0.01019462 0.00000000 0.03707136 0.00000000 2.521008e-02\n", "[17,] 0.00000000 -0.01019462 0.00000000 0.03707136 -4.198322e-18\n", "[18,] 0.00000000 -0.01019462 0.00000000 0.03707136 3.781513e-02\n", "[19,] -0.01019462 0.00000000 0.03707136 0.00000000 -2.099161e-18\n", "[20,] -0.01019462 0.00000000 0.03707136 0.00000000 4.591915e-17\n", " [,6] [,7] [,8] [,9] [,10]\n", " [1,] 0.000000e+00 0.000000e+00 2.780352e-02 0.000000e+00 0.000000e+00\n", " [2,] 0.000000e+00 0.000000e+00 0.000000e+00 2.780352e-02 2.780352e-02\n", " [3,] 0.000000e+00 0.000000e+00 -1.019462e-02 0.000000e+00 0.000000e+00\n", " [4,] 0.000000e+00 0.000000e+00 0.000000e+00 -1.019462e-02 -1.019462e-02\n", " [5,] 4.201681e-02 2.332401e-16 -5.602241e-02 3.731842e-17 -8.403361e-02\n", " [6,] 1.680672e-01 4.201681e-02 -6.530724e-17 -8.403361e-02 -8.403361e-02\n", " [7,] 4.201681e-02 1.680672e-01 1.865921e-17 -8.403361e-02 4.201681e-02\n", " [8,] -1.399441e-16 -1.072905e-16 1.818651e-01 4.201681e-02 5.597763e-17\n", " [9,] -8.403361e-02 -8.403361e-02 4.201681e-02 2.378876e-01 -7.463684e-17\n", "[10,] -8.403361e-02 4.201681e-02 -9.796086e-17 -6.530724e-17 2.378876e-01\n", "[11,] -3.731842e-17 -9.329605e-18 -8.403361e-02 -8.403361e-02 -4.664803e-18\n", "[12,] -3.221629e-16 -8.403361e-02 -2.183711e-16 9.737775e-17 -8.403361e-02\n", "[13,] -1.890756e-02 -1.049581e-16 2.521008e-02 -1.679329e-17 3.781513e-02\n", "[14,] -7.563025e-02 -1.890756e-02 2.938826e-17 3.781513e-02 3.781513e-02\n", "[15,] -1.890756e-02 -7.563025e-02 -8.396645e-18 3.781513e-02 -1.890756e-02\n", "[16,] 6.297484e-17 4.828071e-17 -7.952236e-02 -1.890756e-02 -2.518993e-17\n", "[17,] 3.781513e-02 3.781513e-02 -1.890756e-02 -1.047324e-01 3.358658e-17\n", "[18,] 3.781513e-02 -1.890756e-02 4.408238e-17 2.938826e-17 -1.047324e-01\n", "[19,] 1.679329e-17 4.198322e-18 3.781513e-02 3.781513e-02 2.099161e-18\n", "[20,] 1.449733e-16 3.781513e-02 9.826698e-17 -4.381999e-17 3.781513e-02\n", " [,11] [,12] [,13] [,14] [,15]\n", " [1,] 2.780352e-02 2.780352e-02 0.000000e+00 0.000000e+00 0.000000e+00\n", " [2,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n", " [3,] -1.019462e-02 -1.019462e-02 0.000000e+00 0.000000e+00 0.000000e+00\n", " [4,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n", " [5,] -4.478211e-16 -1.559793e-16 -6.932773e-02 -1.890756e-02 -1.049581e-16\n", " [6,] 7.463684e-17 -9.621155e-18 -1.890756e-02 -7.563025e-02 -1.890756e-02\n", " [7,] -3.265362e-17 -8.403361e-02 -4.198322e-17 -1.890756e-02 -7.563025e-02\n", " [8,] -8.403361e-02 2.128316e-17 2.521008e-02 6.297484e-17 4.828071e-17\n", " [9,] -8.403361e-02 -8.163405e-17 -4.198322e-18 3.781513e-02 3.781513e-02\n", "[10,] 9.796086e-17 -8.403361e-02 3.781513e-02 3.781513e-02 -1.890756e-02\n", "[11,] 1.958707e-01 1.137046e-17 -2.099161e-18 1.679329e-17 4.198322e-18\n", "[12,] 2.040851e-16 1.958707e-01 4.591915e-17 1.449733e-16 3.781513e-02\n", "[13,] 2.015195e-16 7.019070e-17 7.703081e-02 2.100840e-02 1.166201e-16\n", "[14,] -3.358658e-17 4.329520e-18 2.100840e-02 8.403361e-02 2.100840e-02\n", "[15,] 1.469413e-17 3.781513e-02 4.664803e-17 2.100840e-02 8.403361e-02\n", "[16,] 3.781513e-02 -9.577423e-18 -2.801120e-02 -6.997204e-17 -5.364523e-17\n", "[17,] 3.781513e-02 3.673532e-17 4.664803e-18 -4.201681e-02 -4.201681e-02\n", "[18,] -4.408238e-17 3.781513e-02 -4.201681e-02 -4.201681e-02 2.100840e-02\n", "[19,] -8.582488e-02 -5.116705e-18 2.332401e-18 -1.865921e-17 -4.664803e-18\n", "[20,] -9.183830e-17 -8.582488e-02 -5.102128e-17 -1.610815e-16 -4.201681e-02\n", " [,16] [,17] [,18] [,19] [,20]\n", " [1,] -1.019462e-02 0.000000e+00 0.000000e+00 -1.019462e-02 -1.019462e-02\n", " [2,] 0.000000e+00 -1.019462e-02 -1.019462e-02 0.000000e+00 0.000000e+00\n", " [3,] 3.707136e-02 0.000000e+00 0.000000e+00 3.707136e-02 3.707136e-02\n", " [4,] 0.000000e+00 3.707136e-02 3.707136e-02 0.000000e+00 0.000000e+00\n", " [5,] 2.521008e-02 -1.679329e-17 3.781513e-02 2.015195e-16 7.019070e-17\n", " [6,] 2.938826e-17 3.781513e-02 3.781513e-02 -3.358658e-17 4.329520e-18\n", " [7,] -8.396645e-18 3.781513e-02 -1.890756e-02 1.469413e-17 3.781513e-02\n", " [8,] -7.952236e-02 -1.890756e-02 -2.518993e-17 3.781513e-02 -9.577423e-18\n", " [9,] -1.890756e-02 -1.047324e-01 3.358658e-17 3.781513e-02 3.673532e-17\n", "[10,] 4.408238e-17 2.938826e-17 -1.047324e-01 -4.408238e-17 3.781513e-02\n", "[11,] 3.781513e-02 3.781513e-02 2.099161e-18 -8.582488e-02 -5.116705e-18\n", "[12,] 9.826698e-17 -4.381999e-17 3.781513e-02 -9.183830e-17 -8.582488e-02\n", "[13,] -2.801120e-02 1.865921e-17 -4.201681e-02 -2.239105e-16 -7.798967e-17\n", "[14,] -3.265362e-17 -4.201681e-02 -4.201681e-02 3.731842e-17 -4.810578e-18\n", "[15,] 9.329605e-18 -4.201681e-02 2.100840e-02 -1.632681e-17 -4.201681e-02\n", "[16,] 1.141022e-01 2.100840e-02 2.798882e-17 -4.201681e-02 1.064158e-17\n", "[17,] 2.100840e-02 1.421134e-01 -3.731842e-17 -4.201681e-02 -4.081702e-17\n", "[18,] -4.898043e-17 -3.265362e-17 1.421134e-01 4.898043e-17 -4.201681e-02\n", "[19,] -4.201681e-02 -4.201681e-02 -2.332401e-18 1.211050e-01 5.685228e-18\n", "[20,] -1.091855e-16 4.868888e-17 -4.201681e-02 1.020426e-16 1.211050e-01\n", "\n", "$est\n", " [,1]\n", " [1,] 4.360866999\n", " [2,] 3.397261592\n", " [3,] 6.799897620\n", " [4,] 5.880295937\n", " [5,] 0.150915567\n", " [6,] -0.015392510\n", " [7,] -0.078391896\n", " [8,] -0.010238959\n", " [9,] -0.270331441\n", "[10,] 0.275808258\n", "[11,] -0.316117562\n", "[12,] 0.243755523\n", "[13,] 0.279597973\n", "[14,] -0.007610071\n", "[15,] -0.170341439\n", "[16,] -0.012670709\n", "[17,] -0.477830262\n", "[18,] 0.517238387\n", "[19,] -0.478983695\n", "[20,] 0.391961543\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "results" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 20
  2. \n", "\t
  3. 20
  4. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 20\n", "\\item 20\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 20\n", "2. 20\n", "\n", "\n" ], "text/plain": [ "[1] 20 20" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "dim(results$C)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 20
  2. \n", "\t
  3. 1
  4. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 20\n", "\\item 1\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 20\n", "2. 1\n", "\n", "\n" ], "text/plain": [ "[1] 20 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "dim(results$est)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 0.150915567312792
  2. \n", "\t
  3. -0.0153925096914866
  4. \n", "\t
  5. -0.0783918961641912
  6. \n", "\t
  7. -0.0102389585292879
  8. \n", "\t
  9. -0.270331441389241
  10. \n", "\t
  11. 0.275808257580568
  12. \n", "\t
  13. -0.316117561639436
  14. \n", "\t
  15. 0.243755523005399
  16. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.150915567312792\n", "\\item -0.0153925096914866\n", "\\item -0.0783918961641912\n", "\\item -0.0102389585292879\n", "\\item -0.270331441389241\n", "\\item 0.275808257580568\n", "\\item -0.316117561639436\n", "\\item 0.243755523005399\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.150915567312792\n", "2. -0.0153925096914866\n", "3. -0.0783918961641912\n", "4. -0.0102389585292879\n", "5. -0.270331441389241\n", "6. 0.275808257580568\n", "7. -0.316117561639436\n", "8. 0.243755523005399\n", "\n", "\n" ], "text/plain": [ "[1] 0.15091557 -0.01539251 -0.07839190 -0.01023896 -0.27033144 0.27580826\n", "[7] -0.31611756 0.24375552" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 0.27959797274487
  2. \n", "\t
  3. -0.00761007081989917
  4. \n", "\t
  5. -0.170341438593946
  6. \n", "\t
  7. -0.0126707086241202
  8. \n", "\t
  9. -0.477830261602723
  10. \n", "\t
  11. 0.517238387038388
  12. \n", "\t
  13. -0.478983695055083
  14. \n", "\t
  15. 0.391961542524088
  16. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.27959797274487\n", "\\item -0.00761007081989917\n", "\\item -0.170341438593946\n", "\\item -0.0126707086241202\n", "\\item -0.477830261602723\n", "\\item 0.517238387038388\n", "\\item -0.478983695055083\n", "\\item 0.391961542524088\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.27959797274487\n", "2. -0.00761007081989917\n", "3. -0.170341438593946\n", "4. -0.0126707086241202\n", "5. -0.477830261602723\n", "6. 0.517238387038388\n", "7. -0.478983695055083\n", "8. 0.391961542524088\n", "\n", "\n" ], "text/plain": [ "[1] 0.279597973 -0.007610071 -0.170341439 -0.012670709 -0.477830262\n", "[6] 0.517238387 -0.478983695 0.391961543" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "(est_1 = results$est[5:12])\n", "(est_2 = results$est[13:20])" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "est_1_new = 10*(est_1-mean(est_1))/sd(est_1) + 100\n", "est_2_new = 10*(est_2-mean(est_2))/sd(est_2) + 100" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "scrolled": true }, "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" } ], "source": [ "mean(est_1_new)\n", "sd(est_1_new)" ] }, { "cell_type": "code", "execution_count": 20, "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", "
est_1_newest_2_new
106.98465107.31162
99.41299 99.65949
96.54476 95.32381
99.64762 99.52466
87.80615 87.13134
112.67075113.64312
85.72161 87.10061
111.21146110.30535
\n" ], "text/latex": [ "\\begin{tabular}{ll}\n", " est\\_1\\_new & est\\_2\\_new\\\\\n", "\\hline\n", "\t 106.98465 & 107.31162\\\\\n", "\t 99.41299 & 99.65949\\\\\n", "\t 96.54476 & 95.32381\\\\\n", "\t 99.64762 & 99.52466\\\\\n", "\t 87.80615 & 87.13134\\\\\n", "\t 112.67075 & 113.64312\\\\\n", "\t 85.72161 & 87.10061\\\\\n", "\t 111.21146 & 110.30535\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| est_1_new | est_2_new |\n", "|---|---|\n", "| 106.98465 | 107.31162 |\n", "| 99.41299 | 99.65949 |\n", "| 96.54476 | 95.32381 |\n", "| 99.64762 | 99.52466 |\n", "| 87.80615 | 87.13134 |\n", "| 112.67075 | 113.64312 |\n", "| 85.72161 | 87.10061 |\n", "| 111.21146 | 110.30535 |\n", "\n" ], "text/plain": [ " est_1_new est_2_new\n", "[1,] 106.98465 107.31162\n", "[2,] 99.41299 99.65949\n", "[3,] 96.54476 95.32381\n", "[4,] 99.64762 99.52466\n", "[5,] 87.80615 87.13134\n", "[6,] 112.67075 113.64312\n", "[7,] 85.72161 87.10061\n", "[8,] 111.21146 110.30535" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "0.995859246739469" ], "text/latex": [ "0.995859246739469" ], "text/markdown": [ "0.995859246739469" ], "text/plain": [ "[1] 0.9958592" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cbind(est_1_new, est_2_new)\n", "cor(est_1_new, est_2_new)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dokładność oceny wielocechowej" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "
    \n", "\t
  1. 0.0697632637735818
  2. \n", "\t
  3. 0.0201812427467519
  4. \n", "\t
  5. 0.105355689244772
  6. \n", "\t
  7. 0.174767187326834
  8. \n", "\t
  9. 0.172946680277051
  10. \n", "\t
  11. 0.142427497991221
  12. \n", "\t
  13. 0.144248005041
  14. \n", "\t
  15. 0.185780891703133
  16. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.0697632637735818\n", "\\item 0.0201812427467519\n", "\\item 0.105355689244772\n", "\\item 0.174767187326834\n", "\\item 0.172946680277051\n", "\\item 0.142427497991221\n", "\\item 0.144248005041\n", "\\item 0.185780891703133\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.0697632637735818\n", "2. 0.0201812427467519\n", "3. 0.105355689244772\n", "4. 0.174767187326834\n", "5. 0.172946680277051\n", "6. 0.142427497991221\n", "7. 0.144248005041\n", "8. 0.185780891703133\n", "\n", "\n" ], "text/plain": [ "[1] 0.06976326 0.02018124 0.10535569 0.17476719 0.17294668 0.14242750 0.14424801\n", "[8] 0.18578089" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 0.102455353618369
  2. \n", "\t
  3. 0.0308092271822838
  4. \n", "\t
  5. 0.155018905719145
  6. \n", "\t
  7. 0.25687708309551
  8. \n", "\t
  9. 0.253384555790559
  10. \n", "\t
  11. 0.212418627668356
  12. \n", "\t
  13. 0.215911154973305
  14. \n", "\t
  15. 0.271012693707551
  16. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.102455353618369\n", "\\item 0.0308092271822838\n", "\\item 0.155018905719145\n", "\\item 0.25687708309551\n", "\\item 0.253384555790559\n", "\\item 0.212418627668356\n", "\\item 0.215911154973305\n", "\\item 0.271012693707551\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.102455353618369\n", "2. 0.0308092271822838\n", "3. 0.155018905719145\n", "4. 0.25687708309551\n", "5. 0.253384555790559\n", "6. 0.212418627668356\n", "7. 0.215911154973305\n", "8. 0.271012693707551\n", "\n", "\n" ], "text/plain": [ "[1] 0.10245535 0.03080923 0.15501891 0.25687708 0.25338456 0.21241863 0.21591115\n", "[8] 0.27101269" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 0.264127362788451
  2. \n", "\t
  3. 0.142060700923063
  4. \n", "\t
  5. 0.324585411324618
  6. \n", "\t
  7. 0.418051656290025
  8. \n", "\t
  9. 0.415868585345239
  10. \n", "\t
  11. 0.377395678289009
  12. \n", "\t
  13. 0.37979995397709
  14. \n", "\t
  15. 0.431023075604002
  16. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.264127362788451\n", "\\item 0.142060700923063\n", "\\item 0.324585411324618\n", "\\item 0.418051656290025\n", "\\item 0.415868585345239\n", "\\item 0.377395678289009\n", "\\item 0.37979995397709\n", "\\item 0.431023075604002\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.264127362788451\n", "2. 0.142060700923063\n", "3. 0.324585411324618\n", "4. 0.418051656290025\n", "5. 0.415868585345239\n", "6. 0.377395678289009\n", "7. 0.37979995397709\n", "8. 0.431023075604002\n", "\n", "\n" ], "text/plain": [ "[1] 0.2641274 0.1420607 0.3245854 0.4180517 0.4158686 0.3773957 0.3798000\n", "[8] 0.4310231" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
    \n", "\t
  1. 0.320086478343539
  2. \n", "\t
  3. 0.175525574154548
  4. \n", "\t
  5. 0.393724403255812
  6. \n", "\t
  7. 0.506830428344145
  8. \n", "\t
  9. 0.503373177464353
  10. \n", "\t
  11. 0.460888953727854
  12. \n", "\t
  13. 0.464662409683961
  14. \n", "\t
  15. 0.520588795218982
  16. \n", "
\n" ], "text/latex": [ "\\begin{enumerate*}\n", "\\item 0.320086478343539\n", "\\item 0.175525574154548\n", "\\item 0.393724403255812\n", "\\item 0.506830428344145\n", "\\item 0.503373177464353\n", "\\item 0.460888953727854\n", "\\item 0.464662409683961\n", "\\item 0.520588795218982\n", "\\end{enumerate*}\n" ], "text/markdown": [ "1. 0.320086478343539\n", "2. 0.175525574154548\n", "3. 0.393724403255812\n", "4. 0.506830428344145\n", "5. 0.503373177464353\n", "6. 0.460888953727854\n", "7. 0.464662409683961\n", "8. 0.520588795218982\n", "\n", "\n" ], "text/plain": [ "[1] 0.3200865 0.1755256 0.3937244 0.5068304 0.5033732 0.4608890 0.4646624\n", "[8] 0.5205888" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "C = as.matrix(mme3(y, X, Z, A, G, R)$C)\n", "invC = ginv(C)\n", "\n", "invC22 = invC[5:20, 5:20]\n", "\n", "trait1 = diag(invC22)[1:8]\n", "trait2 = diag(invC22)[9:16]\n", "\n", "(r2_1 = (20-trait1) / 20)\n", "(r2_2 = (40-trait2) / 40)\n", "\n", "(r1 = sqrt(r2_1))\n", "(r2 = sqrt(r2_2))" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\t\n", "\n", "
0.2640.320
0.1420.176
0.3250.394
0.4180.507
0.4160.503
0.3770.461
0.3800.465
0.4310.521
\n" ], "text/latex": [ "\\begin{tabular}{ll}\n", "\t 0.264 & 0.320\\\\\n", "\t 0.142 & 0.176\\\\\n", "\t 0.325 & 0.394\\\\\n", "\t 0.418 & 0.507\\\\\n", "\t 0.416 & 0.503\\\\\n", "\t 0.377 & 0.461\\\\\n", "\t 0.380 & 0.465\\\\\n", "\t 0.431 & 0.521\\\\\n", "\\end{tabular}\n" ], "text/markdown": [ "\n", "| 0.264 | 0.320 |\n", "| 0.142 | 0.176 |\n", "| 0.325 | 0.394 |\n", "| 0.418 | 0.507 |\n", "| 0.416 | 0.503 |\n", "| 0.377 | 0.461 |\n", "| 0.380 | 0.465 |\n", "| 0.431 | 0.521 |\n", "\n" ], "text/plain": [ " [,1] [,2] \n", "[1,] 0.264 0.320\n", "[2,] 0.142 0.176\n", "[3,] 0.325 0.394\n", "[4,] 0.418 0.507\n", "[5,] 0.416 0.503\n", "[6,] 0.377 0.461\n", "[7,] 0.380 0.465\n", "[8,] 0.431 0.521" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cbind(round(r1, digits = 3), round(r2, digits = 3))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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 }