######################Germination of Orobanche######################
######################Modeling Dispersion######################
#Orobanche is a genus of parasitic plants without chlorophyll that grows on the
#roots of flowering plants. An experiment was made where a bach of seeds of the
#species Orobanche aegyptiaca (Variety) was brushed onto a plate containing an extract
#prepared from the roots (roots) of either a bean or a cucumber plant. The number of
#seeds that germinated (y) was then recorded. Two varieties of Orobanche aegyptiaca
#namely O.a. 75 and O.a. 73 were used in the experiment.
#y: number of seeds germinated
#n: total number of seeds
#variety: 1, O.a. 75; 2: O.a. 73
#root: 1, bean; 2, cubumber
dat<-read.table(file="C:/jenn/teaching/stat579/data/seeds.txt",header=TRUE)
nrow(dat) #21 obsns
## [1] 21
head(dat)
## variety root y n
## 1 1 1 10 39
## 2 1 1 23 62
## 3 1 1 23 81
## 4 1 1 26 51
## 5 1 1 17 39
## 6 1 2 5 6
dat$variety<-as.factor(dat$variety)
dat$root<-as.factor(dat$root)
dat$resp<-cbind(dat$y,(dat$n-dat$y)) #dat$y #of success, dat$n-dat$y #of failure
head(dat)
## variety root y n resp.1 resp.2
## 1 1 1 10 39 10 29
## 2 1 1 23 62 23 39
## 3 1 1 23 81 23 58
## 4 1 1 26 51 26 25
## 5 1 1 17 39 17 22
## 6 1 2 5 6 5 1
#Descriptive Statistics
ftable(xtabs(~variety + root, data = dat))
## root 1 2
## variety
## 1 5 6
## 2 5 5
tapply(dat$y,dat$variety,mean)
## 1 2
## 27.27273 12.40000
tapply(dat$y,dat$root,mean)
## 1 2
## 14.80000 25.09091
##checking interactions
interaction.plot(dat$variety,dat$root,dat$y)
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABUAAAAPACAMAAADDuCPrAAAAaVBMVEUAAAAAADoAAGYAOjoAOpAAZmYAZpAAZrY6AAA6ZmY6kJA6kLY6kNtmAABmAGZmOpBmtv+QOgCQkNuQtpCQ29uQ2/+2ZgC2kDq2///bkDrbtmbb/7bb/9vb////tmb/25D//7b//9v////mR/DlAAAACXBIWXMAAB2HAAAdhwGP5fFlAAAgAElEQVR4nO3dbWPTynqoYUNXm7SFHtjd+OD2LAf4/z/yWG+2/BLHeiyNNJrr+rJCSBxYlm40Gmm0+QNAyGbuPwBArgQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGClh3QDcBoxk/U6K84orn/bwPrMnqjxn7BMU3wDwZQLAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABTfNjF/1/BYgR0DQ/VkJhhQQ01Q9WUFgdAU33oyUUVkZAE/5sCYV1EdCkP11CYU0ENPHPn/tPAIxHQJP/Ceb/MwDjEND0fwQJhZUQ0PQkFFZCQOegoLAKAjoPCYUVENCZGMdD/gR0NhIKuRPQGUko5E1AZ6WgkDMBnZmEwoR+f//8c8KXF9C5GcdDwK+vm7/+fue3/tF9tK93r08/gj/iHx9+iYDOT0JhsPcDut18aT54e9203kvtXcfXuUNAl0BBYaB3A3o46mzCd+pnqKDH17lHQJdBQmGQjwO6PexU36pzoLvqg8E/QECvLbdSxvEwxIcBPXzBpx/NJNI+cggqoNeW3CgJhcd9GNDDCP7wBXVAf3/fDJ+LF9Bryy6UgsKHDuXcbF76Ad325tr37UnPb/2AdraHz+wvvrQ3tu9/4vQ6dwnookgo3LdrZ4X+twvo7mym6BS+3hC+cwjo/3StrTtc6X7//BMCesvi82QcD/ccc/mvbUBP/ayH3L3wVb/z5SKg/1LPzL/0ctkV9OITAnpLBnGSUHhXlbnjgWb1QfWJL/3fOZ27/P397AjzTzPWb8f926a3dWVfbn7COdBrWaRJQeEdu66A1UWebUlfmt86fKYbfH85ffVZQ7fHfh6+uL0/qf3o6hMCekMmZZJQuOVwUNndl7m7vDTpcAh6GdDTSLyJ7Pb4W7tjd6sPv9z4hIDekEuXjOPhhi6Sf7qrlI7q8fp1QLtxfPOp7fGc5umj6oVebnxCQG/Ip0oSCld61y6dPuyOMt8J6OGT2/bKpW13/No7km1KfPUJAb0lpyZJKFzoHXZ2Ae3f8f5eQLuxez+gx7ml+oWuPiGgt+RVJAWFM9dHoN3VR59/3jwH2ga0HZU7An1SbkGSUOi5Pge6O86xXwa0Oat5O6D9U57NPP7VJwT0huxyZBwPJ9WUUNu5Zha+d+S4vxjCN/PqTUB3F0P4/qT79nIWfmsW/j0ZxkhC4eh48VJ73fwpoNUnzgLaXM9ZB7Q6T1p19xRQ14GGZJkiBYVWl8l66qhK6bbNZn3NfBe+prHb7lbOfXcD0imgZzcetb938YmH1sAT0BxIKDSO1yz9x9dNf/GQzTGgzbT8t7O729tw9gL60b3wvde5R0CzYBwPjbaYx9WYtl32/rM9P9pcOl8dTJ6a2HazF9CPVmPqv84dApoJCYVaHbbeeqDNkeKX023x7Recfu9YwX5AP1gP9Px13iWg2VBQGMxz4ceUd4IkFAYS0DFlHiDjeBhGQMeUfX4kFBZEQHMjobAYApofBYWFENAcSSgsgoBmyTgelkBAMyWhMD8BzZaCwtwENGMSCvMS0JwZx8OsBDRvEgozEtDcKSjMRkDzJ6EwEwFdAeN4mIeAroKEwhwEdCUkFNIT0NVQUEhNQFdEQiEtAV0T43hISkDXRUIhIQFdGwWFZAR0fSQUEhHQFTKOhzQEdJUkFFLIL6DbTedl+DeXUxUFhellFtDd5tzQhpYUFQmFqWUV0LfXzaXPPwe9QlFJMY6HieUU0N/fz4v56+vh13/9PeQlCguKhMKkcgro/iqXVVK/DXmJ4nIioTChnAK6vR6wHw5CB50GLTAmCgqTySigh8PNL1ef3A0bwxfZEgmFiWQU0MPR5vVwfT9sGqnMkhjHwzQEtAgSClPIKKCG8M9QUBhfRgE1ifQcCYWx5RTQ25cxXR+V3lF0QozjYWQ5BfTmhfTDbkUqPCASCqPKKaBNMc99+jHoFYrPh4LCiLIKaH8pJouJBEkojCazgP6xnN3TjONhLPkFdICrAb9y1PyPgHEIaJH8n4Ax5BrQekZ+2CVMFdnoSCg8L7OA7tsVQPex9ZQF9MTxOJxpFmwfdmFPVgFtDjsPf8PT0vRWpI+TUDg6XSQ5JCpZBbSbgP/r/7Z/ya0V6Z8iodDoX2Q+oKA5BfRw3Fn9zeoHy7Xd3FuR/jkKCpVdN6eyHXSFZE4B3bX/MmxP5ykOg3qLiTxHQqE+AG2PxXZDDkEzCuhxObvDkehx4G45u6cZx8P+dNg56ElrGQX0uKBy/7DTgsojkFBKt+tdFLldeUAPf0EBHZeCQmc74FKmvALa/hux7Q/hBXQUEgq1Q2gePy+YUUAPI/frv9fWJNJIjOPJwq0btB/y6A/YDbnFMaOAVn+xy1MTeyvSj0dCWb5wPx/dtPtz1I/8cWJ/jTuvOPYLHlVXGrxcfsKK9CNSUErXXm3+qJwC2twBfzoI3Z7/8hH68AEJpWj7gcdkWQW0+tv15sd2g/spoB8yjqdgQ/uZWUAPR529v95u+Jr02vAxCaVUgxfXyC6gT1KGR0goJapuQZr/KWuL3vN04TEKSnGqRTIXsET7onc8WXiUhFKWQz8HPiS9IqDcZhxPSX59jfRTQHmXhFKMQUsw9Qgo71NQCtE9Za21xtWYxqAHA0koJWgftyagH1CDoYzjKUD/iUgC+j4tGE5C4T0CyocUFG4TUB4goXCLgPII43i4QUB5jITCFQHlURIKFwSUxykonBFQhpBQ6BFQBjGOhxMBZSAJhY6AMpiCQkNACZBQqAgoEcbx8EdAiZJQEFDCFJTiCShxEsrK/Pr6+eeQrxdQnmAcz6r8/r4R0Dvs7WOTUFZkuxHQe+zr41NQVqJ+MpKA3mFXn4KEsgZvrxsBvc+OPgnjeLJXP1fu35wDvctuPhEJJXO76mGcJpHus5NPRkLJ2m7zYhb+I3bxCSkouRPQ++zhk5JQJrdpRX99l4DeZ/+elnE8E9tcBHHor+8T0Pvs3VOTUDImoPfZt6enoGRLQO+za6cgoWRKQO+zYydhHE+eBPQ+u3UiEkqOBPQ+O3UyCkp+BPQ++3RCEkpuBPQ+e3RKxvFkRkDvsz+nJaFkRUDvszenJqFkREDvsy+np6BkQ0DvsyvPQULJhIDeZ0eehXE8KyWgpCChrJKAkoaCskICSioSyuoIKMkYx7M2AkpCEsq6CChJKShrIqAkJqGsh4CSmnE8qyGgpCehrISAMgcJZRUElHkoKCsgoMxFQsmegDIb43hyJ6DMSELJm4AyKwUlZwLKzCSUfAkoczOOJ1sCyvwklEwJKEugoGRJQFkGCSVDAspCGMeTHwFlMSSU3AgoCyKh5EVAWRQFJScCysJIKPkQUJbGOJ5sCCjLI6FkQkBZIgUlCykC+uv//D32D4myV2ZDQslAkoB+3Wxexv4xMfbJfBjHs3ypAnrwZeyfFGCPzImEsnRpzoHuNo1vY/+woeyPeVFQli3ZJNK2Sejnn2P/vEHsjrmRUJYs4Sz87+9NQ/+acUrJzpgd43gWLO1lTO3Z0PmmlOyKGZJQFiv5daDtceinH2P/3IfYEbOkoCxU4oB2s0lzHYXaDzMloSxSyoDujvNI9WHoHJc12QtzZRzPEiUL6MWVTL++PjUh35xMHT4dZR/Ml4SyPGmvA+0ddG4HB/QQza6Y2+h5AHtgziSUpUl4J9LZ8eJhED/0+PEU0GM/Bx+E2v/ypqAsS6qAjjDrfgzovjuWHX6Tvd0vdxLKkqQJ6Ci3cHYB7Z8+3Q0ss50ve8bxLEhG64F2Ad33zqX+/j7sENSutwISymKkDOj+ucs/u4CeHXXuhp0FteOtgoKyENMH9O21SVx3K3x4SaYuoNt+QPfDJvPtdyshoSzC1AGtstkE9Dh1Hp1POgW010wBLZRxPEswcUDrw846cW+vzfB9G16O6eYQXkCLJaHMb+KA7o9D9l3vSDQ2iO9NIp1eYOscaLkUlLlNHNDj+crDoWj70aGDsbvgu4D2Z973A2+pt8Oti4Qyr2kDerp56Ow+zNhM/PHRSsdqVicIht0QandbGeN4ZjVpQKvznu8InAftVmPedKdVt8Pn9O1sqyOhzCijgP7pNbQO6G74NVF2tRWSUGYz+RD+ONg+jrWHr8N04dDldjp+8LkAO9oqKSgzmTagx27uT7cgHaKaajH6W4e+iX40aXlrmcXEs/C7ZrR9OGbsrt2sJn5SrUUvoOXw5jKHiQN6mvhpjzr3m1mfDW8fWy8JJb2pb+XsCtpEs5lVemJxu/oFehNQ/XOrD/3h7GErpqCkNv1iIvXjPNpBe53TJ/q5PT+cFVDOSShppV0PNHwXUuP0KI/uIFRAOWMcT1IZLahcj9+rw9fdqaACygUJJaGcAtq/sX5zXGRUQDmnoCST5plIZ0uADlxA6ai/iEi3KJ6AcoOEksgsAY1dx/Tr69kydnVBBZRbjONJI31A+w/VHOQsoFVBXwSU90goKUx/J9KYK4n0Atrc0SSgvEdCmV6yO5H6YpcyXTzCuCroNwHlfQrK1KYewu9H62d9ONu/CL+6wf6fAsodEsq0ZphECquuA+3Hd9+7SfRBdqfCGMczqZwC2pxR7Q3j9wLKRySUCeV0IX1T0P4xaHVMKqDcp6BMZo6A/vqv8Hp2v79fnEHdCSgfklAmkiig53NJ8y0Iaj8qk3E800gT0O1GQJmVhDKFJAG9fDrnEyuCPsk+VC4FZXxJAtrcdtk8RHM7Zz8FtGgSythSBPT39/rmzX33NHfPRGIexvGMLM11oPXU+dtrfTnoIaeG8MxEQhlVooB+O/2nHcrPw85TPAVlRAkD2l3D+fYaW41pDPYdJJTxJAzon21z6Pn26jIm5mQcz1gSTSLVh567ZiHQvetAmZmEMo4klzG1E+/75qFwu+CCymOw19CQUMaQ6kL6Kp2H/xzSud+YRGIBFJTnpbuV88vpjk6XMbEEEsqz0gS0fvzG8QEf8x2ACih9xvE8KdVydvtmTeXtE0/0GIPdhTMSylPyWlD5aXYWLigoTxBQSiehhAkoxTOOJ0pAQUIJmjag7bT7JXcisTQKSoSAQk1CGU5AoWEcz2BpzoHuN81t8E1SPdKDZZJQBkp1L/zp6vl919I52D24S0IZJNFD5V7e/VVadg4+oKAMkHJB5Zb1QFk0CeVhAgoXjON5VKIV6c+H8BZUZtkklMckOgfaOwTdWVCZ5VNQHpFqFr4btlcfmoUnAxLKx9JcB7o7u47edaDkwDieDyVaTGR/yueMx58CyiASygeSrcbU3tU53wR8ze7AIArKXZazg3sklDsEFO4yjud9AgofkFDeI6DwIQnlNgGFBygotwgoPERCuSag8BjjeK4IKDxKQrkgoPA4BeWMgMIQEkqPgMIgxvGcCCgMJKF0BBQGU1AaAgoBEkpFQCHCOJ4/AgpREoqAQpiCFk9AIU5CCzdtQH9/n/EZ8LfY2hmXcXzZpg3or6/VM+B/fZ31QXJ9tnXGJqElE1B4koSWa/Ih/KGdAsrKKWipJp5E2m5ume/ZxjZzpiGhZZo4oG+vAkoRjOOLNPVlTDcLKqCskIQWKMV1oM6BUgYFLY6AwngktDDuRIIRGceXRUBhVBJakmQB/fV15vmjmi2b6SloORIFtM3n7Am1YZOChJYiTUB7/dxs5pxRslmThHF8IZIE9Pf3w9b0rf5wf/hoxgWabNQkIqFFSBLQfS+aVUy/jf0zH2aTJhkJLUCSgG77Jz6bFZpmYoMmIQVdvRQBPRx0fun9cjfjGN72TFISunJp7kQ6G7Tv3QtPMYzj101AYVISumaG8DAxBV0vk0gwOQldK5cxwfSM41fKhfSQgoSu0gy3cs55M7xNmNko6ApZTARSkdDVsZwdJGMcvzYWVIaEJHRdBBSSktA1EVBITEHXQ0AhOQldCwGF9IzjV0JAYQ4Sugr5BXR7vKI0cEu9TZbFUNAVyCygu825oQ21xbIgEpq9rAL69rq5NPDCfNsrS2Icn7ucAlqvSXK+MN7QlUlsrSyLhOYtp4Dur3I5eG082ypLo6A5yymg2+sB+9DVmW2qLI+E5itRQM/PXsZWFLl4Mkhj4PNBbKgskHF8ttIE9GLyPBbQi2fTNQY+oc5myiJJaKZSPdJDQOEeCc1SoofKjbEOqCE8q6agGUr0WONR1lE2icS6SWh2UgT0ULnrQ8eA25cxDXpp2ydLZhyfm5wCevNC+mHHtrZOlk1C85JoCD9KQM8f7tn49GPYH862ycIpaE6STCINnOm5Y3vRT4uJsD4Smo8kAR1rDF+znB1rZxyfjTQX0r+9jnYMOsTVgN9mSR5sq5lINIk0xoX0gT+KgJIrG2sWVh3Qa7ZJsiGhGcgsoPWiJL2zAUOv0bdFkg9DpuXLaTm70wTScfpIQFkzCV26rAJ6moDvDkIFlHVT0GXLKaDV+L1aj2l3KmgmAXUgQZhtZ8nmCOiv/4qdA9129x1V93Q2Bc0poHYEQmw5C5YooOcrgoZXpD+e+9y2Bc0koM2PthsQY9tZrDQBvbgDc4QFlduCZhRQiJPQhUoS0MvnuQ96kObR+Yr023ouXkAphIIuUqoV6V+quZ+X5uNYPy8C2iwFmmlAHU8wnI1mgRItZ1cNt9vHF+2i19H3z4H+aZ8Jn3NA7Q8MY5NZnoQLKr+91pPoh+YFD0F35wevb6+bT//MM6AVewOD2WiWJlFAv53+0w7lA6pTqf1l8fbDZ6RsfeRNQZclYUC7henfXqNL2+3ObuNsCyqgFEVClyRhQP9sm/i9vYYXE6kK2j8GrY5J1xBQIzMeZ2tZkJTPRGqf7LF/YjWmq8crDZySWuiGZ1KJIWwri5HomUh15fbNrZjjPSFpuAVvdnYKHmdjWYhUF9JX6Wwe7LEPPcxoJLY6VkJCFyHdrZxfTnd0Bi9jGoFtjrUwZFmCNAGtL3o/Lk0/3wFoLgG1b/AAm8n8Ui1nt29WotteTKOnlsn2ZlKJh9hI5pbTgsojyGhrs2/wAFvJvAQUsiahc0oW0Or852EYv5vvEqaKTY3VMVaZUaKANtNHVUDnnIPPNqB2Ee6xfcwmTUDb6fdDQLezXsWUdUDtJLzLxjGTJAFtHgLXLGe3jT7RYxQZb2UCyj02j1kkCei+vnapXQ90N+eFTLYx1sq/sHNI9EiPaurotKCye+GfZFfhBglNL+FqTG1A48/0GMFKti8nRLnJZpFawvVAu4A+s5zds9azeWkot9gq0hLQbNlVuME/rEnNMITfOgcK05HQhBJNIlXHnG1AD8ej1gMdmT2GPglNJtVlTC9dQNuV7Way0s3KCVHO2R4SSbag8uefdUCrx8LNeDv8ercqDeWMzSGJpLdyNma8EWnFAbXHcM6/qCkkWkykGrgvoJ+rDiick9DpJVvOrk3orPksJqB2HGo2hKlZUHmNnBClZUOYloCulIZSsx1MSkBXy45DTUInlCig+02fWzkhIQWdTLIFlQV0Rnag0tkCJpIkoLuNgM7KCdHi2QKmkWgxkRlvfz9T7jakoaWzAUwh0XJ28z7M+KTkLcj+UzpbwPgSBXS+pyCdswFRMgkdW8L1QJfA5lMzmiuVd35kiSaRFnIKVEAbx+m8uf8gJOd9H1WSgLbP9FgAm86RHalU3vkRpbmQ/u11IYN4Gw44FTqeRHciHQrqOlBYCgkdiTuRMKgrkLd8HO5EwqRSkbzhY3AnEg37U3G8489LdCH9QvopoNAjoc9yJxKUy7jjSe5E4prdqhze66ckmkSymEhWTCqVxBv9hFSXMbkTKTcCWg5vdViaC+kXs6Cd7QSu+dcyKtEkkutAYckkNEZA+Zi9qwDe5AgB5WMmlYrgHR7Oc+F5jIAWwHs8lIACHf9MDiSgwImEDiKgDGcnWzNv7gACynAmldbNW/swASVGQFfMm/soAQWuSOhjBJTn2dlWyJv6CAHleU6IrpI39WMCyig0dIW8px8SUEZiZ1shCf2AgALvk9C7BJTx2efWxLt5h4AyPidE18W7+S4BZRIauibezPcIKBOxz62JhN4moMADFPQWAWV69r1V8DZeE1Cm54ToOngbrwgoSWjoKngXLwgoidj1VsHbeEZAgSEktEdASc8umDXj+BMBJT0nRDPn/esIKLPQ0Lx5+xoCykzsgXnz/lUEFAiRUAFlCYwH8+R9E1AWwHPmc1X8uyagLEPxu2KmCn/bBBR4RtEJFVDgKSUPHgSU5Sl5j8xSuW+YgLI8JpWyU+q7JaAsk4Bmpsz3S0CBMRT5T16mAf31tXqz/vp76PeV9wZDMgUmNKOAHqLZFXPbnSN7Gfgapb29K1Hgjpmp4t6oLAN67Ofgg9DC3t21MKmUj8LephwDuj+8RV/azww8Bi3qvV0XAc1EWW9UhgE9/Pfzz/Zzu82nH0Neo6B3FmZSUkIzDOi+Pf6s/P4+7BC0mPcVZlROQjMM6NlR527YWdBC3tW1K2f/zFYp71CGAd32A7o/DecfUcabunomlTJQxvuTZUB7zRTQUgno0hXxDmUY0J0jUMhBAQnNMKD7zebb8ZNb50BhsVZf0AwD2p9578/IP2Ll72apCjjSydbK35q8ArqpL6E/VvOQ0mEjeAFdJ5NKC7buNya/gNbqbNZ3dH778Pv6VvxOlm7d+2nW1vzWZBTQP72G1gHdDe6ngMIM1lvQvALaeHvtpuOHLsYkoGVY7/6arbW+JTkG9GGbG+b+M5GAt3p5VvqWCCir5N1enFW+H6sO6LU1voXcJqBMT0ABggSU9XMwykQElPVzQpSJCChF0FCmkFFA+3ci9ViNiYcIKOMTUICgjAJa3YEkoIzAwSgjySmg9THooNXrrthv+OOEKKPJKqBVQYc9xviSvYaGhjKGvAJ6XEckyi5DR0B5XmYBrZawG7iC3Rm7DDCe3ALaPdcjSEC5wXCeoNwC+uQhqN2EG6zVRVB2AX2OXYR3CCgBAgoQJKAAQQIKlwzneZCAwiWTSjxIQOEWAeUBAgoQJKAAQQIKHzGc5x0CCh8xqcQ7BBQeIaDcIKAAQQIKECSgMJThPC0BhaFMKtESUIgQUP4IKECYgAIECSg8y3C+WAIKzzKpVCwBhTEIaJEEFCBIQGFsDkaLIaAwNidEiyGgMAENLYOAwiQEtAQCChAkoDA1B6OrJaAwNSdEV0tAIQENXScBhSQEdI0EFCBIQCE1B6OrIaCQmhOiqyGgMAMNXQcBhVkI6BoIKECQgMLcHIxmS0Bhbk6IZktAYQE0NE8CCosgoDkSUIAgAYWlMZzPhoDC0njOfDYEFJZIQLMgoABBAgoQJKCwdIbziyWgsHQmlRZLQCEHArpIAgoQJKAAQQIKuTGcXwwBhdyYVFoMAYUcCegiCChAkIACBAko5M5wfjYCCrkzqTQbAYU1ENBZCChAkIACBAkorI3hfDICCmtjUikZAYU1EtAkBBQgSEBh7RyMTkZAYe2cEJ2MgEIBNHQaAgpFENApCChAkIBCaRyMjkZAoTROiI5GQKFAGjoOAYUiCegYBBQgSEChdA5GwwQUSueEaJiAAhoaJKDAH+P4GAEFCBJQ4JyD0YcJKHDOCdGHCShwRUMfI6DADQL6CAEFCMovoNvjM1tfhn+zgMJghvPvyiygu825oQ21GcBgnjP/rqwC+va6ufT556BXsAlAiIDelFNAf38/L+avr4df//X3kJewCQDjySmg+6tcVkn9NuQlBBQYT04B3V4P2A8HoYNOgwooPM1w/iijgB4ON79cfXI3bAzvbYenmVQ6yiigh6PN6+H6ftg0krccRiGgNQEFCMoooIbwwLJkFFCTSLBIBQ/ncwro7cuYro9K7yj1bYYJFTyplFNAb15IP+xWpBLfYkhAQMd6xbFf8KQu5rlPPwa9QolvMTCVrALaX4rJYiLA3DIL6B/L2cHSFTSczy+gA1wN+It5W2FGBe1tAgqMr5C9bdUBvVbCWwqkIqAAQQIKTGvFw3kBBaa14hkIAQWmJ6CPvuLYL9i5cR/SZui9nGt8C4G5CChAUEYBvfVUYwEF5pNTQOtj0EGr110RUGA8WQW0KujA5ZcuCCgwnrwCWo3iBz3C45KAAuPJLKB/dpsbT5Z7nIAC48ktoIdB/DOHoAIKjCe3gD55CCqgwHiyC+hzBBQYj4ACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIEFRdQgPGM3qixX3BMc//PBtZl9EaN/YIApRBQgCABBcw3lg0AAAdGSURBVAgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAb326+vnn3P/GWBOb6+bg08/5v5zLJ6AXvn9fSOglOzX103LnvABAb2ytdlQtFM/7QofEdALh+NPWw1F2x12gS/VB4djic3L3H+aZRPQc825HwGlXNUB6Lfmw5194QMC2lePXf7NOVBKtj8ddlbjsW+z/mGWTkD7qrHLN5NIFG3XDuArWwG9T0D7dtU/vQIKra1Lme4T0CsCCo1fXzd//T33H2LRBPSKgEKjP5rnFgG9IqBQe3t1APoBAb0ioFA59NOe8AEBvSKg8Ke+nMmO8BEBvSKgoJ+PEdArAgrVFaDOf35MQK8IKMWrbkFyF/wDBPSKgFK6akkI1y89QkCvCCiFO/TTDUiPEdArAkrZfn3Vz0cJ6BUBpWiWYBpAQK8IKEXbb85o6T0CekVAKVn9TAYBfZCAXhFQStZ/IpKAfkRAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggSUbPz+/vnn3H8G6BNQ5vfr6+avv9/5rX90H+03lU8/Br/67++b98J7enWIEFDm935At5svzQdvr5vWe6l91/sBPb46xAgo83s3oIejziZxp34OL+i7AT2+OgQJKPP7OKDbQzi/VedAd9UHI/1YAeVZAsr8Pgzo4Qs+/WgmkfaBQfw7BJRnCSjz+zCghxH84QvqgN6ZEhpKQHmWgDKjQzk3m5d+QLe9ufZ9e9LzWz+g3fcdK3r8eHs2Tb89fMe+/nUvufv29c5evdfRw2tJKgMIKPPZtbNC/9sFdHc2U3QKaG8If/zO7kxo27/zb60D+j+b84AeZ6LqXx5fvZfvfeQyKQomoMzm2Lx/bRN2imAdxVNA69/50gvooYUv7YfbOnoX31p9+l/qXr4cB/29mfzq16dX33bZPHzlWOdXKYOAMpdq/H480Kw+qD7xpf87p7OUh7QdDx3bX7fRa44fr791ezwWbQNavUT9Pfs6q3/OzrAeL5YygmcIAWUuu65w1aFhW9L2sPLwmW6Y/eX01f2G7rrf2XcHqxffuj3O1rcBPU3fNycETq9+PPDcGcEzjIAyk9NB5CmlnW5i6GyevBtzv7Rf0Z3rPI9e963b47e2Ae19Yfuyx1dvw2kEz1ACykx6M+nNJPtRPV6/Dmg3ju/OcdbRO7zKy59b37o9TjM1Ae1fKtWeQe3f5/TljxE8wwkoM+kV7fRhd5T5TkAPn9xuzsbfp+n4i289HXAeA9rXnjI4nmGtfr0b7QpTSiGgzKR32NkFtH/H+3sB7cbmzfdcX6L0XkD7X3AZ0ObypcPX9Q9m4WMCykyuj0C7g8TPP2+eA20D2l3BVBey+8X1t94I6OUJztOr1xfQH75irLvsKYWAMpPrc6C74xz7ZUCbE5rnAa3PWHYj+OtvvTGEvxyg9/Jc3bdkBM9gAspMqvme9oivmYXvTcvvL4bwu3ruvQnornfx0f9r581vfOtlQHs/rtML6GEM/08jeAYTUOZyvHipvfj97OL484Aejja7Wzmrc5nf2u//9N+vx5hefutlQHsHqVeXMdXf9u9G8AwmoMyla109vVOldNsmrr5mvptq7672bG/l3B9vMGq+r63k9bdeBbQ+Tfqt/eLzV2+/0QieoQSU2RwvPPqPr2eLh2yOFWymzr/9ObsI6Xg9/OluzRvfehXQs6/pjmu7ptYfG8EzlIAynzZpx9WYtl3f/rM9YdlcOt9etXTZz/rbT7NAF996HdDelUz9G+yPH7uNk8EElBnVCeutB9okrndve/sFp9/rX9Z0NrF++a03AtrdUH880uy9ej0PP91flJUSULIx5XPh3cZJhICSjSkDaiEmIgSUbEwY0N4CzfA4AaV4/SclwRACSvH2ZzNL8DgBpXjVBL4ZJCIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCDo/wOXvXs4bB/+3wAAAABJRU5ErkJggg==)
#fit logit model
#y_i ~Bin(n_i ; p_i )
#logit(pi) ~ variety+root+variety*root
fit1<-glm(resp~variety*root,family=binomial(link=logit),data=dat)
fit1
##
## Call: glm(formula = resp ~ variety * root, family = binomial(link = logit),
## data = dat)
##
## Coefficients:
## (Intercept) variety2 root2 variety2:root2
## -0.5582 0.1459 1.3182 -0.7781
##
## Degrees of Freedom: 20 Total (i.e. Null); 17 Residual
## Null Deviance: 98.72
## Residual Deviance: 33.28 AIC: 117.9
pval<-1-pchisq(33.28,17)
pval #0.01 reject sufficiency of the model
## [1] 0.01038509
#residual plots
rp<- residuals(fit1,type="pearson") #pearson resids
resDev<-residuals(fit1,type='deviance') # Deviance residuals
plot(resDev, ylab="Deviance residuals")
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABUAAAAPACAMAAADDuCPrAAAAWlBMVEUAAAAAADoAAGYAOmYAOpAAZmYAZpAAZrY6AAA6OgA6kJA6kNtmAABmAGZmOgBmtv+QOgCQ29uQ2/+2ZgC225C2///bkDrb/7bb////tmb/25D//7b//9v///8EvmTIAAAACXBIWXMAAB2HAAAdhwGP5fFlAAAgAElEQVR4nO3da2Pa2HqAUTKdNm6d1jlDhw5O/P//ZiUQIBzbMS9b+6a1PpwkHocjGfFEl62tzQsAIZvSCwDQKgEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUICgugO6AUgmfaKSv2JCpX/aQF+SNyr1C6a0wD8YwGoJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgLZkqRm0gBABbchycxACEQLajnM4FRTqIKDNmGWz5dWAjghoM+bL3vJ6QD8EtBkCCrUR0GYIKNRGQJshoFAbAW2GgEJtBLQZAgq1EdBmGMYEtRHQdhhID5UR0Ia4lRPqIqAt0U+oioACBAkoQJCAAgQ1FtDnh81m8+c/5z///L754+8b/r6AAum0FdDtdBHl6+kLAgqU01RAt+dhPKedUAEFymkpoOPx+9Pw6+5SUAEFymkpoNvNl78OvxmyORVUQIFyGgroEMvzuc/tVFABBcppKKA/vh0O4I+mgn4c0M0bllo6YH1aDehY0K8CCpTUbEDH86CPDuGBghoK6Pwc6PGPmycBBcppKKDj8KXZLug4qunLvwQUKKalgI7jQB9nf94fTmoKKFBISwE9jKCfH8bvBRQoqKmAHgo63wcd90kFFCikrYCOV44er76wE1CglMYCei8BBdIRUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQCtlMSi8HxAkoZWw2CkrzBJQizuFU0GX5Z2pRAkoJs0+092RJdvSXJaCUMH8jvCnLsaO/MAGlBAHNwo7+0gSUEgQ0Cz/mpQkoJfhkZ+HHvDQBpQSf7Cz8mJcmoJTgk52FH/PSBJQSXN3IQkCXJqAUYXxNDgK6NAGlDCO8M7CjvzQBpRD9zMCO/sIEFDpmR39ZAgo9089FCShAkIACBAkoQJCAAgQJKECQgAIECSjQr4WHcQko0K2lbyQQUKBXi9/KKqBAp5afTEVAgU4tP52fgAKdEtDEBBTWQ0ATE1BYDwFNTEBhPQQ0MQGF9RDQxAQU1sMwpsQEFFbEQPq0BBTWxK2cSQkorIrJRFISUCAdAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIyhvQ54fN5s9/Uv8/3kBAgXQyBfTHt7GbYz83my9/pf6//DwBBdLJE9D9ZvPH3y8/v28Oht+WIqBAOlkCOu55Dvudwy9DO3ebzWPq/89PE1AgnSwBHZr59fjLmM7t4Q+3+/Ft85ab9mcFFEgnR0CHQ/fxDOjwy+H05z54DC+gQF1yBHQo39Pxl0PsogGdrkHdENC3vj30/wzwhowB3R8P5OMBPeyD3nL+VECBJWUM6HZz+OVlFx8KOrzSfYOgBBRIJ9M50MfLKdAhgrGLSKPhKP6ugfgCCqST5Sr8dkzn7ngT0ngi8yn+6ru7/vYaA+rcBSwm2zjQ0dPhOP6ufchh//Wev76+jjj7C8vJcyfS7vAJHsu3u/dGpPt2QVeXkXM4FRTSy3Qv/H5zvgR/zxH43dZWkVk217bqkIHp7Lo2X9+1rTssT0C7JqCwJAHtmoDCkpYNaIrb19Mu3MoiIqCwJAHtmoDCkgS0awIKS3IOtGuxYUyG3sPnCGjfIgPp3bwEnySgnbu9hm5egs8qEdAf/+0caD7hfq7ypwU3yXgrp4tITXDdKcJJj5XKE9Ctq/DNENAAp43XKut0difl5hOxgf+WgN7OaePVyjSh8jgV0+4wH9O2ZD8V4fcE9GZOG69XxscaTw+T2xU8grd9/56A3syPbL0yPVRufJbm88PhmUhDTh3CV0wNbuZHtl55nws/PZUz/lC5e9m8f8vx6M0EdL0yBvT4cM5xT/SuB2vexeb9e66I3EpA1ytjQF+2x13P5wfDmKpmTM6NBHS98j0Xfjx2P+x67o0DrZx+3kZA1yvLMKbpwvt+c7iKtLvrwcT3sXmTntPG65VrIP2YzuGXIZ37jYtI9MVp49XKdyvn4+WOTsOY6IvTxmuVJ6A/vx+qOU1QX24HVEBZhn6uVK7p7PaH05+nXdFibOFAOiZUBggSUHjF8TifJaBwzRUhPi3TnUgmVKYVxiTxeQIKc0bFcwMBhTn3ZXKD3OdAp8mVS/GJ4DcElBvkv4hU8lZ4nwh+R0C5Qf6AmpG+R/1ctxZQblBgGJPZmPrT0cgfAeUGRQLqIlJnehr5I6DcQEC5W1cjf7paGZaWP6A/vjmE70xfO2097U6ztOwBHQeFmlC5L30FtKcTuiytxED648x2RfhMLKGzgHY0pICllQhowQlBfSiW0FtA4bPyB7TcFaQXH+9lCChrZTo77iagrJWAcjcjf1grAeV+Rv6wUgJKAkb+sE4CSgr6ySotG9BfplIufSHeJxxIR0ABggQUICjPOdD9+f7NManl5lMWUCChLAF9fpjdvrl3LzzQhywB3V7Nv7Q1GxPQhUz3ws+P2vfOgQJdEFCAoBwB/fn91SG8GemBHmQ6BzrbBd2ZkR7oQ66r8KfD9vG3rsIDXcgzDnR3NY7eOFCgC5kmE9lX8UQkAQVSyjYb03RXZ9EHeggokJLp7ACCBBQgSEABgpafzu7p10nt3IkE9EBAAYIEFCDIOVCAIAEFCBJQgKC8AR2nEik3l92LgAIpZQroj29jN8d+lr0ZXkCBdLI9lfOPv8eJlUvfDy+gQDrZ5gMd9juHX4Z27jazR3TmJqBAOlkCOk1CP6XTUzmBPmR6JtJ4BnT45XD600PlgD5kfCrn8MuhnAIK9CFjQPfT0+QEFOhDxoCens2581hjoAuZzoE+Xk6BDjl1EQnoQabnwg/p3B1vQhqHNJV7LKeAAulkGwc6Pc94OI4veTOngALp5Hwu/BjOXdkHcwookE7G58JPl+DLHb+/CCiQkunsAIIEFCAoW0DHByON1+KLTgcqoEBC2eYDnWYC3ZUcxCSgQEp5Ajo9l3MI6HZTtKACCqSTJaDjTMp//vP8MN6JtDWhMtCJLAHdHyYCPQbUhMpALzLdynm8ifMQ0Gl20DIEFEgn32Qip4AOu6CmswN6kHE6u1NAzQcK9EFAAYIKHMJvnQMFupDpItK4zzkF1ITKQCdyDWP6egroOCbUhMpAD/LciXQYPX8I6M6EykAvst7KeWRCZaAPmSYTGQ/cK+ingAIJZZvObkpo0XwKKJCSCZUBgvIMYyq833khoEA6mQbSl52H/kJAgXQy3cpZbuj8NQEF0rEHmuqFjxZ6eaBGue5EKvo0+IulCncZ5brM6wM1ynMV/vmhkn3QhQJ3DqeCwppkOgd6rbfp7GbZFFBYEQFN/KoKCushoIlfVUBhPdyJlPhVBRTWQ0ATv6qAwmf0MW5FQBO/ausbBGTRycg/AU38qo1vD5BFLyP/BDTJqxrGBDfo5iMjoGleto9/TiGPbg7aBDTR63ZxQgfyEND3XzH1C6ZkMhGogIC+/4qpX/CV7XlnLzBHXtPvFfRCQN9/xdQvOLd7dU/TrQ1t+r2CXgjo+6+Y+gUvnh82r914V2jT7xX0QkDff8XUL3h2eLDnrJiHe+xvmyav6fcKemEY0/uv+PaXx9p9+etld8e8oPtfcjkm9aapmtt+s6AXvYz8yxTQ44RMY0DvmJx+++sB+62PW2r7zYJudDLyL09ApwnthoCO19CDBR12Nx9/+eLutmP4xt8t6EYX/cwT0PFI+89/nh+GgI4FDU4HOlT41/TuP3q1Xy45Nf92ATXJ9VC5x/Ea+hjQ8SD+1x3JzxBQoC5ZAro9HGdPAQ0/5NghPFCXTM+FH8s3BXRoXvAY3kUkoCqZnok0HnufAvrhUfdH3h7GdNP5AAEF0mkpoG8OpL/ttQQUSKfAIfw2eA70jcd7HkZG3bRwAgokk+ki0rifOAX01tOWr17JZCJALXINY/p6CujNd1++Zjo7oBJ57kQ6jJ4/BHR38/wfSQkokE7WWzlDM9AlJaBAOpkmEzlcQC/fTwEFEso2nd2U0KL5FFAgpZYmVE5AQIF0BBQgKFdAd8exSz++lT2GF1AgnXxX4Q8BHR8LFx9Gfz8BBdLJNqHycdKPfXAAfCoCCqSTJaC72T3ru/gzPRIQUCCdTJOJzM58Dkfx5XZBBRRIJ9N0dvNJO2+cRD4pAQXSyTgf6El4PtAEBBRIR0ABgjKdA50ftMcnVL6fgALp5LoKfzkJuo8+1jgFAQXSyRLQcRz9NI5p/G3Bm5EEFEgnz51I+6vncJQbBiqgQEKZ7oUfb+GMPQYuLQEF0sk2G9M0KX3RfAookJLp7ACCBBQgSEABgvJfRCr6YCQBBdLJE9Dd1TAmAQW6kCWg18NABRToQ5aAbos/zvhEQIF08k+oXJSAAukUmFC5JAEF0hFQgKBMh/ACCvQn03yg5aZQviagQDq55gOtZBdUQIF08gykf36oZB9UQKFZp4HkpZdjJtNFJAPpgftcClJ6SS4EFGjBOZw1FVRAgQbMslnRx9h0dkAD5p/dej7HAgo0QEArUM8PHriFgJ78+G/nQIGbrDqgexeRgDusOaBbV+GBe6w4oNdPRNpsnlL/f35aPT944BYrHsY07IB+HacU+Xr8fbl+1vSTB26x2oH0P78f7oTfHw/ddyWnp6/oJw/cZMW3co6TMT0/fPnr5ZBTh/DAzerrZ66APl1+mQ7ly6jqZw80LmNATxPTPz+Um9pOQIF0Mgb0ZXvc9Xx+MIwJ6EHOZyJNT/bYGwcKdCHTM5EOydxvDleRSj4hSUCBdHINpB/TeXywx37jIhLQhXy3cj5e7ug0jAnoQZ6A/vx+qOY0NX25HVABBRLKNZ3d/nD687QrWoyAAumYUBkgSEABggQUIGjZgI5XjZ481hjok4ACzNwy6ZOAAlzcNO2oc6AAZ7dNfC+gACc3PnpJQAFObnz4Z575QEvefHRFQIEPVBnQohOIzAko8IFaA7o5TgZamIACH6gwoIeZQAuPX5oIKPCBKgP6Mk5Hf1RuNvqRgAIfqDWgL+OzPDbmAwUqVvcwpp07kYCK1T6Qfi+gQLVqvpVzaw8UqFo9k4lcOz1TruCwegEF0skW0AquIL0IKJBSnoCe6ll2DNOLgAIpZbwTqfgo+hcBBVLKFdAa7uN8EVAgpTwBrWMmkRcBBVIyHyhAULaATsfxO/fCA73IFNDjdaQxoGVnBhVQIJ08Ab3MCLotO7eygALpZAnoz+/jENDnh/FS/LboeCYBBdLJEtD94fbNY0DHg/hy93IKKLW45Y5rapUloNvDLUhTQIfd0XIXkmyuVOKmOX+oVY6ADskc9zmngA67oGZjYu1um3WSWmW6E2m8bnQK6F5AWbsb5z2nVgIK+d345B1qVeAQfuscKGsnoJ3IdBFp3OecAjrsj5abE9S2ShUEtBO5hjF9PQV0HBNabiS9bZUqCGgn8tyJdBg9fwjoruy0yrZVqiCgnch6K+em+MTKtlWqIKCdyDSZyHjgXkE/bavUwTCmTmSbzm5KaOHnethYm9bRvTsG0vfBhMo0o6u7H7tamfUSUFrR2U6bfvYg23R2FTwU/kVAW+a0IfXJENDdZq7s8+V88trlwjX1WTygs+vvR0UfiuSD1y4BpT5LB/TYz9Nu5650QX3w2iWg1GfpgL5+BtK+7JlQH7x2CSj1WTigzw+b4xRMH3wlJx+8dgko9Vk4oG88AMkzkQgRUOqzbEB/fv/11iPT2RFiGBP1WTygv1wy8lA5YjobSE8Plg3om3ubW4/0IMTdj9RGQGmHflIZAQUIElCAIAEFCBJQgCABBQhaPKBvEVCgBwIKECSgAEGeiQQQJKAAQQ0FNMX5AAEF0hFQgKCGAjpOZi+gQD1aCuhhH/S+2ewFFEinqYCOBb3vgUoCCqTTVkDHo/i7prMXUCCdxgI6PpLu6fff9S4BBdJpLaDDQfznd0Hfuua05MIB69JaQG/aBRVQYEnNBfQ+AgqkI6AAQQIKECSgAEHNBnQXGlIvoEA6AgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQc0GNEZAgXQEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAoRWbSenl4ExAoRGbjYLWRkChDedwKmg9BBSaMMumzbgaAgpNmG+7tuNaCCg0QUBrJKDQBAGtkYBCEwS0RgIKTRDQGgkoNEFAaySg0ATDmGokoNAGA+krJKDQCLdy1kdAoRX6WR0BBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABffO7zNkA/J6AvvVNCgp8goC+8T3mrQU+Q0B//RZPTgA+RUA/+hYFBT4goB99i4ACHxDQj75FQIEPCOhH3yKgwAcE9KNvEVDgAwL60bcIKPABAf31WwxjAj5FQN/4HgPpgc8Q0Le+ya2cwCcI6JvflaOfKg2tE9BS7OdC8wS0EGdaoX0CWoZr/dABAS3DaFPogICWIaDQAQEtQ0ChAwJahoBCBwS0DAGFDghoGQIKHRDQMgxjgg4IaCEG0kP7BLQUt3JC8wS0GP2E1gkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIDSt82k9HLQJQGla5uNgrIcAaVn53AqKEsQUDo2y6a3ngUIKB2bv9/ee9ITUDomoCxLQOmYgLIsAaVjAsqyBJSOCSjLElA6JqAsS0DpmGFMLEtA6ZmB9CxKQOmaWzlZkoDSN/1kQQIKECSgAEGNBfT5YTgY+/Of859/ft/88fcNf19AgXTaCuh2Op/19fQFAQXKaSqg2/MV1dNOqIAC5bQU0PH4/Wn4dXcpqIBCBVY71qGlgG43X/46/GbI5lRQAYXy1jvatqGADrE8n/vcTgUVUChuxfd7NRTQH98OB/BHU0EFFEpb84wDrQZ0LOhXAYX0bj0cX/OcV80GdDwP+iigkNrNJzQFNOkrpn7Bk/k50OMfN08CCmndfkJTQJO+YuoXPNtt5rug46imL/8SUEgpcEJTQJO+YuoXPBvHgT7O/rw/HGZ8ENDNGxZbOuhBoIYCmvQVU7/gxe7qNs6poAIK6QjoTZoK6KGg833QcZ/UITykEwqoYUwJXzH1C879/H4V0DGpAgrpRHYnDaRP+YqpXzCl1b2/cJvQ8fh6T5EJKHARO6G51n62G9DdaWaRm6zwHYZbrPmEZoCAAjMrPqEZIKDA3HpPaAYIKHBFPz9PQAGCBBQgSEABggQUIEhAAYKaDWiMgALpCChAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQNDqAgqQTvJGpX7BlEr/sIG+JG9U6hfkXms/cWH9Sy9BWW2tf1MLuw5tbUDpWf/SS1BWW+vf1MKuQ1sbUHrWv/QSlNXW+je1sOvQ1gaUnvUvvQRltbX+TS3sOrS1AaVn/UsvQVltrX9TC7sObW1A6Vn/0ktQVlvr39TCrkNbG1B61r/0EpTV1vo3tbDr0NYGlJ71L70EZbW1/k0t7Dq0tQGlZ/1LL0FZba1/Uwu7Dm1tQOlZ/9JLUFZb69/Uwq5DWxtQeta/9BKU1db6N7Ww69DWBpSe9S+9BGW1tf5NLew6tLUBpWf9Sy9BWW2tf1MLuw5tbUDpWf/SS1BWW+vf1MKuQ1sbUHrWv/QSlNXW+je1sOvQ1gaUnvUvvQRltbX+TS3sOrS1AaVn/UsvQVltrX9TC7sObW1A6Vn/0ktQVlvr39TCrkNbG1B61r/0EpTV1vo3tbAANRFQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgS0Jr8/L45eyq9MLn9+PbH3/M/jT+Ex3KLk918/de2ITw/jGv65a/LV1p5/wW0JsetZj2fm7mhGbOAbqefwvwz1ber9V/XhnBZ2/NPoJn3X0Brst+s6XNzbTv7+LzsNr98pHp3tf6r2hDm/1pMP4J23n8BrcluBR+Xtx2OWc8flvGI7s9/jh35WnS5crle/3VtCLvTsfr29HY39P4LaE221f+Du5DjObD5AdzX09erP4hL4dX6r2pDGHdAp38tdtPPoKH3X0ArMuyGjP/urs7hGO4/LucAhz+ffrtr4ULCvV6v/7o2hP1lL3PcD39q6/0X0IoMG07lByzLGI/hnmYXUYb9jtPPYfht/yl5vf7r2hDmjdweAtrS+y+gFRn+LX467I1UftiS2m78vMwCMvtIXV+b79Tr9V/thjAEdFzjlt5/Aa3IbvPlP6eLj3X/s7uE64A+Xb66kojMU7HWDWH4R2Nc35befwGtyGnw2+o+OKNZQLaza9Dbyj9AycwDutYNYdr1bOn9F9B6jKfQj1vLvombMNK6CujlQ1P7ByiZ2fqvdUM4nfBs6f0X0HoMBzCnjWXYlOo+9ZPeuwFdyYjIn1ejENa4IZzXtaX3X0DrVPt2k5490DdTuZ4NYX8eCNvS+y+gdap+/FtyAvpmQFezIVz62dT7L6B12q/lc3PmKvybAV3LhrCdXS5r6f0X0Dqt5XNzMQvIvqFxgMmsOqDjZbPLrQMtvf8CWqdVzSZxsOI7kQ7eP4Tvf0MYpwJ4vPpjM++/gNbj+gaMuo9c0vu53nvhD969E6v/DeH1lCEtvf8CWo/ZkJV99dN4JfdqIHkrs/Ekc70HvqYNYTZqa9LQ+y+g9RhPBB0/OLv13QV9FdCG5oNM5tVA+vVsCNMUTHMNvf8CWpHjrJBHdR+4LODqHGA7M5In8/ofkNVsCPPZ908T8Lfz/gtoTS4fnP4vHLzmmUiz9V/RhjB/ft5lhZt5/wW0Lsd/jis/bFnEq6vQrTyVMZlX67+aDWH+RKTZvxitvP8CChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECSuV+fNv8+U+C74H0BJTKCSj1ElAqJ6DUS0CpnIBSLwGlcgJKvQSUygko9RJQKneK48/vmz/+ftltBpdYDv9xs/k6D+h+/IbN0+H3zw/jXzn+3S9/5V90uiegVG4e0P992GxmgTzmdAjq/50C+nz6hmM4h//+OP8V0hJQKjcL6MVxf3J3+uO/T99z7udU0GnP0xE+CxFQKncd0HHXc+zm4/G/HP7T/nRYP37LIa37w4H98TfDr1sH8CxDQKncPKDHA/OxoF+Pv1wO3KeSTnuaw186NnNs536qKaQmoFRuHtDpROYQzOFLsytDU0pne5r76ZuHv/0f30/hhcQElMrNAnrK4zGgw384hfH8hfOpzuErx73O3eyaEyQmoFTuo4Ceenn87WFM08XxP447rg7gWYiAUrl3A3r839n3zK7Bz0aLbgWUxQgolfv8HuisqBd7h/AsR0Cp3E3nQH+5WjT8rX/7H6NAWYiAUrl3Azqe3bzckXT9hbPhvzzt3YfEQgSUyr0b0Ms40NOI+t15pOhpGNPhO90Jz1IElMq9H9Cxm2MwDxePTl847oNup1uSjiNDz2OaIC0BpXLvB3SaeWnwX99efWG62fN0D9LWdSQWIaBU7oOAnoL5xmxMh37+PN2DZDYRliGgVO6jgB5nGLmaD/Q4Q9P5LqTHl1e/g4QEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAARMXVJoAAACJSURBVIIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCDo/wHm12F0qczNAAAAAABJRU5ErkJggg==)
plot(predict(fit1),resDev,xlab=(expression(hat(eta))),ylab="Deviance residuals")
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABUAAAAPACAMAAADDuCPrAAAAV1BMVEUAAAAAADoAAGYAOjoAOpAAZmYAZrY6AAA6AGY6kNtmAABmAGZmZjpmtv+QOgCQkDqQkNuQtv+Q29uQ2/+2ZgC2///bkDrb////tmb/25D//7b//9v///9A0DGuAAAACXBIWXMAAB2HAAAdhwGP5fFlAAAgAElEQVR4nO3dYWPaSHeAUWW7jeu8ddq+y7I48f//nUUCDLblAa7HM1f4nC/JJtpkgsVjCY1GwxMAIUPvAQAslYACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKEBQ7oAOANXUT1T1P7Gi3q82cFuqN6r2H1jTJ3zDAL4sAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIEtDZrT5noSrgtgjo3EYKClxAQGe2GV79BGCOgL7dZDj+9FMHAyycgJY2UVCgQEBLmwgoUCCgpU0EFCgQ0NImAgoUCGhpEwEFCgS0tImAAgUC+nYT05iAiwjozDYm0gOXENC5jdzKCVxAQGe30k/gPAEFCBJQgCABBQhaWEAf74Zh+POf5//+/XP44+8r/n8BBepZVkBX+6s73w+/IKBAP4sK6Op5ftHhIFRAgX6WFNDx/P1h++P6WFABBfpZUkBXw7e/pp9ss7kvqIAC/SwooNtYPn/2udoXVECBfhYU0F8/phP4nX1BywEdZnzW6ICvZ6kBHQv6XUCBnhYb0PFz0Hun8EBHCwro6Wegu/8cHgQU6GdBAR2nL50cgo6zmr79W0CBbpYU0HEe6P3Jf2+mDzUFFOhkSQGdZtCfnsZvBBToaFEBnQp6egw6HpMKKNDJsgI6Xjm6f/ELawEFellYQD9KQIF6BBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFODEsHfZttX/9tp/YE0CChQNwxUFFdAMrvmWB3yi5/fhRW9IAU3gqm95wOc5eRcK6BspE3Xdtzzg85y+By94Pwpod1d+ywM+j4CWZCzUlV8x4PMIaEnGQAkopCGgJRkDJaCQhoCWZAyUgEIaAlqSMVACCmmYxlSSMVACCnmYSF+QMVCmMUEibuV8X8pCmUgPiVhM5F05C+VWTlgmAc1AP2GRBBQgSEABggQUIEhAAYIEFCBIQAGCBBQuYaoZMwQULuBmB+YIKJzndltmCSicZcEX5gkonGXJQeYJKJwloMwTUDhLQJknoHCWgDJPQOEsAWWegMJZAso8AYWzTGNinoDCeSbSM0tA4QJu5WSOgMIl9JMZAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQW0D+ng3DH/+U/tvvIKAAvU0CuivH2M3x34Ow7e/av+VlxNQoJ42Ad0Mwx9/P/3+OUy2P+1FQIF6mgR0PPLcHnduf9i2cz0M97X/zosJKFBPk4Bum/l998OYztX0H9f79WOYc9XxrIAC9bQI6PbUffwEdPvD9PHnJngOL6BALi0Cui3fw+6HKXbRgO6vQV0R0LnNQ38zwIyGAd3sTuTjAZ2OQa/5/FRAgc/UMKCrYfrhaR2fCrr9kz42CUpAgXoafQZ6f/wIdBvB2EWk0fYs/kMT8QUUqKfJVfjVmM717iak8YPMh/ifvv7Q/33hP9fpPnCJZvNARw/TefyHjiG3x68f+d8v+uf6wBS4SJs7kdZTjsbyrT96I9LHDkEv+ec+h1NBgaJG98JvhudL8B85A/+wC/65J9kUUKDEcnalTRQUKBDQ0iYCChQIaGkTAQUKPjegNW5frzs4AQWqEdDSJgIKFAhoaRMBBQp8Bvp2E9OY4Au75j4aAZ3ZxkR6+LKuuhNRQOc2cisnfFHXHUD1COivf2X+DPTJYiLwZV35EV7DWzkXchGpB8WGHK68iNwmoKsFXYXvwGcGkETGgL5+mFG/9URSJspVK8giY0BX01JM62k9plXPfqYMqHlTkEbCgO4fa7x/mNy64xl8ykKZuQ9pJAzorx/TszQf76ZnIm1z6hT+lIBCGjkD+nD8YX8q30fGQAkopJFwGtO+nLuHc45Hoh96sOaHZAyUgEIe+SbSHw49V7tDz8c705hOCSgkku5WzsOh53r3QM2NeaAvCChkkm4xkf2F980wXUVaf+jBxB+TMVCmMcFStZpIP6Zz+8M2nZvBRaSXTKSHhWp3K+f98Y5O05hecisnLFObgP7+OVVzv0B9vwPQpAG1mAgsU6vl7DbTx5+HQ9FuNAqox4LKAEECCnAi3TSmPAQUKEo3kf7Nw41NpAdySnkrp4ACC5BzMREBBRYg4XJ2p/aLK/cioEBB8oD2vRVeQIGS9AG1Ij2Q1TYRzxfhUwbUakxAVi8v11ywefUBnNug51PlBBQoOGZTQN8SUKAg4TSml379cAoP5JT9CHScFGpBZSClhJ+Bvp5Iv1vZrgsBBQoSXoV/HdCOC4IKKFCQcB7oy4D2u4L0JKBAUcKAJiKgQIGAlggoUJB+GlNXAgqU5FsPNBEBBYqumMQkoAAvXNHPTw7om6WUe1+IF1CgHgEFCBJQgKA2n4Funu/fHJPabz1lAQUqahLQx7uT2zc37oUHbkOTgK5erL+0shoTcBMa3Qt/eta+8RkocBMEFCCoRUB//3x1Cm9FeuAWNPoM9OQQdG1FeuA2tLoKfzhtH3/qKjxwE9rMA12/mEdvHihwExotJrJJ8UQkAQVqarYa0/6uzq4P9BBQoCbL2QEECShAkIACBH3+cnYPbxe1cycScAsEFCBIQAGCfAYKECSgAEECChDUNqDjUiL91rJ7ElCgpkYB/fVj7ObYz743wwsoUE+zp3L+8fe4sHLv++EFFKin2Xqg2+PO7Q/bdq6Hk0d0tiagQD1NArpfhH6fTk/lBG5Do2cijZ+Abn+YPv70UDngNjR8Kuf2h6mcAgrchoYB3eyfJiegwG1oGNDDsznXHmsM3IRGn4HeHz8C3ebURSTgFjR6Lvw2nevdTUjjlKZ+j+UUUKCeZvNA988z3p7H97yZU0CBelo+F34M57rvgzkFFKin4XPh95fg+52/PwkoUJPl7ACCBBQgqFlAxwcjjdfiuy4HKqBARc3WA92vBLruOYlJQIGa2gR0/1zObUBXQ9eCCihQT5OAjisp//nP4914J9LKgsrAjWgS0M20EOguoBZUBm5Fo1s5dzdxTgHdrw7ah4AC9bRbTOQQ0O0hqOXsgFvQcDm7Q0CtBwrcBgEFCOpwCr/yGShwExpdRBqPOfcBtaAycCNaTWP6fgjoOCfUgsrALWhzJ9I0e34K6NqCysCtaHor544FlYHb0GgxkfHEPUE/BRSoqNlydvuEds2ngAI1WVAZIKjNNKbOx51HAgrU02gifd916I8EFKin0a2c/abOvySgQD2OQGe32vns0QDL1upOpK5Pgz+66J97nLL66eMBlqzNVfjHuyTHoJf8c5/DqaBAUaPPQF9KvZzdSTYFFCgR0NImCgoUCGhpEwEFCtyJVNpEQIECAS1tIqAcmJnBDAEtbeLtwp65bcwR0NIm3i3smNvGLAF9u4lpTLxip2CegM5s42CDl5yWME9A5zbycRcvCCjzBHR2K/3klIAyb3kBXT0fHwbWyLPzEyGgzFtYQNev7mm6tqF2fiIElHmLCujj3fDalXeF2vmJEFDmLSmg04M9T4o53WN/3TJ5dn4iTGNiXrOAjrX79tfT+gPrgm7e5HJM6lVLNdv7CTG3jVmNArpbkGkM6AcWp1+9PWG/9nFL9n5izG1jTpuA7he02wZ0vIYeLOj2cPP+zS+urzuHt/sTpJ/MaBLQ8Uz7z38e77YBHQsaXA50W+G36d2U/rQ3l5zs/0BFrR4qdz9eQx8DOp7Evz2QvISAArk0CehqOs/eBzT8kGOn8EAujZ4LP5ZvH9Bt84Ln8C4iAak0eibSeO59CGjxrLtkfhrTVZ8HCChQz5ICOjuR/ro/S0CBejqcwq+Cn4HOPN5zmhl11eAEFKim0UWk8ThxH9BrP7Z89SdZTATIotU0pu+HgF599+VrlrMDkmhzJ9I0e34K6Prq9T+qElCgnqa3coZWoKtKQIF6Gi0mMl1A799PAQUqarac3T6hXfMpoEBNS1pQuQIBBeoRUICgVgFd7+Yu/frR9xxeQIF62l2FnwI6PhYuPo3+4wQUqKfZgsq7RT82wQnwtQgoUE+TgK5P7llfx5/pUYGAAvU0Wkzk5JPP7Vl8v0NQAQXqabSc3eminVcuIl+VgAL1NFwP9CC8HmgFAgrUI6AAQY0+Az09aY8vqPxxAgrU0+oq/PFD0E30scY1CChQT5OAjvPo9/OYxp92vBlJQIF62tyJtHnxHI5+00AFFKio0b3w4y2cscfA1SWgQD3NVmPaL0rfNZ8CCtRkOTuAIAEFCBJQgKD2F5G6PhhJQIF62gR0/WIak4ACN6FJQF9OAxVQ4DY0Ceiq++OMDwQUqKf9gspdCShQT4cFlXsSUKAeAQUIanQKL6DA7Wm0Hmi/JZRfElCgnlbrgSY5BBVQoJ42E+kf75IcgwooQYdJzL3HQSqNLiKZSM+yHffe3iMhEwGF857DqaCcElA46ySbdiFOWM4Ozjrdb+xDHAkonCWgzBNQOEtAmdcjoL/+5TNQFkVAmdcooBsXkVgwAWVem4CuXIVnyQSUeU0C+vKJSMPwUPvvvJidnwjTmJjXakX67+OSIt93P+/XT3s/MSbSM6vRcnbjnfCb3an7uufy9PZ+YtzKyZyGCyo/3n3762nKqVN4Fkc/mdEooA/HH/an8n3Y/4F6Ggb0sDD9412/pe0EFKinYUCfVrtDz8c705iAW9DymUj7J3tszAMFbkKjZyJNydwM01Wknk9IElCgnlYT6cd07h7ssRlcRAJuQrtbOe+Pd3SaxgTcgjYB/f1zquZ+afp+B6ACClTUajm7zfTx5+FQtBsBBeqxoDJAkIACBAkoQNDnBnS8avTgscbAbRJQuITVmJghoHAB64Eyx2egcJ4V6ZkloHCWZyIxT0DhLE/lZF6b9UB73nz0gp2fCAFlXqMFlbs+ivOEnZ8IAWVes4AOu8VAO7PzEyGgzGvzGejjXe/5S3t2fiIElHnNLiJt9rPo+q1GP7LzEyGgzGt5FX69b6j1QFkY05iY13ga09qdSCyRifTMaj8PdCOgLI9bOZnTOKArR6Ask34yo2VAD8+U6zit3v4P1NMsoAmuID0JKFBTm4Ae6tl3DtOTgAI1NbwTqfss+icBBWpqFdAM93E+CShQU5uA5lhJ5ElAgZqsBwoQ1Cyg+/P4tXvhgVvRKKC760hjQPuuDCqgQD1tAnpcEXTVd21lAQXqaRLQ3z/HKaCPd+Ol+FXX+UwCCtTTJKCb6fbNXUDHk/h+93ImDaj7rGGRmgR0Nd2CtA/o9nC034WknI2y0g8sU4uAbpM5HnPuA7o9BLUa0wvWmoSFanQn0njd6BDQjYC+YLVzWCoB7c7zdmCpOpzCr3wG+oKAwlI1uog0HnPuA7o9Hu23JmjGQAkoLFWraUzfDwEd54T2m0mfMVACCkvV5k6kafb8FNB132WVMwZKQGGpmt7KOXRfWDljoAQUlqrRYiLjiXuCfqYMlGlMsFTNlrPbJ7Tzcz1SFspEekjkmvsCLaicgFs5IY2r3o4CmoF+QhLXnRA2W84uwUPhn9IGFMjhyksSDQK6Hk71fb6cgAIFV06K+fSAnlx/3+n6UCQBBQqSBXTXz8Nh57p3QQUUKEgW0NfPQNr0/SRUQIGCXAF9vBt2SzAVfqUlAQUKcgV05gFInokEZJUqoL9/vr31yHJ2QFappjHNPUDOQ+WAtDJNpJ892lx5pAeQVaJbOQUUWJg8i4kIKHDDBBQgSEABggQUIEhAAYI+PaBzBBS4BQIKECSgAEGeiQQQJKAAQQsKaI3PAwQUqEdAAYIWFNBxMXsBBfJYUkCnY9CPrWYvoEA9iwroWNCPPVBJQIF6lhXQ8Sz+Q8vZCyhQz8ICOj6S7uH8Vu8SUKCepQV0exJ/+SHo3DWnzxwc8LUsLaBXHYIKKPCZFhfQjxFQoB4BBQgSUIAgAQUIWmxA16Ep9QIK1COgAEECChAkoABBAgoQJKAAQQIKECSgAEECChC02IDGCChQj4ACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECWgGw17vcQBXEdAEhkFBYYkEtL/ncCooLIuAdneSzYzDA94loN2djinj+ID3CGh3AgpLJaDdCSgslYB2J6CwVALanYDCUglodwIKSyWg3ZnGBEsloP2ZSA8LJaAJuJUTlklAM9BPWCQBBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABnd3K2h7AeQI6t5GCAhcQ0JltrG8MXEJA327iCRvARQS0tImCAgUCWtpEQIECAS1tIqBAgYCWNhFQoEBAS5sIKFAgoKVNBBQoENC3m5jGBFxEQGe2MZEeuISAzm3kVk7gAgI6u5V+AucJ6NxGCgpcQEBntvEZKHAJAX27iavwwEUEtLSJggIFAlraRECBAgEtbSKgQIGAljYRUKBAQEubCChQIKClTQQUKBDQt5uYxgRcREBntjGRHriEgM5t5FZO4AICOruVfgLnCShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAKcGPYu27b63177D6xJQIGiYbiioAIK8Ow5nBcVVEABDk6yKaBvCChQcJqIC3IhoAAHAloioECBgJYIKFAgoCUCChQIaImAAgUCWiKgQIFpTCUCCpSYSF8goECRWznfJ6BAmcVE3iWgQD0CChC0sIA+3m2PrP/85/m/f/8c/vj7iv9fQIF6lhXQ1f7Die+HXxBQoJ9FBXT1fHnscBAqoEA/SwroeP7+sP1xfSyogNLINZdm+TKWFNDV8O2v6SfbbO4LKqC0cdXkQL6MBQV0G8vnzz5X+4IKKE1cd3sKX8aCAvrrx3QCv7MvqIDSwpU3SPNlLDWgY0G/30xAnR0md+USPXwZiw3o+Dno/Y0E1Odr2Qko8xYU0NPPQHf/OTzcREB9vpaegDJvQQEdpy+dHIKOs5q+/fsGAurztfwElHlLCug4D/T+5L8301lvIaDDjE8bXZg3Z36+RsxbUkCnGfSnp/EbAaUJXyPmLSqgU0FPj0HHY9IbOIWf/zl5+JiFecsK6Hjl6P7FL6wFlAZc6GPWwgL6URn3fgFdgswfAtGPgHYnoIugn8xYbEDXh5VFrpJx//f5GiyVgPbn8zVYKAFNwOdrsEwCmoF+wiIJKECQgAIECShAkIACBAkoQNBiAxojoEA9AgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKEPTlAgpQT/VG1f4Da+r9YgO3pXqjav+BnyHzqXzisSUeWuaxGVpI4rF95tDy/qtPJP7aZB5b4qFlHpuhhSQem4AmHmXisSUeWuaxGVpI4rEJaOJRJh5b4qFlHpuhhSQem4AmHmXisSUeWuaxGVpI4rEJaOJRJh5b4qFlHpuhhSQem4AmHmXisSUeWuaxGVpI4rEJaOJRJh5b4qFlHpuhhSQem4AmHmXisSUeWuaxGVpI4rEJaOJRJh5b4qFlHpuhhSQem4AmHmXisSUeWuaxGVpI4rEJaOJRJh5b4qFlHpuhhSQem4AmHmXisSUeWuaxGVpI4rEJaOJRJh5b4qFlHpuhhSQem4AmHmXisSUeWuaxGVpI4rEJaOJRJh5b4qFlHpuhhSQem4AmHmXisSUeWuaxGVpI4rEJaOJRJh5b4qFlHpuhhSQem4AmHmXisSUeWuaxGVpI4rF9+YACZCSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChCUN6C/fgxb97O/txl/a3hoPKKj0thGj3fdBlcY2nZUW9/+aj2kp+Kozr2Wny3nCzZJvJtlfnue/Yp+r/iXpQ3oatiZ2X13r8LWH393GNhTcWyT8SvYafd5f2i73arPy1Z4wc69lp8t5ws2SbybZX57vj+03z8PX9F6L1vWgK4P/9S3X4bnL9Aw/PlPsrHtrLp9/31/aMcctN+zCy/Y2dey29B6vmBnhrbXbzfL/PYsDG31/Fv1XrekAR2/CuPLP54MvDrgHnfs6SszvlI9zvwKY9vZ1P0ed4XC0J5frdV7w+4xqrOvZb+hdXzBzg1tp99ulvntWRja4ZfGIVb7lpg0oKvDv377erw6FF8/f2fb9PkeVxjbZHfo0mXPfn9oJ6d769ZHVIUX7Nxr+dlyvmBnhrbTcTfL/PYs72y7otf87CNnQLf/wsMuu371fez3z+fX5eSnDRXGNtmO6o//6bNnF4a2OX4/Hj8Jajm6wqjOvZafLecLNkm8m2V+e5ZettXzeCrubDkDuv3u8f340xffx05+q4/C2CbbL87Dus+eXRja6S6zaju6y76Ys6/lZ8v5gh3Gk3U3y/z2LL1sLwJ620egJ3vv+J329Pxp03PKy6gwttH2y3Zf8wt0jTNDO1i1PTIojOrCAX+anC/YJPFulvntWXrZnn+v5rFx2oAe9ozX/9jpt3bzEbpd6H5vbLtf2n7b6xfQwtAOtmc5TY/1zn0x53+riZwv2CTxbpb57Vl62caPPr+f/FhFzoCenjO9+vY//udhpkKX2SWFsT3/Sqc9uzy0g9afNha/mJcM+PPkfMEmiXezzG/P4st2nJlW72OGrAE9/tPffoX+u+v8vMLYns9gugW0MLSD5h82lr+Y5wf8iXK+YJPEu1nmt+eZr+iq8izQZQT09N+7OzuYduhOcwffH9vxZC9FQOeHsM1B4/26MKpLBvyZcr5gk8S7Wea3Z/ErepxIX+979RICevqvnb5C+xOqPnMHS9/jVvv3WoqAzr42m/bHBcs+Au3wgk0S72aZ357ll203tnGQN3gR6fkDiu2/rfieO55QNTtouXBszzt0wz37wqHt9cjBogPaq5/ZdrMXsr09X/z97w7tJOjrenP8cwa0cC1tdfKJfrM5E5eN7TgHrVNAz15UXvW4OWTJV+G7vGCTZLvZC9nenifKO9v9O7/1ATkDurlo6mCfgL4/tuMyBjtNxnbhy7b/tR4fGhdGdWbAny7nCzZJtptdOLY+b88ThaGdHpDW+9aTJ6CnCvcTnNxi12fW7vtj675nl+9eGZdZ6LT4ygLvRHrq94Id/vKsu1nmt2dhaF8poIU7Wk8vivb4kKUwtu57dvH+6U7LdSz1XviOL9gk8W6W+e1ZGNrLU/jbDuhFC/j0m2BydgWhjjOc3xvads/qlYNFrsbU8wWbJN7NMr893x/aZji9iFRraEkDesESktPvpVwPtN+e/f7QuqwodHZUidcD7fmCTRLvZpnfnuWv6PM0phtfjWlmXenjAcHJKUySFenfHKx027PfHdrm5Wlf29Fd9MVMsiJ9ihesPLTTLbKsSJ/n7fn+0E7Wyq83tKwBfb5r4LDLnOw9m65foPLYdvrt2e8M7fg0mC49KLxg2Z6JlOMFKw3tqN9ulvnt+f7QjvfC1xta2oC+frbei71ndfINpofS2EYd9+z5oZ0+4KdHDwovWLKnciZ5wQpDO+q4m2V+exaGtjsIrTm0vAEFSE5AAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQFm89TA89B4DX5OAsnS/fgzD996D4GsSUJZuMwzDt796j4IvSUBZuN8/h//8Mdz3HgZfkoCycI93w8N6+OPv3uPgKxJQFm41/PnPNqIOQelAQFm2qZ3b0/g//+k9Er4gAWXZ1tMFpI2ZTPQgoCzarx/Tsef2BzOZaE9AWZPxJhUAAAGjSURBVLTDoefKTCY6EFCW7PfP/fX3xzuHoLQnoCzZczefSwoNCShAkIACBAkoQJCAAgQJKAu22i1kNy5o5xISHQgoy7V6mKaBjuvZKSg9CCiL9et/pzuQNuMU+u2xqOVEaE5AWazH/xunf/7Hf42HnvtbOqEpAWXJjvPnVwJKewLKkm0Dur+DU0DpQEBZMgGlKwFlyQSUrgSUJRNQuhJQlkxA6UpAWTIBpSsBZckElK4ElCUTULoSUJZMQOlKQFkyAaUrAeU2CCgdCCi3QUDpQEC5DQJKBwLKbVhZUZn2BJTbIKB0IKDchrWA0p6Achs2w/hkD2hKQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGC/h8QDkgFBphRngAAAABJRU5ErkJggg==)
par(mfrow=c(1,2))
plot(jitter(as.numeric(dat$variety),amount=0.1), resDev, xlab='Variety',
ylab="Deviance residuals", cex=0.6, axes=FALSE)
box()
axis(1,label=c('O.a. 75','O.a. 73'),at=c(1,2))
axis(2)
plot(jitter(as.numeric(dat$root),amount=0.1), resDev, xlab='Root',
ylab="Deviance residuals", cex=0.6, axes=FALSE)
box()
axis(1,label=c('Bean','Cucumber'),at=c(1,2))
axis(2)
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABUAAAAPACAMAAADDuCPrAAAAaVBMVEUAAAAAADoAAGYAOjoAOpAAZmYAZrY6AAA6OgA6ZmY6kLY6kNtmAABmAGZmtrZmtv+QOgCQZgCQkGaQkNuQ29uQ2/+2ZgC225C2///bkDrbtmbb/7bb/9vb////tmb/25D//7b//9v///8PCQbyAAAACXBIWXMAAB2HAAAdhwGP5fFlAAAgAElEQVR4nO3d64IbN5MeYNrrSMnKiZxkP8aziSlb93+RYfPYPMwMCXajC6jn+WNphqQFVPElG31a/QSgyGrpfwBAqwQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClBIgAIUEqAAhQQoQCEBClAodoCuaMzSHVPX0rPNs6ZvgclfcUJLzzbPW7pnalp6rnne5D0w9QtOKdfbsQe5KpZrtD0QoISWq2K5RtsDAUpouSqWa7Q9EKCElqtiuUbbAwFKaLkqlmu0PRCghJarYrlG2wMBSmi5KpZrtD0QoISWq2K5RtsDAUpouSqWa7Q9EKCElqtiuUbbAwFKaLkqlmu0PRCghJarYrlG2wMBSmi5KpZrtD0QoISWq2K5RtsDAUpouSqWa7Q9EKCElqtiuUbbAwFKaLkqlmu0PcgdoHNcUJpJ5SrQh6PVrQGlDtCZrsnPhHKV56PR6taIMgfooRv1ZGS5qvPBaHVrSLkD9PZnBJOrOh8G6GePYAGJA/T4Zy0ZWa7qvD9a3RpT4gD1md6CXNXxDbQ1uQPUqlJ4uapjDbQ1mQPUfs0G5CqPvfCtSR2gjqyLL1eBHAfamtwBSni5KpZrtD0QoISWq2K5RtsDAUpouSqWa7Q9EKCElqtiuUbbAwFKaLkqlmu0PWgsQH98Xa1Wv/11+vs/f6x+/fOJ52vQ1rRVMf2ZTVsBuj4cC/fl+AMN2rumKqY/02kqQI/9ef6Q16C9a6li+jOflgJ02D76vv3v27lDNWjvGqqY/kyopQBdr3751+4P27Y8dKgG7V1DFdOfCTUUoNtmPK0trQ8dqkF7107F9GdGDQXo37/vNpD2Dh36cYOu7pjrX8c82qmY/syo1QAdOvSLBu1fOxXTnxk1G6DDOtM3m0jda6di+jOjhgJ0vMa0/+vquwbtXTsV058ZNRSgw+Eho4/44aiRX/5Dg3auoYrpz4RaCtDhOLtvo79vdotGGrRrDVVMfybUUoDujlAebyZtNGj3WqqY/synqQDddej4M374zNegXWuqYvoznbYCdFiZ/3bxgzcN2re2KqY/s2ksQF+lQVuTq2K5RtsDAUpouSqWa7Q9EKCElqtiuUbbAwFKaLkqlmu0PRCghJarYrlG2wMBSmi5KpZrtD0QoISWq2K5RtsDAUpouSqWa7Q9EKCElqtiuUbbAwFKaLkqlmu0PRCghJarYrlG2wMBSmi5KpZrtD0QoISWq2K5RtsDAUpouSqWa7Q9EKCElqtiuUbbAwFKaLkqlmu0PRCghJarYrlG2wMBSmi5KpZrtD0QoISWq2K5RtsDAUpouSqWa7Q9EKCElqtiuUbbAwFKaLkqlmu0PRCghJarYrlG2wMBSmi5KpZrtD0QoISWq2K5RtsDAUpouSqWa7Q9EKCElqtiuUbbAwFKaLkqlmu0PRCghJarYrlG2wMBSmi5KpZrtD0QoISWq2K5RtsDAUpouSqWa7Q9EKCElqtiuUbbAwFKaLkqlmu0PRCghJarYrlG2wMBSmi5KpZrtD0QoISWq2K5RtsDAUpouSqWa7Q9EKCElqtiuUbbAwFKaLkqlmu0PRCghJarYrlG2wMBSmi5KpZrtD0QoISWq2K5RtsDAUpouSqWa7Q9EKCElqtiuUbbAwFKaLkqlmu0PRCghJarYrlG2wMBSmi5KpZrtD0QoISWq2K5RtsDAUpouSqWa7Q9EKCj32nfeHLVJNdom/BJLAjQ869WIjSeXBXJNdoWfBYLAvT0m9XHD2AulRs0sidH6xN/bp/GggC9/I2GrO7jz/hcBXlutDaaZvdpLAjQy19ox9o++YzPVZCnRmujaXafx4IAvfyNbqztk4nPVZAnA/T55/Ac30AvWQON5rPP+FwFeWa0NpoqsAZ6yV74cHwDHfENNBp74S84DjQca6Aj1kDDcRzomGaLx174M3vhWyNAWZrjQE8cB9oaAUpouSqWa7Q9EKCElqtiuUbbAwFKaLkqlmu0PRCghJarYrlG2wMBSmi5KpZrtD0QoISWq2K5RtsDAUpouSqWa7Q9EKCElqtiuUbbAwFKaLkqlmu0PRCghJarYrlG2wMBSmi5KpZrtD0QoMzp5bO1c1Us12iDKepVAcqMXr9eUK6K5RptLGW9KkCZzwRXrMxVsVyjDaWwVwUo85ngmum5KpZrtKEU9qoAZTZT3LUnV8VyjTaS0l4VoMzHN9An5RptKL6BPkKDVmUN9Em5RhuKNdBHaNC67IV/Tq7RxmIv/AM0aGWOA31KrtEG4zjQz2nQ1uSqWK7R9kCAElquiuUabQ8EKKHlqliu0fZAgBJarorlGm0PBCih5apYrtH2QIASWq6K5RptDwQooeWqWK7R9kCAElquiuUabQ8EKKHlqliu0fZAgBJarorlGm0PBCih5apYrtH2QIASWq6K5RptDwQooeWqWK7R9kCAElquiuUabQ8EKKHlqliu0fZAgBJarorlGm0PBCih5apYrtH2QIASWq6K5RptDwQooeWqWK7R9kCAElquiuUabQ8EKKHlqliu0fZAgBJarorlGm0PBCih5apYrtH2QIASWq6K5RptDwQooeWqWK7R9kCAElquiuUabQ8EKKHlqliu0fZAgBJarorlGm0PBCih5apYrtH2QIASWq6K5RptDwQooeWqWK7R9kCAElquiuUabQ8EKKHlqliu0fZAgBJarorlGm0PBCih5apYrtH2QIASWq6K5RptDwQooeWqWK7R9kCAElquiuUabQ8EKKHlqliu0fZAgBJarorlGm0P6gboj6+r1W9/Tf1/fIIGbU3ViulPnlQpQP/+fejLoT9Xq1/+NfX/8nEatDV1KqY/KVMnQDer1a9//vznj9XO9o9L0aCtqVIx/UmhKgE6fLJvP9e3/9n25ttq9W3q/+fDNGhralRMf1KqSoBue/LL/j9Da653f3ne37+v7nnq+4IGbU2NiulPStUI0O2m0bDCtP3PbnlpU7iNpEEzqlAx/UmxGgG67azv+//smqm0QQ9r/E806L2HF/2fWUyFiulPilUM0M1+Q6m8QXef8c+sT2nQ9tULUP3J8yoG6Hq1+8/Pt/JD7bav9NpBJhq0NfUCVH/yvEproN/OS0zbJitbpB9st5JeOtBZg7amzhqo/qRMlb3w66E13/YneQwLRd/LX/3tpWdr0ObUqJj+pFS140AH33fbSS99Rm+/H7zydA3amlrHgepPStQ5E+lt159DZ729eqLHax/xGrQ1VSqmPylU6Vz4zeq0i/OVLZyXPTtc+0WXVmf+W+3PD15I41bhcnYfP9yhJQvLNftTjVbj1iJAP3z06vknMalckz/RaDVuNQL080frwwXlmvzJAnTKV+MD8wboFKcHT/uPe2a4xwfrwwXNOvlt9+enr6Jx5ydAP3+0PlyQAC1/GY07PwH64aMtJS1NgBa9jMatxBroxw+3M3NhuWbfXvjWCNBPHq8Ll5Vr/h0H2hoBSmi5KpZrtD1YIkD//u+trzGNX1DPz2qB6e2pP/nIBO/eiqdydrRIP3o9a03zqncqZ4/9yUemePfWCdB1Z3s5zy9nb+fMqsxtr/3JRyZ591a9nN3RctdrmDxA53hVRmpezq67/uQjk7x7K11QebjUzdvuejfrJftz4uE64+MVj2091bmgcvz+tFI0tWnevZVu6TFcavFws663BbeQfAON48EFqDq39Ajfnxbbp9fMN9C/f9/dq/DH1909Z7bt2s0mkjXQYo9OXYW5baA/NdoMmlkDPd93+3DXw/Kbdr3KXvgoHv34rxOg0fvTps4cWtkLf+jM/c0Ph0/6l25c+BLHgQbx8AJUvQAN3J8W2+fRyHGgx4/29f6j/cdXh4kQ7xto5P70DTSqSjuRdh/tb/sbFm4cZ0ekNdAG+tMaaFRVDmM67NjcrHar9G8v3fj1NTowjDh74VvoT4vtQdU6kH5oze1/tq25WQVcpKe+OMeBttCf4jOmeqdyfjufMRfvMBGiqnYqp/6kQJ0A/eePXVceLgC+3Ae8Bm1OlYrpTwrVupzdZre8dPyoX4wGbU2liulPirigMqHlqliu0fZAgBJarorlGm0PBCih5apYrtH2oNKZSC5YS5k6ZyLpT8oI0NIX0utVCFBCeOcdL0ALX8eJIXUIUCJ47x1few30cPHapUw0XKcmV1N5jvvoT6b27ju+/k6kJU81ni5Ap3w1PlB9jnvoT6b27ju+foCGvOJ32ato9/lVn+MO+pOpvf+OX+AwppBXuyl6Ge0+v/pz3EF/MrVA30AXvWuXNdDWLBGgzfcnUwu0BhqvQQv2p9sLX4sApYJP381R9sLvDhqJtYlUFIbis5Lq0xyuP5nfAxmw4HGgF4aD7kJdsNbmeGi16xKuP5nfCxmwxIH0+yuHLeL97+BaN6YFDqSP1Z/M74UMWCJAF7zg4vuLwFo3pgUCNFR/Mr9XMqB+gC63Qv/TN9D2VA/QaP3J/IJ/Aw3EGmhrctUl12jDCL4GGshke+GpJFdhco02jvIMEKAOSQotV2lyjTaQ4gwQoISWq2K5RtsDAUpouSqWa7Q9mDdAby5Vu/SOTg3amlkrpj95kQAlNAFKZAKU0AQokdVZA92czo8bWna569Vq0OZUqZj+pFCVAP3xdXR63Ma5xjyuRsX0J6WqBOj64vo2a1e74WE1KqY/KVXpXPjxVtHGGhMPq3MuvP6kjAAlNAFKZDUC9J8/rjaRXPGbR1WomP6kWKU10NFH/JsrfvO4Omug+pMytfbCHzeLhj/ay8nDKu2F158UqXMc6NvFccqOs+NhVSqmPylU6WIim3N7Lvj5rkHbU6di+pMy1a7GdDhrbtEbJmjQ9tSqmP6khMvZEVquiuUabQ8EKKHlqliu0fZAgBJarorlGm0P5r+c3ffbi4Y504NHzX45O/3JCwQooQlQIhOghCZAicwaKKHlqliu0fZAgBJarorlGm0PBCih5apYrtH2oG6ADpdqWO5aYT81aHuqVkx/8qRKAfr370NfDv257MnGGrQ1dSqmPylT7a6cv/45XLh26fONNWhrat2VU39Sotr1QLef69v/bHvzbTW6BWJtGrQ1ta4Hqj8pUSVADxf5PrSmux7yuBoV05+UqnRPpGGFafuf3fKSm3bxuDr3RNKflKl4V87tf3adqUF5XL27cupPnlcxQDeHu3VpUB5XL0D1J8+rGKDHex++uW0sD6sXoPqT51VaA/12XmLatqtFeh5VZw1Uf1Km0n3ht635tj/JYzhkZLnbHmrQ1tS5L7z+pEy140AP94vdbictebKcBm1NreNA9Sclat4XfmjMt2VvfKhBW1PxvvD6k6dVvC/8YRfncttHPzVoe+rdF15/8jyXsyO0XBXLNdoeCFBCy1WxXKPtQbUAHW48M+zrXPRyixq0ObUqpj8pUe16oIcrLb4teZCIBm1PreuB6k9K1AnQw30Ptw26Xi3aoRq0NVUqpj8pVCVAhyvV/vbXj6/DmR5rF6zlCTUqpj8pVSVAN7sLLe4b1AVreUaNiulPSlU6lXN/ktyuQQ9XX1yGBm1NnVM59Sdl6l1M5Nig2494lwvjUdUuJqI/KVDxcnbHBnW9RR5X73J2+pPnCVBCE6BEtsAm/NoaEw+rvwmvP3lcpZ1Iw2f6oUFdsJYn1NmJpD8pU+swpi/HBh2OuXPBWh5V6TAm/UmROmci7Y5O3jXomwvW8owqFdOfFKp6KueeC9byuJqncupPnlXpYiLDhlGA/tSgzalTMf1JmWqXszu06KLtqUHbU6ti+pMSLqhMaLkqlmu0PahzGNPCn+tnGrQ1VQ5j0p8UqnQg/bLX+T7ToK2pcyC9/qRMpVM5lzs0+ZIGbU2dUzn1J2V8AyU030CJrNaZSIvebftMg7am0plI+pMidfbC//ga5DNeg7amSsX0J4UqrYFeavByYcM/e9J/CY+pswbaSn9qw2gE6IPPW+ndRQjQi19rw2AE6GNPW73wbF4gQMe/1YbROBPpiafp3PpyzflnAfrAg6hKgD7xLJ1bX645/3C02jAgAfrE03Rufbnm3DfQ1gjQx55m8WkhuebcGmhrBOiDz7P7cy4fT2yuSbcXPrib6Regjz5R387jk1DINe2OA43ttlkFKIv6bLM0V8VyjbY5d5pVgLKoz3aM5KpYrtE2506ztheg69PxzgXXINOgwXx6aE5zFdOf3brXrI0F6NvVOSPP9qgGjaavb6D6s2utfwP98XV17cmz7jRoND2tgerPzjW+Brq7ceKoI3fnMD93GTINGk4/e+H1Z/cW3As/dNMv//r59sJ1Fzc37Ti07FOXwtWg8cQ4DlR/8oCljgPdX/BmaNAXLv69vt0gevZ2Nhq0NXUqpj8pUydADxcM2zbosI+ysEO3H+ffbn749tw2kgZtTZWK6U8KVQnQYUvmt79+fN026NChhZdb3Hb5bWtvPnq1myV953E0p0bF9CelqgToZjV8Nu8bdNhIuv2gfoQGzahGxfQnpaoE6Hq3HXNo0OKbyNpEyqhGxfQnpWoE6KGzDg267anCbSSL9AlVqJj+pFiNAD1s2xwb9MOtmo/cP0zkqe0tDdqaChXTnxRrKUDvHqj83Gtp0NY0FKD6M6EFNuHXhWtMd26fuDvy5Kl/nAZtTP1NeP3J4yrtRBo+hw8N+uyy0NUrXXKxht7V2YmkPylT6zCmL8cGffrstmsuF5ZKpcOY9CdF6pyJtDs6edegw/W+Xjjd+FUatDVVKqY/KVT1VM69wiX6SWjQ1tQ8lVN/8qxKFxPZ7aBcvj81aHPqVEx/Uqba5ewOLbpoe2rQ9tSqmP6kREsXVJ6ABm1NrorlGm0PBCih5apYrtH2oFaAvu2PDfn792W3kTRoaypVTH9SpN5e+F2DDrfdKj9M+XUatDXV9sLrTwpUu6Dy/qIKm8IDjKeiQVtT64LK+pMSVQL0bXRO8HCk8ktnerxEg7amRsX0J6UqXUxktLK03Upa7iNeg7amzsVE9CdlKl3ObnxRxCcv0j0pDdqaOpez05+UqXg90KPi6y1OQIO2pt71QI/0J48ToIQmQIms0hroeKOo/IK1r9OgramzBqo/KVNrL/x5kWlTetvYKWjQ1lTaC68/KVIlQIfjlA/HiQx/XPBkDw3amhoV05+UqnMm0ubiPgfLHWanQZtTpWL6k0KVzoUfTpEru83WtDRoa+pUTH9SptrVmA4X/V60PTVoe2pVTH9SwuXsCC1XxXKNtgcClNByVSzXaHsgQAktV8VyjbYH9XciLXrjGQ3amuo7kfQnT6gToG8Xh4loUB5WpWL6k0JVAvTyMDsNyuNqVEx/UqpKgK4Xv13skQZtTY2K6U9K1b+g8qI0aGuqX1B5UfqzNQtcUHlJGrQ19S+ovCT92RoBSmgClMgqbcJrUMrU2YTXn5SpdD3Q5S5Re0mDtqbO9UD1J2VqXQ80yEe8Bm1NpeuB6k+K1DmQ/sfXIJ/xGrQ1VSqmPylUaSeSA5UpU2cnkv6kjAAlNAFKZAKUGQxVnuqVpnmdD+jPrCZoU5ezY3qHIJrmpaZ4lVbkGu3SpmhTAcrkDk05yWTnqliu0S5skjYVoEzuMMsC9Gm5RruwSdp0iQD9+79bY+rZcZJbDVD9mcI0bVopQDcW6RNp7xuo/kyooW+ga3s5M2luDVR/ZtTOGujlHWdWq+9T/z8fpkGraGwvvP7MqZm98NsP+C/DJRu+7P+8XH9q0EqaOg5Uf2bVyHGg//yxO9N4s980elvy8t8atDV1LmenPylT8YLKP77+8q+fu3a1icSj6l1QWX/yvEoB+v38n8Om0jI0aGvqBKj+pEzFAD1e+PvH1+UuHaZBW1MvQPUnz6sYoD/X+4/2H18dJsKj6gWo/uR5Ne+JdLhzwsZxdjys4j2R9CdPq3RPpF1Lbla7Vfol70CjQVtT555I+pMytQ6kH1pzf+OEzcoiPQ+rdCC9/qRIvVM5v53PmHOYCI+qdiqn/qRAnQD9549dVx4u/b3cB7wGbU6ViulPCtW6nN1mt7x0/KhfjAZtTaWK6U+KuKAyoeWqWK7R9kCAElquiuUabQ8EKKHlqliu0fZg3gAdVuW/u20s5WatmP7kRQKU0AQokQlQQhOgRGYNlNByVSzXaHsgQKnrydso5KpYrtEuZro7zghQ6nr2Rl65KpZrtEuZ4mZy59ea5GXGr3jzk79/X/LkjgsadGGrZ28lW+V6oPozk6d78OMXm+JVLl7x5ie7RfoFb3U4okEXdihAsADVn5k83YMPvNiE3gnQ1f5iiwvToMs6zn+8ANWfWTzfg4+82nTuveCPr0sfH3KgQRcW8Ruo/syluW+gO5vDyu1yV/seaNCFRVwD3dGfaTS3Bnr0duhR11tMLPBeeP2ZRGt74UfenOmRXejjQPVnCk0fB7rRoDyuesX0J0+oHKBrn/A8pW7F9CfPqRmgx3t2LXjYsgZtTcWK6U+eVi1AA6zQ/9Sg7alVMf1JiToBeuzOZY8R+alB21OlYvqTQhXPRFr8KOWfGrQ99c5E0p88r1aARjhP7qcGbU+lANWfFKkToDGu1PBTg7anSoDqTwq5Hiih5apYrtH2oFqAHraT3pxrzDNqVUx/UqJSgO7X6YcGXfbKixq0NXUqpj8pUydAz1dcXC977VoN2poqFdOfFKoSoP/8MRxi9+PrsKtzvejxIhq0NTUqpj8pVSVAN7vT4/YNOmwkLXeunAZtTY2K6U9KVQnQ9e4Uj0ODbj/ul1uo16CtqVEx/UmpGgG6bcnhM/3QoNuPeFe74VEVKqY/KVbpTKRhXf7YoBsNyntuLnVb50yk1vpzwisC8xIBSiC3N1sQoHcfNt09KXjJApvwa2tM3Hfndl/1N+Hj9+ekd0XjJZV2Ig2f6YcG3X7eL3fNRT0X2p0bztbZidRWf056X15eUuswpi/HBh2OuVvuSGU9F9mxOrUDtLH+vDdNLKTOmUi7o5N3Dfq27GVr9VxoC30Dba0/fQONo+qpnKvFL1yr50JbZg20uf60BhpHpYuJDBtGAfpTzwW3yF74n831p73wYVS7nN2hRRe+b4KmC26J40B32upP8RmFCyoTWq6K5RptDwQooeWqWK7R9qDa5ewC3HT758OHiejjMGpdzq6h/mRWz739KwTo8abbe8vev+uR4Vqhj2T+QrTWn8zqybf/7AE62r+5t+hNZx4YrmNEQpm7Ds31J7N69u0/d4Du+/P4sf62dIc+FKCPPpIKZq5De/3JrJ59+88doNf3mNksu9L0+XCdJxfLzHVorj+Z1dNv/5kD9MfX1f4SNx/8pCbfQFszbx0a7E9mFewb6J0bzES/54w10FDmrUOD/cmsYq2B/vPH7akd4S8XZi98JLMWosn+ZFah9sLfu0FX/Jt2ic9A5g7QBvuTWUU6DvTup/k6/i0TCGPWiulPXiRACU2AEpkAJTQBSmQClNAEKJEJUEIToEQmQAlNgBKZACU0AUpkswfoPRqUR80doPqTVwhQQhOgRCZACU2AEpl7IhFarorlGm0PBCih5apYrtH2oKEAnWJ7S4O2pp2K6c+MBCihtVMx/ZlRQwE6XCxcg2bTUMX0Z0ItBejuM/61q4Vr0Na0VDH9mU9TATp06Gs3rNGgrWmqYvoznbYCdNhKeuly4Rq0NW1VTH9m01iADrf8+v75o96lQVvTWMX0ZzKtBeh2I+nxj/h7a/pz/uOYXmMV05/JtBagT33Ea9D2tVYx/ZlLcwH6Gg3amlwVyzXaHghQQstVsVyj7YEAJbRcFcs12h4IUELLVbFco+1BswH6VnTIsgZtTasV0585CFBCa7Vi+jMHAUporVZMf+YgQAmt1YrpzxwEKKG1WjH9mYMAJbRWK6Y/cxCghNZqxfRnDgKU0FqtmP7MQYASWqsV0585NBugZTRoa3JVLNdoeyBACS1XxXKNtgcClNByVSzXaHsgQAktV8VyjbYHApTQclUs12h7IEAJLVfFco22BwKU0HJVLNdoeyBACS1XxXKNtgcClNByVSzXaHsgQAktV8VyjbYHApTQclUs12h7IEAJLVfFco22BwKU0HJVLNdoeyBACS1XxXKNtgcClNByVSzXaHsgQKlhtSqc+lwVe2C0xTPJHAQoFaz2ip45+T8msM9HWz6TzEGAMr/DO75k9nNV7NPRvjCTzEGAMr/DtAvQzzwQoI89jkoEKLM7zroA/cxno31lJpmDAGV+voE+yDfQ1ghQ5mcN9EHWQFsjQKnAXvjH2AvfGgFKDbd7IdcAABOMSURBVI4DfYjjQFsjQAktV8VyjbYHApTQclUs12h7IEAJLVfFco22BwKU0HJVLNdoeyBACS1XxXKNtgcClNByVSzXaHsgQAktV8VyjbYHApTQclUs12h7IEAJLVfFco22BwKU0HJVLNdoeyBACS1XxXKNtgcClNByVSzXaHsgQAktV8VyjbYHApTQclUs12h7IEAJLVfFco22BwKU0HJVLNdoeyBACS1XxXKNtgcC9Jkn6+/qck35i6PVodUJ0Cee62409eWa8NdGq0PrE6CPP9X9EBeQa75f3UJ6+TV4kgB98qnas65c8/1igL7+GjxJgD75TO1ZV675fm2J6fXX4FkC9Mmnas+6cs23b6CtEaCPP9UK0wJyzbc10NYI0Ceeax9nfbkm3F741gjQZ56sOavLNeWOA22NACW0XBXLNdoeCFBCy1WxXKPtgQAltFwVyzXaHghQQstVsVyj7YEAvfitBo4mV0UKRqtpFyVAx790FEg4uerx/Gg17bIE6Oh3jkOOJ1c5nh6tpl2YAL3+nV4MJVc5CgK07HlMRIBe/0ovhpKrHM+OVtMuTYBe/04vhpKrHL6BtkaAjn5nOSmeXOWwBtoaATr+pR2a4eSqh73wrRGgF7/VidHkqojjQFsjQAktV8VyjbYHApTQclUs12h7IEAJLVfFco22BwKU0HJVLNdoeyBACS1XxXKNtgcClNByVSzXaHsgQAktV8VyjbYHApTQclUs12h7IEAJLVfFco22BwKU0HJVLNdoeyBACS1XxXKNtgcClNByVSzXaHsgQAktV8VyjbYHApTQclUs12h7IEAJLVfFco22BwKU0HJVLNdoeyBACS1XxXKNtgcClNByVSzXaHsgQAktV8VyjbYHApTQclUs12h7IEAJLVfFco22BwKU0HJVLNdoeyBACS1XxXKNtgcClNByVSzXaHsgQAktV8VyjbYHApTQclUs12h7IEAJLVfFco22BwKU0HJVLNdoeyBACS1XxXKNtgcClNByVSzXaHsgQAktV8VyjbYHApTQclUs12h7IEAJLVfFco22BwKU0HJVLNdoeyBAi19Js9eQa5ZzjXYxE755BWjpC61EaA255jjXaJcy5ZtXgBa+zmrKl+NduaY412gXMumbV4C+8jr6fXa5pjjXaBcy6ZtXgL7yMvp9drmmONdolzHtm1eAvvI6+n12uaY412gX4htoOWugrck1xblGuxBroOXshW9NrjnONdql2AtfzHGgrck1y7lGuxjHgZbSoK3JVbFco+2BACW0XBXLNdoeCFBCy1WxXKPtgQAltFwVyzXaHjQWoD++rlar3/46/f2fP1a//vnE8zVoa9qqmP7Mpq0AXR8OQPhy/IEG7V1TFdOf6TQVoMf+PH/Ia9DetVQx/ZlPSwE6bB993/737dyhGrR3DVVMfybUUoCuV7/8a/eHbVseOlSD9q6hiunPhBoK0G0zntaW1ocO1aC9a6di+jOjhgL07993G0h7hw7VoK377Ky6diqmP9vz+jmdrQbo0KFfNGjzPr2uQzsV05/NmeCqIs0G6LDO9E2DNu7zK4u1UzH92ZoprmvXUICO15j2f11916Bt+/zatu1UTH+2ZoorKzcUoMPhIaOP+OGokV/+Q4O27LgF1UWAztCfrpk4p0nu7dFSgA7H2X0b/X2ze/t90KCrO2b711Hg87I0VLHJ+1PLzivbN9DdEcrjzaSNAG3boSCdBOjU/TnFEh0fSLYG+nPfoePP+OEz3yZ8ux74WGuqYtP2pzsXzi3XXvjBP39cNOjQsgK0Wfvlz44CdNL+nGSJjg+9vk3aWIC+SjOG8sB3rFwV8w20NQKU5TywCJWrYtZAW9NsgL4dr9zwFM0Yy+eLUK1WbIr+tN8zPgHKkj4NiFYrNkl/is/wBCihtVox/ZmDACW0ViumP3MQoITWasX0Zw4ClNBarZj+zEGAElqrFdOfOQhQQmu1YvozBwFKaK1WTH/m0GyAltGgrclVsVyj7YEAJbRcFcs12h4IUELLVbFco+2BACW0XBXLNdoeCFBCy1WxXKPtgQAltFwVyzXaHghQQstVsVyj7YEAJbRcFcs12h4IUELLVbFco+2BACW0XBXLNdoepAtQWrN0z9S09FzzvMl7YOoXnNLSk83zlu6Zmpaea543eQ9M/YJLyvX2rcOcvswUzm+xOe6qtjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgJ0Cjp1eub0ZaZwfgIUoDUCFKCQAAUoJEABCglQgEICFKCQAAUoJEABCglQgEICFKCQAAUoJEABCglQgEICFKCQAAUoJEABCglQgEICFKCQAAUoJEABCglQgEICFKCQAAUoJEABCglQgEICFKCQAAUoFD5A//59Nfj1z6Jnv60ufBt+9s8f5x98n/Tf2orp5/Tnzx9fR3/JZDwdX5b+x7Rq3z2zzd/25ed66eABuhl1Z0nY3XuzH+IjbYDOMafnD6V0M3o5HemGP4XRFM6Tc2kDdH3RnL/99fQL3HuzvxogjZtjTjN/qb+ajmzDf924eYoa8nNZA3R9/kx6e3ly30YvlbjLZ5nTt0OQDl/uZ3kHBDZqp2FDNNvwX7bLz1/+Nfxxt204R9IlDdC3i+lcr15aYtucentduvjXg1nmdPse2L8DhrfA4U9pjD+PEw7/ZRc9uJnnO3zOAB0+j8bv7rfVC915bu3tuz3vt4R55nTUn/m+3l+MeJ1u+K8avrWPOvJtlu/wOQP0ei6H7/o307Bf0fs8BNanMm3f9nl3ls40pxc/S5YgVwF6mrbN9ZLoxbRuJ367HfT6Ikrz1pcTcPhQPm/TDOF3fMBmNF23j9g+dfuft9Mcv53nfxeg1wsElwVab5+8eajtL8UN0OG9ffle3NwceXNewP+kC88b8MMfv+/mMuHG1lxzerKd2GzrI1eb8IcpORyXc57dq2kdAvT/HB+T7CNn7LYjjz++CdAfF9N1P0D/3x/HPZvHY22+HB7x5fD00ZMuCrQN0P987HvDpbgBuh3h1XvxZoVpvAP0w6W80Wxvn/TLf51zj19kc83p+MnZ0uByJ9JhSk5vz+Mb9HpaL3Y9J/woP7rtyJ07AXo++nD3m/sB+j+PD/mPi+Pqto/4L79fTPZNgdarf/t6CtwnxA3Q2y84159WpwW9T3f+vo0mZnwYT7YEnWtOd9Y5v01dHMZ0CINhVndv1M3hLXkzrfsA/f7zfARDTvc2Y37eDdD1fr6GCfxy9xG7hB0KsI/a4U9DTH75eYjL4THHetwUaPfyJXHQUoBeL7ptTm/hdz7HjsZfs05Tt5u7ZJ0705yeXil9gB6Gf57nwzTdTOvQhofpvfkkyuThAD3127urpOeP/M1q9KfhD8cgPW0k3BToei32YXED9E5f3dtrsfPJ0ttm9FKjN/4nEdGhmeZ057hRmu1b/dWB9IevM+N9SRfze5jWYbbOJ8Fmm7Ozmzbau43H8wP3ayb3A3T/+XV+jx9m+7y2cuj32wK9+z74RNwAvTO17wxy99b94M1+b7Xu+Hq5vjDNPqfHDaxERmugu0n7Nt6XdH0EzWla7+9mzufhb6DXx8fdDdBDw9786XJP/pd7BSrNgtABemdz8/vNgy5Xn+55t0Pfsm3Dzz+n+Q4lv3hnv+0mbXy1hfHW5GhaBejeO4O/nZ71VV89G6DHT7HzEU2XBbp+/UfFDdAH9hiPdqV99GZ/NyevN6+6V2FO030oXQTofjNyPIer8yqcAL31zgfurAH62193CtRfgH5+zOLxc+TXPz9cr3t3Cz5fgJrT6V0E6H6CbyPxZloF6N5NR+5XmWb+BnpnyvsL0M/Pmnlb3c7YHe83aMrzDmeeUwH6/c7U3UyrAD247sj9mvz8a6DXvd1hgF7vkLg+b3s0hbfn01w+79vdv73/Napbs89pvh1z15vw2/m7+Vp1O60C9ODqXPhj153z7O10BOf5QLAv9x7xcYAee3nXn3fOf+owQD+7ctDFNYA+eLNfTs1oLt85hKJrc89pykPDrncijb5w7r+R306rAD26OJHgdGz2+timx6M7z0cirU/HHV0+4uMAPTx4c79AfQbo/tDs/eTeu+zC+jALu9+9O/z9ZRvGfz2fXZftC+jPCnOabAv+3mFM+yXP3U/X+1m8mVYBenQ+J2s/bfu5ODbS+fohhzORLqby4hGfBOguQY/fH24K1GeAvnP19OOeu/Gl5fc/ubdPb3zM12C8By7Ze31n7jnN9qX+6kD6/Tt3PI3fft6ZVgF6cnlM0Tn5Dv7H79fnwl/fmefwiI8D9HQu/P0CdRqgF+/M28vWnrLg3w9rGvfe7Df9eX7NXIt1R3PM6bmbs+XnVYCOljLGb8+baRWgI6PP9FP7HBPu2+nD+tRj3+4/4pO98IeKvF+gLgP07h0kz2/p060gN6dv5vfe7Ndv6k3Ot/rJHHO6f1rCJBgH6HhG3y6b7GpaBeiFtzufvkOsbudotLWzuZrjy0d8dhjTrusvd32O/pfdBihAVAIUoJAABSgkQAEKCVCAQgIUoJAABSgkQAEKCVCAQgIUoJAABSgkQAEKCVCAQgIUoJAABSgkQAEKCVCAQgIUoJAABSgkQAEKCVCAQgIUoJAABSgkQAEKCVCAQgIUoJAABSgkQAEKCVCAQgIUoJAABSgkQAEKCVCAQgIUoJAABSgkQAEKCVCAQgIUoJAABSgkQAEKCVCAQgIUoJAABSgkQAEKCVCAQgIUoJAABSgkQAEKCVCAQgIUoJAABSgkQAEKCVCAQgIUoJAAncJ6tfp28YO31erLh8/454/Vr3/e/9Xf/3uqfxa8Y9ugI7/99cRT9eeYAJ3C5qoHt/H4y78+fMb7Abq+ymKY3mWAfvZxP6Y/LwjQKWzjcPV99PcfXz/7UH83QDfXX2ZhetcB+nDT6c9LAnQSV5vsb+VdpkGp4KJDtz337oLSNf15SYBOYvuVc9SBHyxwfkqDUsHlR/zfvz/cdfrzkgCdxOU2/OaZNaUrGpQKrraRPt3peaI/LwnQaVxk5vqcpuvdAtNxj9J69dtfm93fR19SN7uHfB/9efu3UZ9uvx5oWaZ1FaCj9h23480PTv1Z7R8anQCdxjblTlvt2z8fdiGdl+r3P9gG6H+uLgN0u/G/t/vrqUHPr7H92Sd79OFZVwG6PgbosDF/bsebHwjQawJ0IqNvnadvj+Ndnd/2D/q3r/uDRo4BesrPfYeeG3R9jM3tI586TA8+dxmgQxfu/nqKy2OCXv1AgF4ToBPZ9uBxI+iYfael+eEPuxBcn76LHgJ0WDvdPXhzPBbvmL7b1/v28+IPMJmLAH075eX68FH/dmzHmx9YA70kQCdyXtQ8HQR6Xlg67qRfnw64Pzz8fAT+NmSPSfrt8ID9b95swTO5m+NAd98pt316aLbDn25+IECvCNCpvI023K+3cI4rpOdTPg8Buj6n46EzR9v/u1/ZgmcG1wH6/fjTL+cHfLvzAwF6RYBO5fjF8/Yg0GE7/RSg308/2/5ktKvouAZwatDDprsteGbwdic+Lxby9+148wMBekWATubwbXK0GHpedD8F6PEL5ylAx4YwPTXo4avnW/kx+fCe8fbSKSPHl3DYfR+4+YEAvSZAJ7M5LrePe+7jAB0/4DpA94cvbR9XfEw+vOe8E2nowUOLjTeedhtHNz8QoNcE6GT2HTbaKj9+v/z1z/Ma6E2AXi9wnht0dwD99hEOGWFyo73wm9Voo8c30CcJ0OnsvnuOGuztdDjyewE6Pvz+YPT84bwlW/DMYXwY03kr/vJg5i93fiBArwjQ6ez296xPiTf6+N68swl/fRm8nxcNut2G/w9b8MxhHKBDF+6bdrTTfX29F35tL/w9AnQ6w26f//v7qePOATpsy98N0NGX1JvDmHZP+2+24JnDxYH0p2VQx4E+TYBOaLsN/79GXynXh3TcHTKy676bAN0tk34/PPjYoKdl0fMZIjCpy1M530ZdeDrx6HTu3MUPrm++kJ0AndAuDs/ttRnvYb8foBePOZ7CeVqSGu0ghSldBuhpI/6zc+Ev+xMBOq3j5/Xorzv/fljrvA3Q0ZFMpz2hd/eKwoSursZ0+qj+5GpMl/2JAJ3WZnWZePt0/HbahXknQI8nhZy+ae469PC3tc0lZnF905nznvgPrwf687I/EaCROY0TYhOggbkQE8QmQOO6OKseiEeABnWzFAWEI0CD2qys1UN0AjSo021qgLAEKEAhAQpQSIACFBKgAIUEKEAhAQpQSIACFBKgAIUEKEAhAQpQSIACFBKgAIUEKEAhAQpQSIACFBKgAIUEKEAhAQpQSIACFBKgAIUEKEAhAQpQSIACFBKgAIUEKEAhAQpQSIACFBKgAIUEKEAhAQpQSIACFBKgAIUEKEAhAQpQSIACFBKgAIUEKEAhAQpQSIACFBKgAIUEKEAhAQpQSIACFPr/1Vo60or7ECkAAAAASUVORK5CYII=)
#Nothing in the plots is shows an indication that the model is not
#reasonable. We doubt that the big residual deviance is because of
#overdispersion.
#Fit of model with overdispersion, quasi likelihood
fit2<-glm(resp~variety*root,family=quasibinomial,data=dat)
summary(fit2)
##
## Call:
## glm(formula = resp ~ variety * root, family = quasibinomial,
## data = dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.01617 -1.24398 0.05995 0.84695 2.12123
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5582 0.1720 -3.246 0.00475 **
## variety2 0.1459 0.3045 0.479 0.63789
## root2 1.3182 0.2422 5.444 4.38e-05 ***
## variety2:root2 -0.7781 0.4181 -1.861 0.08014 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 1.861832)
##
## Null deviance: 98.719 on 20 degrees of freedom
## Residual deviance: 33.278 on 17 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
#compare to fit without taking care of dispersion
summary(fit1)
##
## Call:
## glm(formula = resp ~ variety * root, family = binomial(link = logit),
## data = dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.01617 -1.24398 0.05995 0.84695 2.12123
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5582 0.1260 -4.429 9.46e-06 ***
## variety2 0.1459 0.2232 0.654 0.5132
## root2 1.3182 0.1775 7.428 1.10e-13 ***
## variety2:root2 -0.7781 0.3064 -2.539 0.0111 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 98.719 on 20 degrees of freedom
## Residual deviance: 33.278 on 17 degrees of freedom
## AIC: 117.87
##
## Number of Fisher Scoring iterations: 4
#Note that the standard errors shown in the summary output are bigger
#than without the overdispersion - multiplied with the overdispersion sqrt(1.8618)
#calculate dispersion parameter phi
phihat<-sum(rp^2)/fit1$df.res
phihat
## [1] 1.861832
#Model reduction
fit2<-glm(resp~variety*root,family=quasibinomial,data=dat)
drop1(fit2, test="F")
## Single term deletions
##
## Model:
## resp ~ variety * root
## Df Deviance F value Pr(>F)
## <none> 33.278
## variety:root 1 39.686 3.2736 0.08812 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit3<-glm(resp~variety+root,family=quasibinomial,data=dat)
drop1(fit3, test="F")
## Single term deletions
##
## Model:
## resp ~ variety + root
## Df Deviance F value Pr(>F)
## <none> 39.686
## variety 1 42.751 1.3902 0.2537
## root 1 96.175 25.6214 8.124e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit4<-glm(resp~root,family=quasibinomial,data=dat)
drop1(fit4, test="F")
## Single term deletions
##
## Model:
## resp ~ root
## Df Deviance F value Pr(>F)
## <none> 42.751
## root 1 98.719 24.874 8.176e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#confidence interval
par<-coef(fit4)
par
## (Intercept) root2
## -0.5121761 1.0574031
std<-sqrt(diag(vcov(fit4)))
std
## (Intercept) root2
## 0.1531186 0.2118211
par+std%o%c(lower=-1,upper=1)*qt(0.975,19)
## lower upper
## (Intercept) -0.8326570 -0.1916952
## root2 0.6140564 1.5007498
confint.default(fit4) # same as above but with quantile qnorm(0.975)
## 2.5 % 97.5 %
## (Intercept) -0.8122830 -0.2120691
## root2 0.6422414 1.4725649