Yesterday, I received an email from a phytools user with the following question:
“I've run across a problem entering starting values for phylosig
to estimate Pagel's \(\lambda\).
Could you tell me what the function maxLambda( )
does? Could you also give me a pointer about how to
enter starting values for \(\lambda\) estimation? I'm having some problems with convergence and I'd like to
find out if there are any 'attractors' among a set of starting values.”
I'll try to answer both of these questions, but to do so I will use some
new features
of phytools, so users following along should probably first update phytools from
GitHub using devtools.
library(phytools)
packageVersion("phytools")
## [1] '0.7.11'
I'm also going to use the files mammalHR.phy
&
mammalHR.csv
.
First off, what does the (internal) function maxLambda
do? This one's pretty simple. maxLambda
compute
largest value of \(\lambda\) for which the correlation matrix implied by the tree is not exactly singular.
What does that mean in simple terms? Well, remember that \(\lambda\) also corresponds to a tree transformation
in which internal branches are multiplied by \(\lambda\), while terminal edges are adjusted to maintain a
constant total height above the root for each tip. maxLambda
gives us the largest possible value of
\(\lambda\) such that all the terminal edges are non-zero and thus such that the implied correlation between
all pairs of species are < 1. In addition, \(\lambda\) > than this value implies negative terminal edge lengths.
Let's see what this means using our mammal tree:
mammal.tree<-read.tree("mammalHR.phy")
par(mfrow=c(1,3))
plotTree(mammal.tree,fsize=0.6,ftype="i",
lwd=1,mar=c(1.1,1.1,3.1,1.1))
mtext(expression(paste("a) ",lambda," = 1.0")),line=0,
adj=0)
max.lambda<-phytools:::maxLambda(mammal.tree)
max.lambda
## [1] 1.014493
plotTree(phytools:::lambdaTree(mammal.tree,max.lambda),
fsize=0.6,ftype="i",lwd=1,mar=c(1.1,1.1,3.1,1.1))
mtext(expression(paste("b) ",lambda," = maximum")),line=0,
adj=0)
plotTree(phytools:::lambdaTree(mammal.tree,max.lambda+0.05),
fsize=0.6,ftype="i",lwd=1,mar=c(1.1,1.1,3.1,1.1))
mtext(expression(paste("c) ",lambda," > maximum")),line=0,
adj=0)
![plot of chunk unnamed-chunk-2](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAA2sAAANrCAMAAADva/CiAAABC1BMVEUAAAAAADoAAGYAOjoAOmYAOpAAZmYAZpAAZrY6AAA6ADo6AGY6OgA6Ojo6OmY6OpA6ZmY6ZpA6ZrY6kJA6kLY6kNtmAABmADpmAGZmOgBmOjpmOmZmOpBmZgBmZjpmZmZmZpBmZrZmkJBmkLZmkNtmtrZmtttmtv+QOgCQOjqQOmaQOpCQZgCQZjqQZmaQZpCQkDqQkGaQkJCQkLaQkNuQtmaQttuQ27aQ29uQ2/+2ZgC2Zjq2Zma2ZpC2kDq2kGa2tma2tra225C22/+2/7a2/9u2///bkDrbkGbbtmbbtpDb25Db27bb29vb2//b/7bb/9vb////tmb/25D/27b/29v//7b//9v///9WpRH2AAAACXBIWXMAABM5AAATOQGPwlYBAAAgAElEQVR4nO2dAV/jVpZn5QzsQrpJbbvjJITqgDsOotJddrp3ZlIGJz0TqmQTaheQaiV9/0+yepIMsi1xn5B0dZ/0P/lVMGBbqsM7YChfnhUCADiw2j4BAHoCWgOAB7QGAA9oDQAe0BoAPKA1AHhAawDwgNYA4AGtAcADWgOAB7QGAA9oDQAe0BoAPKA1AHhAawDwgNYA4AGtAcADWgOAB7QGAA9oDQAe0BoAPKA1AHhAawDwgNYA4AGtAcADWgOAB6K15WfOS+5162a/H1vW3n++5I6EElwM5uvLL1RE4lqjRu63BaAohmgtuDjIvHJl7WlJu7E25C6tGKO8PE+2tZcpojFsIT0HFMVQjyFdy368rJrZfyDvMvjV2mjNO7TOVX9Py9N4sq29RFHfeF5RMO3SY55iqNaCiycvi+vw09dkMN6xtXecbW2ZfPZZmvVJ6Fk2WiuvqHc8ryi46NZ3GEXktfbp7R+jzz1fXMevLDOfkiLcjWAiSylPLi8H399fZFpbr0vvsDuf8KO/07urlyoKLj77EN14GAZvo/+rNz4Jd+NrRbex4wdIO1e9jF16hwc59yKN4OrQGkySy88qCv9f9BfZu87etJOKclqLHvPFrBM5yL7T3Xgtr7X/86D+ktnWklf8rxr6DrkFotb++GJFkZC/qNfPL9LvYrPCL9UbltZBuF5Im1fdWEib75JG+hdPzuxZRerK/zpcf+ZKbttFRTmtXVpfRqvi01fJ2flfZb8c6QWz0dr6DjbeaDjROhr8ZxhcWvGaKakouvFn19H3r9b+dbhSBWaF+18N5m787W6ykLavmllIW++SxlJ9qfIOk09HGopu/mj9aV1bNxUVf7+Wfn+18a1J9Eadn3D0oTX1mChdMyUVJTdWK2Zz1SXCXWvvq/jO04W0edWNhZR7L1JIpZRZRavDpy/9XVSU21rw8b+n366/6l5mvERn/F8aX4t70Fryd7lMvg0ppyhZd8ldrGvNCL+0kq+WyULaump2Ie3ciyg2v5JpKNr8utZFRTmtqe8kMw+1s5aizyv+OHvVnO/Xwq2suvn9WvJ3We62RitKbpxZApvC3fSTfvrNiFO4kLbeJYxnWttWFOZ9v9ZBRbnfrw1O//d//9/lbmvqhIMfM+es11oXfw6Zfl1L/mrlFO2ujqxw/6v0moYtpG2KW9tRlPtzyA4q2m0tlRStkZ1H2vGn8aW9c5NtNh8udvLf14q+X9NQtL0ENoVfWl9ept+V5C2k+Ljqh3DSFtI2qZT0B5DPKtr997VuKsprbRD9xT8dWzs/h0wuutFnoJvJs3e62Vonnzeifg756WLn55A6inYXUka4+gFb8t187kJaqk/ovx8KXEg7xD+H/JR8Unpe0e7zRrqpKPcxZEq8kDL/MpJ8xogfFD3/te3pW7T0Z1Hdez5k8i83yV+zpKKcz8SPwv34B2zxP9fmLqTk35n+l8AHSDukD56TxqAozP/ZyFV0tv/2/W3y/dXy8cvR2pd3bA2fv9Ot1rr4PP/4GQl/Sj5Xl1S0+43/k/DkX+zix0q5C0nd82Ai8ZuRXeLnjQyh6JEyz4cEuUARCRQpyjzPH+QCRSRQpCgzvwbygCISKIppZi67T0ARCRTF4PeNAMADWgOAB7QGAA9oDQAe0BoAPKA1AHhAawDwgNYA4AGtAcADWgOAB7QGAA9oDQAe0BoAPKA1AHhAawDwgNYA4AGtAcADWgOAB7QGAA9oDQAe0BoAPKA1AHhAawDwgNYA4AGtAcADWgOAB7QGAA9oDQAe0BoAPKA1AHgwuDX/2NpTO8MushsOXc2Lrt5L4KgYdjfmtha8OQ+9Vw+hf3p8ELwd2N7p598N3k0/XFgj/+vDkT8OZvcXVr/3joWjYvjdmNuaNwnjz0neq/vZr+fB7GbkjtwfxsuRf3Y7vBt7E3/sf3ve9lm2CxwVw+/GzNYsywpdO/rc9LeH6KX/+u+Wtf+v+Wq+/Otk4UTabG/o2pHMT29HbZ9qi8BRMdYaRjeGthb98V4/BL9GKlZzb6Ieci+chXP1D/tq7g1Xc3cUiYsu++O2T7VFLDh6hsSNxejG3NbCm0NrGCpH0aegw8G72f3s/sd/zL3DvUjayl4OTqPLgz7/GMCCo+eI3ViMbgxuDRDAEg2nIzM/HmaeNTewRIPWKMw8a25giQatUZh51tzAEg1aK8B6+kktKASWaFpxZNTHw9q5AHaBJZpWHBn18cAq0gGWaNAaBVaRDrBEg9YosIp0gCUatEaBVaQDLNGgNYpNQ+v5I8XV5oVeT2hlLcFRPspR7CaVxeLG3NbW80euHa7eTT9Mj88v9qMX778+HAXT98f9ndDKWIKjAqz1+rEY3RjammU9zR+pQaOxf3K9PPde3Y294d1YPTv7Y28ntDKtwVEBmX9f43NjZmuK9WyWP1ajRtGqWjh+5CZ6ZRj9//DowO3phFbGEhwVYLXhxtzW1rNZwT9nD64dmbua34xcW80eufbCuRn1dUIrYwmOCrDacGNua+vZrPAyesA9X8zD6HNRZErNHkWGrKNJXye0spbgKB/liN2NSa09/+S1YPaw6NV6KeA5S3AUk/s0yObddKe1cJl8muo7z1qCI4WVu+wbd9Oh1kAMLJHkt8ZwWHPAKtIBlkjQGglWkQ6wRILWSLCKdIAlErRGglWkAyyRoDUSrCIdYIkErZFgFekASyRojQSrSAdYIkFrJFuraGlZatgh+MXZuaZ6Xmlf2bAER3kkrbG7Mbi1SzsMfnQiZX8+tkbB9PiHdOYouPhs+nMPZ7JSNizBUR5Ja+xuzGotQxhcWNZgkuyj9fHs7uR6PXO0Ove/+a2HM1kpcERiha24Maq1jde8Vw+uegKba6tho8jYeuZo4fivb3o4k5WStQRHucSK+N1Ib23za1kW9yAM3szV3mJq2Ch6hL2eOVrMb4b9msmyiizB0SM7ivjdiG+t4HL0re0o+XN5EA8bzcP1zFH0+WnSr5ksq+gVOHpkZyHxuzG3NfBIYWvgEQELSfpHRoAi+aA1GgELSfpHRoAi+aA1GgELSfpHRoAi+aA1GgELSfpHRoAi+aA1GgELSfpHRoAi+aA1GgELSfpHRoAi+aA1GgELSfpHRoAi+aA1GgELSfpHRoAi+aA1GgELSfpH5jlF3h/W0xC93t/o2dbgKCVnIXG7Mbi14M3pPN5FaxRMnStr/8PbgR30dEak6BU4WrO7kNjdGNXa5tNH1R4H6s0fz+7Grnpy23kwu+vpjEj2FTjKY3shteDGpNY2UfNHkatkIEL9JvZLNWbb3xmRPODokR1H/G7MbW05isf8koGIq7maO1IPvPsyI5KlUBIcPbLjiN+Nsa35J06ovKhJiNXcO9xz1AxEf2ZEshRJgqMnth214MbY1sATkEQjwJGAU3jCyqHtcxIHJNHIdCTgFJ7IORlR5ycCSKKR6UjAKTwhU5EwIIlGpiMBp/CETEXCgCQamY4EnMITMhUJA5JoZDoScApPyFQkDEiikelIwCk8IVORMCCJRqYjAafwhExFwoAkGpmOBJzCEzIVCQOSaGQ6EnAKT5RTdGlZe8kTanq1vVEpSXBU/CZ2N+a2pn4Z+6ezhzDewufYGnmnn383UBcaPEEJlJEER8Vv4ndjWGuZp9z43zjqWdqh2sLHCT+e3YzigaSPZx2fg9RYR3BEvsmy+N2Y1trTRe9V9Dlp9hDGW/gcHh2s5qv58jS60Ow5to5Oa+sLcFT4JovfjbmtuZGUZfRAW23ho4aQFs7CufoiutDsObZOmdbgqPBNFr8bc1tbJttChvEWPtbR97P72f2PP1lHk0ZPsX3KtAZHhW+y+N2Y21pfKdNaX5HpSNSHBa1pIHMdyUKmI1EfFrSmgcx1JAuZjkR9WNCaBjLXkSxkOhL1YclTJG6SvW00JPXekkxHoj4s9MmIOt120FDQe0syHUn6sGh8rpF0uu2g8wm555a0vmihNfIqDKchG6HrSBIWvq6RoDUN0BoJWqMp2Vrw1hrkD0Nszkh0azOkcq310pFma+xuDG7t6jz8dKJ+B7u1//7YGvlfH47iDX4uPpu+ezuw1etqi5/76YcubYJUrrVeOtJsjd2Nua15w+SlGz9Z9OPZ7fBurF7+69z/5qfzYKZe9789D/34/52hVGv9dKTXGr8b01p7+ucRL32SqNrOR23s49reUL1UMxJ/t6z9m+j18NPbUXS9Lm2CpNVavx3ptGZZ/G4May2D//o6/H2sHk578UzEau6O4tmIeEYiDNXraosfteFPhzZBKiWpn460vq614Mbc1sKbQ2sYPYxW2/mojX0WzspWL9WMhNraR72uXq66tQlSOUm9dKTXGr8bg1vrJ5BEotkaO5LOCstIA0giQWs0WEYaQBIJWqPBMtIAkkjQGg2WkQaQRILWaLCMNIAkErRGg2WkASSRoDUaLCMNIIkErdFgGWkASSRojQbLSANIIkFrNGWXkfcHJ+etnRrF2qWkpD460m2N243BrQVvTiMt/no7HzVzFEyPz6fO1kBSQ2fbEuUk9dKRZmvsbmS1VgI1eJTMz6bb+RydB7O7k2t/vD2Q1JURyAQ4IrF0JLXgRlRrpa4dXFiWnYwexdv5RP72byehN9keSGrobFui3Je1XjrSU8TvxtzWlqN4FDKznU/8iyLUuNHmQFJDZ9sSpST105GeIn43xrbmnzihGuF73M7nNzVzNA/VuNHmQFJTp9sOpUZF++lIS1ELboxtra9AEolURZLOS9K5iAWSSKQqknReks5FLJBEIlWRpPOSdC5igSQSqYoknZekcxELJJFIVSTpvCSdi1ggiUSqIknnJelcxAJJJFIVSTovSeciFkgikapI0nlJOhexQBKJVEWSzkvSuYgFkkikKpJ0XmXPJZ4/CmZbz8N+mj9a5I0nmc5LZvx65khXEbcbg1tL54++PRz542B2r7bKikeQ4tmjk+ORf3p8sH5Pd2ZGXjTj1y9HmorY3RjX2vb8kTf0X99Okk2ylmoEaZzMHt3Pbl9Ffybd2lqsnKSeOtJRFNlhd2Nea+sL6/kj2z+7sZNNsi7jEaRk9iiY3dj+6/V7mjtpbkpJ6qcjvTVtsbsxt7V0/khNGa3my3iTrEU8gqRmjpa2N1zNvcn6Pc2dNDelJPXTkWZr7G6MbW09f/TrH/ec5eA03iQr2S4rnj36y56zcFx7/Z4Gz5qZMpJ66khvfs1id2Nsa8+y88Ol7lCbpO460vy69sz7mnHTzdY69IBom9okdddR9daacdPN1joMJJFUb60ZJH1UsIw0gCQStEZDncv6H416Dfn3hyQ9RWiNfL+kE24DeiHpXKnTaCnC1zX6/ZJOuA3QGglao0FrGqA1ErRGg9Y0QGskaI2mZGuXlrWn/kF/YxOf4JcODolkKddaLx1ptsbuxtzWgjfz8NPZQxhM34+D2V26gc/S+nPXRkQ2KdVaPx1ptsbuxrzW1vMi/jdOPM7nj71J9CfdwOd24nVtRGQTvdZ67UhHUWSH3Y1xrT3ivXqIn7fmTVw7/jNMhiHcro2IbKL3STuln470FPG7Mbc19yAMl/EAkhp8WG/gE2+l1a0RkU1KtdZPR3qK+N2Y29rSsgZqv57IzeB0vt7A5/KgcyMim5RqrZ+O9BTxuzGntb4/8yiFkABLchVJ+sgQrYk619agFhIsiVUk6SOD1jSQupAEIVWRpI8MWtNA6kIShFRFkj4yaE0DqQtJEFIVSfrIoDUNpC4kQUhVJOkjg9Y0kLqQBCFVkaSPDFrTQOpCEoRURZI+MmhNA6kLSRBSFUn6yKA1DaQuJEFIVSTpI1OyteCtNbC3J46uOvq8o0fKLaReOtJUxO7G4NauzsNPJ064tP6sNvFRU0fT4/Pp+6+f9vjp4nhWudZ66UhTEbsbo1rLEnrD5O1qAwQ1d+SP/ZPrZBCpy+NZ9ELqvSMtRS24Mam1TZQjhWsnc0feJHpLOojU4fEsciFl6acjPUX8bkxo7XHnvg3819fh72M1E5HMHS1t146HkZ72+Gn6hNug8AOWZ6mfjorXdFYRvxsjWst/782hNYweUV8eqLmjeAZpHotbdW+XowzFreW9u5eOnmkt+252Nwa31k/KtdZLNFtjR9LHBq1pgNZI0BoNWtMArZGgNRq0pgFaI0FrNGhNA7RGgtZo0JoGaI0ErdGgNQ3QGglao0FrGqA1ErRGg9Y0QGskaI2mbGveHzb27+n6qEhCydb66Ei3NW43BrcWvDmN7Hinn383ePd2YAfT98fWKH7155Pjzs6MlGutl440W2N3Y3Br7si14/+P3KPzYJbs6XOjXv3r8H5229GZkXKt9dKRZmvsbmS2Zm2S+7bgwrJs9TTt1XwZvb5/O1GzEfGrp3Ywu+nozEihpBxLPXWUWdN5jlpzI7S1grdnWY7iEaSFs3CuvnDUHJKajYhf/c72hl2dGSmSlGupn46s3Ivbr/K7MbY1/0RtCRlGD6pn9z/+djiIPz0dfR+/evqXPaerMyNlWuupI63WWnBjbGvPoHaL7Cylvq4V02VHel/XimnKjZzWsjPF1Vrr3IOiDIWSSlrqsCOramtNuRHUWm1f17pMXa11mMqtNYWcjxBa0wGtkaA1ErSmA1ojQWskaE0HtEaC1kjQmg5ojQStkaA1HdAaCVojQWs6oDUStEaC1nRAayRojaR0a8tR9MfefFvn57NKttZHR7qtsbsxuLVLOwzezMNAzR5dfPbP9/GuPh3f7qhsa310pNsauxsjWssjnomwPnPC5Xkw+/XcexUPHHV9u6PnWoOjhOdaa9WNCa3l4716iP9cWtb+1PHHya4+Hd/u6JnW8uilo2day8LvxtzW3IPkz8JRf9a7+nR8u6OSrfXSkWZr/G7MbU19axs95A7VkJH7+VG8kU/ntzsq2VovHWm2xu9GYmvbO/bp0PUfrq15UlNeUk8cZT8daTvicCOytfJnFUwd+kpdoEJrfXH0ktZY3HSktd5Q5etaT3jR1zUO5JwMWtMBrZGgNRK0pgNaI0FrJGhNB7RGgtZI0JoOaI0ErZGgNR3QGglaI0FrOqA1ErRGUro1/9jay/8XSHdjUqJT/4RbsrU+OtJtjd2Nua0Fb87j546GV9a+2s3HV9MQ0cvg4rOp2uNHvR5cWPv30w8dmh0p11ovHWm2xu/G3NbUxgfqiaNu/Hzsj2e3w2RXn3+d+9/8dB7M1OtqMqJbsyPlWuulI83W+N0Y1lpm/kh9oQ/+9hAuoi/xajefeBoierlw/Nd/t6z9m+h1NRnRrdkRrdb67UirNcvid2Naa08XvdcPwa8j9XDaGyYzEe4o3tVnfjNUcxLqdTUZ4XZqdkSvtfWFXjrS+7pm8bsxt7Xw5tAahmomYs9ZWkeTeBoiehl9bpqoiQj1uhfv9dOl2ZFyrfXSkWZr/G4Mbq2XlGytj+i2xo6cDwta0wGtkaA1ErSmA1ojQWskaE0HtEaC1kg2Wiuk1VMUQLY1WMplozVJjuR8WLKtPXOlnpNp7ZkrcZyJWPR+/RFaSy89c6Weg9ZI0BoJWtMBrZGgNRK0pgNaI0FrJGhNB7RGgtZISrdWNH/UpVmsHUq21kdHuq2xuzG3tfX8kZo58k4//26gho6C6fH51MlOJHVs76NyrfXSkWZr/G5Ma+3pn0fW80dq5uhm5I7co/Ngdndy7Y+zE0ld2/tIr7VeO9JqzbL43RjWWob1/JGaOVrNV/NltLj2byfR+tqYSOrY3kdarT3SS0ea24vwuxHVWql/0l/PH8UzR87CufrCiQWqeaOniaSu7X1U7okPvXRk6a0kfjeyWltf0iKdP1paR9/P7mf3P/6mho7moZo3ykwkdWzvo6f1o2Wpj44szUdI7G4Mbq2XlGytj+i2xo6cs0FrOqA1ErRGgtZ0QGskaI0EremA1kjQGgla0wGtkaA1ErSmA1ojQWskaE0HtEaC1kjQmg5ojQStkaA1HdAaCVojQWs6oDUStEZSujV3/yEMLkZhMHtIho4eL3SYkq310ZFua+xuDG5t9aeHcDWdh/44mH6I/jv+Ib7gNHeKAijZWh8d6bbG7kZ6a8W/3i+Ynjr+6YkTet+/mftj/+Q6udDKybOR2xocZclvTYAb8a0VXt8fr+aL/ziIHgv82/6DN/Em6YVmz7Jt8lsrunYvHRW0tn01fjfmtuZNVj//+3IUhsvzN3PXdu30QrNn2TblWuulI83W+N2Y25prryYPl7Ya9luOVmryKLnQ7Fm2TbnWeulIszV+N+a21k/KtdZLNFvjR8AppKA1HdAaCVojQWs6oDUStEaC1nRAayRojQSt6YDWSNAaCVrTAa2RoDUStKYDWiNBayRoTQe0RoLWSNCaDmiNBK2RlG/N+0POc7K7PS9SurUeOtJujduNyNY2no9dSPDmVHm5svY/JLv5nNrh6t3Tjj4nx6PojT9EF1nOn4dMazqS+ugo29pzjtjdSGxt863rC9szEaE7Uk8UdeNnj8a7+cT792R29LlXb1QXO7KxmCI3rCJJ/XSUaS371qd3t+XGnNa23xNcWJZ65mj0ueky2c0n8mRnd/QJZtEb1cVmT5wVorWtd/TSEdVa+pLfjbmtLUfxVn7JBj7xbj7BP2cPTzv6LOPds+Itf5o9cVbKtdZLR5qt8bsxtjX/xAnVEK3awMdLdvMJL+3sjj5/2XOiN6qLDZ85J6Va66cjvdZacGNsaxTxL2rpHuW+rhF005Hm1zWCBtx0trWO/lKNWlvrpqN6WmvATWdb6yi1ttZN6mmtAeR8dNCaDmiNBK2RoDUd0BoJWiNBazqgNRK0RoLWdEBrJL1rrehX8T5D/1p7gaS+tfYyQzIdNdbaC26Ru4yeWWLG84K/U35r3ZX0IkW5rbXvSHprFe5QPnW1VuUOhVNba1XusSYMbi24sPac6MUvzqJob5GruehRreZbM95Rg62xuzG3teDNf4beMAyX1p9Pjw8icSPv9PPvBu8urP376BX/68OR9O2PGm/NfEcvU6TjiN+NpNY0v/Snj7rXO4p4E+/V/ezXkX92M3JH7g/fnofL6JXb4d1YjSOVPg9GXrKQdB8fdcTRixRpOHr6zo3PjYmtJS9cW30uUi9d2389dYLZv+ar+dL+9Ha0iF65sb1hvNOPYBptTf3PfEdNtaacWMxuzG3NGz6oIYhwFRnyJsnU0cK5Opr7Y/XKaq7mbmVvf9R4a+Y7aq41b2gxuzG3tXB1OPjyIQwvDxaOa3uHex9m97P7H387HMQjSAtnZSuNpc+DkcZbM99Rc62FK4vZjcGtmU/zrRlPg62xK0JrLYLWSNBaE/crVlFzoDUStNbE/YpV1BxojQStNXG/lCIBT2irm/pb65ylBlprTZExrVW5a6nU3lq1+5ZI/a1Vuu9KoLUWQWskaK2J+xWrqDnQGglaa+J+xSpqDrRGgtaauN/SitKZiC2CXxzJMyIbNN+a8Y4abC2wmN2Y29p6JiK42J9+SPbuUXMQS+tLyTMiGzTemvmOmmsteGMxuzG3tfVTsJfn3qt07x41B6E29il97JZovDXzHTXXmjdJrsXnxsTWkn8bWc+LLBw/EpRu46O2FhE9I7JBo611w1FTrcVriNmNga2lrOdFruY3o3TvHjUHEf2RPCOyQZOtxZjvqMGva+xuzG1tPS8SfUKy07171BzE5YHoGZENGm/NfEcN/myE3Y3BrSUEs4eFKQtnm+ZbSzDYUZM/849hdGN8a9EnpOFLbiYBrtYMdtR4a4xuzG/NYNhaM5fmW+MDrbUIWiNBa03cr1hFzYHWSNBaE/crVlFzoDUStNbE/YpV1BxojQStNXG/YhU1B1ojQWtN3K9YRc2B1kjQWhP3K1ZRc6A1ErTWxP2WVeR/ZVnW/kPp4wii8dbMd9Rca/xuzG0tDL1XD/E05CiYHv91Eq5+Hgczo9YVw9c10x01+XWN242s1srsiRyG7kH0P7V/z93Jdbx9z+LUkAnIlJe11itHL2xNpBthrZW6xXIUqvGjYHY7CYPZyg7dUemjtsoLWytzfdMdvbQ1nStyuzG4teCNen622r/HtcNgevbgj02Z7E9pvjXjHTXYGrsbg1vzv1Fm4v17ImmXatzPmHHjhOZbM95Rg62xuzG4tQ2Cv5n0HX8Kw2PILCY6avIxZBYONx1pzf/alKH+LLytGemIqTUWNx1pzUyYv66ZCNfXNQ7QWougNRK01sT9ilXUHGiNBK01cb9iFTUHWiNBa03cr1hFzYHWSNBaE/crVlFzoDUStNbE/YpV1BxojQStNXG/1otuZTRNPom9I7xMkUxHBre2tNRv0YzHIEzZTWyL5lsz3lGDrbG7Mbe14M1cefK/PRwF0w8X1v79RbKDVukjt0bjrZnvqLnW+N3Ib61o/sj/yhrYarsR//Xd2P/2PJ5EUjtoVT53NmpsrauOamytdTc6f5cS04kbs3jlzyRXUcH1vVcPalc6144MTcJPb0dqEunG9tr5zfUvclRna7nX7oCjFxwm1FpI/G60WmvofndvUaY19yDeglXtmeXaV3N/rCaR1Gulj1wHL3p80Hhr5jtqrjV+N+a2dmlZe9G3tL/+cc9Zzb3DwTyeRHJW7TyZXWZr5jtqrjV+N+a2JguZrclCVmv8oLV6QGs0aK2e61S/jVhFWqA1GrRWz3Wq30asIi3QGg1aq+c61W8TK6rnp74twNmaqZY4W5OoSFZrdR2dH8bW6js+M4yt1XX4WkFr9YDWaNBaPdepfhuxirRAazRorZ7rVL+NWEVaoDUatFbPdarfprQi99D6MvPqcvDdoL2xEaGtGe+oQUXsbsxtLfjRCRbXYXhl7b8/tkZLy7L+59uBHUyP23gGkszWzHfUnCJ+Nwa0VvQz7psvIlPp3iIfz+7G/nh5HszuTq5fcLqVab21jjpqTFELbuS3VkxwFXlS+x64h0cH3sSbXFrW/p0DXCkAACAASURBVG07W0O03VoRpjtqUhG3G3Nbcyfhp9cP8Z4+C+dmtLSX9sJR80ilj1wHMlsz31FzivjdmNuaf2Gp72XVJMTSOpos5ot4MmLRzvf+Mlsz31FzivjdmNuaLGS2JgtZrfGD1uoBrdGgtXquU/02YhVpgdZo0Fo916l+G7GKtEBrNGitnutUv41YRVqgNRq0Vs91qt9GrCIt0BoNWqvnOtVvI1aRFmiNBq3Vc53qtxGrSAu0RoPW6rlO9duIVaQFWqNBa/Vcp/ptxCrSAq3RoLV6rlP9NqVvEVxYe07eO3J2+Eme5XY1b2zzH6GtGe+oQUXsbsxtLXjzn/EvZF8PIHmnn383eBdPIH1Qu/ucHI/86M3B9Pj84rPpz/HFD/HmPw+lz63+s3/pjcrdxHxHzSnidyOotXKEXjr8sB5Auhm5I/dITSCNk9197mcPai7p5Hp17n/jpCNKavOfJmBrrWeOGlPUghs5rZXFtaMHAZkBpNV8NVfDtfu3k2R3n+hPPJcULhz/9U06oqQ2/2nibLhaK4f5jppTxO/G3Na84UN4aYePA0gLZ+FcfRFPIKk3Le3kzdEj7cX8Jr2YbP7TxNnIbM18R80p4ndjbmvh6nDw5cPjANL3s/vZ/Y+/HQ7mq2R3n7+s55LUZ62JuriaJ5v/NHEyMlsz31GDitjdGNzas8RbjjMitLVnMcFRW4qacNPV1pp5FFSMia2Z4KgtRU246Wpr3JjYGjd9X0dorR7QGk3f1xFaqwe0RtP3dYTW6gGt0fR9HaG1ekBrNH1fR2itHtAaTd/XEVqrB7RG0/d1ZHBrwVtrkPxC6LzfC/34tsZGRDYQ2prxjhpUxO7G4NauzsNPJ04YXHw2facGH7JjEY9TEMfn0/fjYHbf0JjII0JbM95Rg4rY3ZjbWjx7FKEmHt6fXKs5iMexiH+tpyBOrv2xN2luTOQRma2Z76jJ5x4nL/nciGut9GyWmni4VYMPkaL1WMT0cQoiupprNzcm8nTabDdKbtkXR40pasGNvNZ0r+i/vg5/H4fxxEP04DoZjUjHIp6mIKLH3ZG/xsZEHuFuTe9q5jtqThG/G3NbC28OreFDPOenBh+8w70Pj2MR6ymI6O2RqcFpY2Mij8hszXxHDSpid2Nwa6IQ2poohLXGDlqrB7RGg9bquU4dt6l2w3ZBazRorZ7r1HGbajdsF7RGg9bquc7ObV7MCw4mgJe11i9LvOtInqKmWnsx8hRpwXzaRlriPWl5itBaPaA1GrRWz3VqQ54iLdAaDVqr5zq1IU+RFmiNBq3Vc53akKdIC7RGg9bquU5t6B/sMmfqiGsUawehrfXYEXk0djfPnFE7Pz3VPljwRlkJLvbjAaO7Y2vkf304CqbOlbWvNvUJ3g5yZdZKWz9h1jxcnx1RR+N381xrGtepH+1/PVETRmG4PPde3U7Uk7A/nt0O1a49ao8fNYz06znD79Bux5H2TE2fHclzI641bZJZv4UT2bG9idrVJ3o59NRzs+NhpL9bDY8ZK+CIRqgjfjfmtraMh/eu5jcjNWCkpo1Wc3eUbOiTDCMxnAQc0Qh1xO/G3NYuLesgnj6y1YBRPG3krOxkQ5/kT8PjWAo4ohHqiN+Nua3FRA+pF+38UC0BjmgkO2J1Y3hr4dIatnl4OKIR7YjTjemttQwc0cBRClqrBBzRwFEKWqsEHNHAUQpaqwQc0cBRClqrBBzRwFEKWqsEHNHAUQpaqwQc0cBRClqrBBzRwFGKwa1dWpY1CrefjL1+GtvVXP3X9DnAEY1UR+xuzG0tmT8K/W8PR77aIWt6bF9Z++9Pjw/U5FEw/RD91/hTa+GIRqgjfjcGtbY9f/SVZX3mhN5Qbenjj/2TazV55L26n6nJIzWJ1PC2K/E57VxoGzgiacuNSa1tvuq9ir/4u7Z/dmN7E28SqmeRRq++VpNHt5P4TWznJMTR1onAUR7JefC7Mbc1N9l8Tg0dqQGkdAut1dybqIfcrq3+4zsnIY62TgSO8kjOg9+Nua2pb233H8Jf/7jnqAGkZAstZ+G4tpo8WsX/8Z2TEEdbJwJHeSTnwe/G3NZEIM2RoBN5RJyjts4DrVVCmiNBJ/KIOEdojUTIaWwgzZGgE3lEnCO0RiLkNDaQ5kjQiTwizhFaI7AEbqglzVH8L0dtn8MO8hy1dWSNd4lwZMk4jS1kOVJnIeNEsohzhNYI0JoGaI0GrZGgNQ3QGg1aI0FrGqA1GrRGsqsomYlY8zgbwTEmkjmr7QutstMaHO3wuJDY3ZjbWjoTEapBCLWhj5qNiF6qYQi1yc/9BcO2EMIc7bQGR7usFxK/G5Nay52pSZ5Cqjb08V8/RC/VMISnZiS+PWc5q+0LrWKFuTM1cJTBCttyY1BrW6QzEfEgRLyhj392E71UoxBqk5/w09vR7r00dlYiHO2cBRztsj4LfjfmtpbORCTb9yQb+qiXahRCzUhczRnGIIU52jkLONplfRb8bsxtLZ2JUIMQakOfeDYieqlGIdSMBMt2R8Ic5f74CI42WZ8Fv5sKrWlsIlsnNf/F60GYI5GS4Gj9F9V4V5Gjus/lWUQuI2GOZEqCI/rIshyJXEbCHMmUBEf0kWU5ErmMhDmSKQmO6CPLciRyGQlzJFMSHNFHluVI5DIS5kimJDiijyzLkchlJMyRTElwRB9ZliORy0iYI5mS4Ig+sixHIpeRMEcyJcERfWRZjkQuI2GOZEqCI/rIshztHi2dP1pkNxNhHMuKkeVo93BwVHw4djfmtrbe0yfdxMc7/fy7wbvphwtr5H+dbvPT/HiWLEc7h4Oj4sPxuzGnte1nta339Ek38bkZuSP3h/Fy5J/dDu+4xrNkOdqWBEd5h2vNjTGt7fC0p0+8ic+/1POz/zpZOJE22xsyjWfBEY1QR/xuzG3tcU+fZBOfhbNwrv5hJ/v6JNv8MIxnwRGNUEf8bsxtLZ0/SjfxeTe7n93/+I95sq/PymYaz4IjGqGO+N2Y25oI4IgGjlLQWiXgiAaOUtBaJeCIBo5S0Fol4IgGjlLQWiXgiAaOUtBaJeCIBo5S0Fol4IgGjlLQWiXgiAaOUrRak/aL9uQARzRwlKLT2gvf3zSXlrU3Tzfz2ZmFYBocgSMaqY7Y3ZjbmpqJ+HT2EKrNfIKpc2Xtq3mIYHr8w9fqDR949zt64fsbBo6K4Xdjbmv+N/EzRkNv6L++G6tnkqp5iLuTa2+o9vVh3u/ohe9vGDgqht+N4NaKfyN7/CBfzUSor/9qM5/bidraR81D3E6iN3hD9v2OXvj+Gs4AjsgTEOJGcmvPv9s9iD4T2Woowh258SyE+uPaj29g3e/ohe9v+AzgqPgA/G7MbW1pWYNJ9FJt5rOKZyHieYh5qEYiVtz7Hb3w/Q2fARwVH4DfjbmtiaBtRyZYat2RFEVorRJtOzLBUuuOpChCa5Vo25EJllp3JEURWqtE245MsNS6IymK6mrt+R/Qv4wa/npNU24dNSFJvqXW15EURbW1VsO5lDy+BEq21sIZtE/b60iMIrRWCbRG0/Y6EqMIrVUCrdG0vY7EKEJrlUBrNG2vIzGK0Fol0BpN2+tIjCKDWwveWgM7vcy9z9Ea6a0Z54hREbsbg1u7Og8/nThhcGHt309/U1v6/DwOZrfqAtuqkt6acY4YFbG7kdza8/9m4g2Tq6lBI3/sqi19/jpJL9jP33ON51jq/Y0spG45YlPUghvBrRF4k/TCp7cjb7JSW/qc2ukFmeuIX5J5jvgU8bsxtzX/9XX4+1g91I4+Gdnxlj6nkaX4gszHR/ySzHPEp4jfjbmthTeH1vAh+vx0OJivki19fhqcJhcc+tb1ILw18xwxKmJ3U6U1ic85Y6aco35KwjpKqdDay6/cHcr9tfspCesoBa1VAq3RYB2loLVKoDUarKMUtFYJtEaDdZSC1iqB1miwjlLQWiXQGg3WUQpaqwRao8E6SkFrlUBrNFhHKWitEmiNBusoxdzW/K8sy9p/iLfPUqRPYlO/iV3Oc/2qXLsGzHPE+HxIdjfmthbGO43E22f5x2rjrA/T4/OLz6bvogtinutX5dq1YJojTkXcbmprjZsw2Wkk3j7rIfx4djf2T65X5/4378ccm6+s/9rlrg1H5JU77Kau1tpgqTbIUttn3RweHXgTbxIuHP/17eRxMql54IhGrCNuNwa3pjZhTbbPWjg3I9d27XAxvxnGF7iAIxqpjtjdGNxavAlrvH3W0jqarOaLeegeHsWDtWznAEc0Uh2xuzG4NQnAEQ0cpaC1SsARDRyloLVKwBENHKWgtUrAEQ0cpaC1SsARDRyloLVKwBENHKWgtUrAEQ0cpaC1SsARDRyloLVKwBENHKWY21pwYVl7TntbHcXAEY1QR/xuzG1NPT3btcNg+uHtwI7Ejdo4CTiiEeqI342JrSUzEd6rhyD6nOSPl+fB7NeRf/bQxqm0cExN4Og5ksEaXjdGthb/341cDdXWPpeWtT91HudrWzgVkcQTWnBUCL8bc1tbJl/01XY+6kH3euu6Nk5FJOrU4KiYyE386YjPjbmtXSZTRqu52tTHO9xjm+nfPRWRqFODo2IiN+rMGN2Y25oIBJ3KNmJOTcyJ7MB8ZmitEoJOZRsxpybmRHZAaySCTkzQqWwj5tTEnMgOaI1E0IkJOpVtxJyamBPZAa1RSNrnVdCpbCHHkpgT2YZbEVqrhKBT2UKOJTEnsg1aI5GziuQ6kmRJzIlsg9ZI5KwiuY4kWRJzItugNRI5q0iuI0mWxJzINmiNRM4qkutIkiUxJ7INWiNJFaXzR8EvzuMEEv+YllRHqSU4KobfjbmtpfNHS+vL6YeL/XhHn+iCtc/6VHapjlJLcFSM/9qyeN0Y2Vp2fs2bqAkk79XdybU/9r89Zz4X3sOVAI6e5XF+jdGNia0lL9L5I9f2JgvHH3uTSFv46S3v6LFUR+mZwVEx/G7MbS2dzVrN1a6rak+fePtVxj38MucikPjM4KiY5Sg+MUY35raWzmZdHkS2Do/sxTwdRWrjXAQSnxkcFXNpJ4r43JjQWs72qxsEs4dFS78oSowjyhIcKfIVcbkxorVnX41YWq0M94eCHJGW4CjcPpXH15jcdKK19hB0LmItyTmTotZaObruu3gRu4pEnYtYS3LOBK2RiF1Fos5FrCU5Z4LWSMSuIlHnItaSnDNBayRiV5GocxFrSc6ZoDUSsatI1LmItSTnTNAaidhVJOpcxFqScyZojSR/Fa339NnAtVnOaPtcJJBnCY42yZ4KvxtzW0vnRfyT45F/bI2C6fH5xWfTn8fB7J5vakSMo3xLcLRJ9lT81xazG3NbW8+LDO9nD+HHs7uT69W5/81vE9apETGO8i3B0SbZU/FeWcxuTGwteS7b40xEMLs5PDrwJuHC8V/f2KxTI2Ic5VqCoy2eTiVZRLxuDGwtJZ0XWdrecOHEMxGL+c1wNV9yTo2IcZR/KnC0SfZU+N2029r2c9MLyL1tOi+y+Mues7SOJot56B4eTZaDU86pEY51pCkp71TgqNgRv5uWW6t8rXa2ynyEZR1VvRocFV+N043prTHPGG8jeh2tgaPiq3G6Mb21lhG9joQAR/TR0RoN1hENHNFHR2s0WEc0cEQfHa3RYB3RwBF9dLRGg3VEA0f00dEaDdYRDRzRR0drNFhHNHBEHx2t0WAd0cARfXTxrXl/UONH8a/RXP/zf/DL9kRSs4hfR3BUfDVuNwa3Frw5jTQF/7wbB7O7bw/jKaSl9WfW4Syx6ygFjoqvxu7GiNbyn0Xqqudmh2qHkejPUI1Ffjy7nXisw1mS1hEc0Vdr1Y0JreWjZtgjV95E7erj2v6ZmkJybZd1OEvQOsoDjorhd2Nua8uR2jArdG01crSauyM1hbSa8w5nCV1Ha+CoGH43xrbmnzjxk7QjN4PT+a9/TKaQLg94h7OErqMUOCqmBTfGtpYSTHl/qLaFzHW0BRwVw+nmZa3pjsGWmZN9EcFFWzsdJXA4qioJjophdfPC1po/uhlwPCwwXRIcpaC1SmAd0cBRClqrBNYRDRyloLVKYB3RwFEKWqsE1hENHKWgtUpgHdHAUQpaqwTWEQ0cpaC1SmAd0cBRisGtBW+tQf52WXwTWtLXERwVw+7G4NauzsNPJ5GUK2v/g9o8yxqpi+/HwexX68/HFsuz2KWuozVwVAy7G6Na23hSjpc+u8YdqSeSXi9H/tnNSM1IqImkMPx4xjEJKXEdwZHG9dtwY1Zr2VeUEIWaYY83zwpm07makVDTSO7h0UFN56h/Rtrvqu0Y1PXhqPj6/G7Mbc1/fR3+Po4eAsy9oWvHL9SfZBpJzSLVdI76Z6T9rtqOQV0fjoqvz+/G3NbCm0NrGH2d9w73nOhzk3qh/qj5o8sDNYtU0zmWOCPdd9V2DPL6cFR8fXY3BrcmAanrSBJwRB8UrdFgHdHAEX1QtEaDdUQDR/RB0RoN1hENHNEH5WiNZ9K9OVjWkeGW4Gh9knW/q0mkrSKJjuRZgiP6oPIcSVtFEh3JswRH9EHlOZK2iiQ6kmcJjuiDynMkbRVJdCTPEhzRB5XnaPOwy1H0xw7DRYu/aVSeo80Dw1HxgdndGNzapR0Gb+ahf3p84E3C1c/H1ih4O7D9r5PtfVo4I913NUv2wHBUfGB2N2a1lv2hrdpnxPrMCb1X97O7sfrV7B/Pfj0PZrfDuzHXuIhAR1lLcFR44DbcGNXaBt6rh/iPa/uv72ereAzi0rL2b2xvyDUuAkcaCHXE78bc1tyD5M9q7k2C6dmDGoNQD73X2/uwnAQc0Qh1xO/G3NbUt7aX8be2rq0uqDEItZnPwlnZXOMicKSBUEf8bsxtLUvwN56tn3eAIxoDHPG46URr/tf5vw+peeCIRr4jJjedaK094IgGjlLQWiXgiAaOUtBaJeCIBo5S0Fol4IgGjlLQWiXgiAaOUtBaJeCIBo5S0Fol4IgGjlLQWiXgiAaOUsxtzX/9ELqTMFwO5qH6TdGtnAQc0Qh1xO/G3NaCH51g6oRL63+Mg9n99L/sZASJ9yTgiEaoI343BrRW9EvHFo7aakTt4OOPk/+zjWQ9nduL3tUEcETSthsTWit4u9pNRG3no3bxUab4djl6QoqjosPBUfHR2N0Y3Jr7VzX3oHbwWdquHfxz9sA2kvWIFEdFh4Oj4qOxuzG4Ne8P6peyLOZqF5/VfD2CxHlmchwVHQ6Oio/G7sbg1iQgxZFkS2Icta0IrVVCiiPJlsQ4alsRWquEFEeSLYlx1LYitFYJKY4kWxLjqG1FaK0SUhxJtiTGUduK0FolpDiSbEmMo7YVobVKSHEk2ZIYR20rQmuVkOJIsiUxjtpWhNYqIcWRZEtiHLWtyODW1C+IXqbPqdEciQh+qXn/HymOCp+jBUeFR2N3Y3Br6Z4+4ZW1/2H6/uvDeEefYHp8rl5cWPHGPurVi/3ph3Svn6X1pbqkhigurP0antItxVHR4eCo+GjsbupoLX/opT7yzyHeO+sbJ3RH8WDE8G68PA9mdyfX6sWvI//sIfx4Fr/qvVrv9eNN1Dtv1RDFt+clTcl2lG8JjooV8buppbXSRy1FgTflSe3ps5ingxFDtaPP7SRUL6ZOMLs5PFKb2C0cf/y414+dXpqEn97W8ZRuKY7yLcHR1v2368aI1nLfHHkKLkfqobZykOzkEz0It+M9WdUb1YRE9OrV/Gmvn2RkSQ1RXM3VaGD1c9N9F0NrOW+Eo837zx6B381zInS/GrfTWnBlWcPoIbN3uOes4p184h195qF6od6oJiSiV93DI3u918/lgbqkhijUy1rOTYajAktwtHn/2SPwu3nZ309Ca5oEs4dFK7/TRkJrmvTEkfWSI9TopvOthUtrWNOJlMOg1nri6EWt1eim+621hUmttYUJrdV5/Oq3Qmt5oDUatFb2VmgtD7RGg9bK3gqt5YHWaNBa2VuhtTzQGg1aK3srtJYHWqNBa2VvhdbyQGs0aK3srdBaHmiNBq2VvVVLrQVvrYHd7KGrIKI1ONq4/40jsLsxuLWr8/DTiROqMaN4oiiZPKpx8KoaIlqDo+z9bzpid2NAa+v733rWqvf41JmPZ/FEUTJ5VOPgVTWYW8uzBEe5x2vLjUmtbb7VS7c5UBv5JNv6hDUPXlWjrdayb4OjvOMl/+d3Y25r/uvr8Pex2rLuZpRs6xPWPHhVDQmtwVHe8ZL/87sxt7Xw5jCeP1JjRmqiKJ08qnHwqhoSWoOjvOOlx2F3Y3BrwhHRmnBabY0dtNYUaI0GrZW9FVrLA63RoLWyt0JreaA1GrRW9lZoLQ+0RoPWyt4KreWB1mjQWtlbobU80BoNWit7K7SWB1qjQWtlb4XW8kBrNGit7K3QWh5ojQatlb1V463l/ULqeE+flHj3rNr3DasI+zrKswRHm0fIKmJ3Y0Jr+YeJ984K/ZPjUTB1rqz9n6w/SxnLimnF0daB4Kj4qPxuzGstu6dPGHrD+9ndWO2hpaYihIxlxbTaGhwRR1V6uN0Y2Fr6ajLr59pqvE/tbuDarpSxrJh2W0tegaPio3rxC0435raWbHW8jHegU9tlreZixrJiJLQGR8VHXcYvON2Y29qlZR2E4eIvavcstV3W5YGYsawYCa3BUfFR+d2Y21pMMJPxTX4OElqLgaP8o6oXrG4Mb03Kg6EcxLQGR/lHVS9Y3RjemmDEtCaY1ltjBa01BVqjQWtlb4XW8kBrNGit7K3QWh5ojQatlb0VWssDrdGgtbK34mot/ynIUmmrNZMktdRaS44Maq2Vo72cNh2ZIqlfjgxubb2nj8hpESHrCI6eOTi3G4NbS/f0Cabvj63R0vpyevyD2t7n68NRIGFoRERrcPTMwbndGNja1n5H6p/+4119Tq7jC8M7GUMjLbcGR88f+fEbNj43JraWvFzv6eNN1rv6rLf3GYoYGmm7NfUCjorx0oXE58bc1tZ7+rh2sqtPvKdPfGEkYmhEQmtwVIwfLSReN+a2tt7TZzVXu/pcHizm8fY+C2dlixgakdAaHD13cIvZjcGtCUdEa8JpuTXe46G1xkBrNGit7K3QWh5ojQatlb0VWssDrdGgtbK3Qmt5tP19P+8BXwZaK3urrj+P7WW06sgQS+06MuN4aI0GrdGgtbK3Qmt5oDUatFb2VmgtD7RGg9bK3gqt5YHWaNBa2Vu11NrGbJY8RLQGR8UEFrMbg1t7ms0aB7O7k+ORj/2Otg8IR8VcWcxuTGxte35N7eWj9va5xX5HmQPC0bNHzsyvcbl5YWub1HxOemRms2z1J5jdyNrvCI5I2nTE76aOv1876+hpNkvt5aP29pG139EmcETD64jfjbmtZWazBqdztbePrP2ONoEjGmZH7G4Mbi2L4H2PEuCIpkVHLG460prMB0UZ4IimRUcsbjrSmnjgiKbrjtAaD3BE03VHaI0HOKLpuiO0xgMc0XTdEVrjAY5ouu4IrfEARzRdd4TWeIAjmq47Mrg1d/8hDC6yz14TOjmigCMaZkfsbgxubfWnh3A1nYfe6effDd69HdjB9IP6v5iJkSxwRMPsiN2Nua0F01PHPz1xQncU/Xd0Hszuxkv1fzETI1ngiIbXEb8bSa1ZZQj98Wq++I8D9eRRtf2BZe3fTi6j/z+ImRjJUts6giP6foS6EdVaqWt7k9XP/76MtCychXP1hZPs7qMedEt84l99rZW5cj8d6d0PvxtzW3Pt1eTh0g6D2f3s/sffDgfRpyc1ESFzYqSd1vrpSO9++N2Y25pZtNOaWXR9HaE1HtAaTdfXEVrjAa3RdH0doTUe0BpN19dRC60983PY7lL279ZHSV1fR220VtP9GEXp1uq6I4Po+jpCazygNZquryO0xgNao+n6OkJrPKA1mq6vI7TGA1qj6fo6Mri13fkjRfyUtne/OMLmtFpqrcuOqipid2NwazvzR77a1+f0+CCY/mR9Of0gakarpda67KiqInY3olorQ8780a3a1+fV/ewu3uVHzs5HivpagyPq6lLdSGqtHLvzRzfxvj7+69t4lx85Ox8pamutFJ12VFERvxtzW9udP0r29VF7aa3mrrCdj9pprdOOKirid2Nua7vzR2pfn4WjTF0erITNaLXTWqcdVVTE78bc1nYQvedRO63t0CVHdStq3E2HWpP0cGgHIa11yVHdihp306HWRCOkNdF0fR2hNR7QGk3X11EzrRH/tNFDyjrqo6Sur6OGWmv6kMZR0lEvJXV9HaE1HtAaTdfXEVrjAa3RdH0doTUe0BpN19cRWuMBrdF0fR2Z21pwYVl7jugNxTK001qnHVVUxO/G3Nb81w+ha4eBGjQahVfWvnoZTI/teu6+ZtpprdOOqj7Pn92NUa1t/OOK9+ohiD4n+ePlyD+7GYWhenl3cv3ie28Uxtb64qjK0m3FjVmtZV9xo8U0VKMRCyeYTedqNiKY3U5efOfNwtna08VOO6rUWhtuzG1tmUzyqUEjb7j+48p8dNRWa512VLE1fjfmtnaZaFGDRnvO+s9C6g8B2mmt044qtsbvxtzWzKKd1syCtzV+0BoPaI0Grb3kPtDaNmiNBq295D7Q2jZojQatveQ+nnf0cmo42ZYo3VoPLfGto3YctdBanQcyBjZHBlvquiO0xkPX11EddN0RWuOh6+uoDrruCK3x0PV1VAddd2Rwa8vR4xNt1sidHWlpHcFR8aHY3Rjc2uVnjnqijf/14cg/jschzqfOlbX/PnqlmeNXoKXW4Kj4UOxuzGot+0Pb4Mf/GvnfOKE3vBuH4Uc1DuGP3djSxzNxv0mbcx3BEX2oNtwY1doG/uv7H397peb9vKF7eHTgTdR8RPQgQL3CcPxywBENryN+N7W01so/FboHoXsab5/ljhbOzci10/kI9QrD8csBRzS8jvjdNPO34VhH0be18W7HC2dlL60j9SkpmY9QrzAcvypwnwZv5QAABDRJREFURNOkI3435rZmOnBE0y1HaK0t4IimW47QWlvAEU23HKG1toAjmm45QmttAUc03XKE1toCjmi65QittQUc0XTLEVprCzii6ZYjtNYWcETTLUdorS3giKZbjgxuzT3YmT+KkTugtQEc0TTqiN2Nwa1d2mHwZh6q2SPv9PPvBu/eDuyNIaR4MmkczO4vrH1xAyRwpEGjjtjdGNXa5vya8vSNE6pxo5uRO3KPzoPZxhDS7fBu7E38sf/teTPnUwmW+TU4KrjjNtyY1Vr2lXjOT80fHR4drOar+TJaXPu32SGkeDLJ9ibhp7fyBkhY5rLhqPiO+d2Y25raq+4ynom4GS2chXP1haN2+MkMIanJpMhh9Lo/buaEqsDRGhwV3zG/G3NbC64saxg9jF5aR9/P7mf3P/52OJhvDCHFk0mD0+j1gcAfBXC0BkfFd8zvxtzWTKdbv4+tGbrlCK21RbfWUTN0yxFaa4turaNm6JYjtNYW3VpHzdAtR2itLbq1jpqhW47QWlt0ax01Q7ccobW26NY6aoZuOUJrbdGtddQM3XLUUGsN0cjJtkRTjrpkqVuOjPrAbJ7s7kxE8IvDej4yyVqCo3yUI3Y3BreWmYkIpsf2lbX/k/XlRfSKyPEQPrKW4Cgf5YjdjbmtZWci7k6u1SCEN1mO/LM7keMhfGSfDwlH+VhtuDGxteQBd3YmwpuE8SCEvXCC2YPI8RA+nizBURHpd228boxsLf5/dibCtUM1CLGaq//LHA/hI2MJjgpIZmosXjfmtpaZiVDTfWoQ4vJA/V/meAgfGUtwVIDVhhtzWwMFwBKNlfk/80ENAatIB1iiQWsUWEU6wBINWqPAKtIBlmjQGkX3noTUBLBE04oiEz8iJp4zP7BEgtZITDxnfmCJBK2RmHjO/MASCVojMfGc+YElErRGYuI58wNLJGiNJD3nvD19rt794hiy31HTxJbgqBj3wOJ1Y3Br6fxRvJ/Pz8kQ0vn0J+vL6W/xW8RudMREbAmOirm0LV43RrYW/9PIev7IVfv52MkQkh/v4RO/5a9iNzpiAo6eI/n3NV43JraWsJ7NivfzOU2GkLyJ2sMneYvYjY4YgaNi+N2Y29rTbFa8n088hOTaq7ka+Ivecip2oyNG4KgYfjfmtpbOHwXxfj4/JUNIq/nlwepd8haxGx0xAkfF8LsxtzUAzAKtAcADWgOAB7QGAA9oDQAe0BoAPKA1AHhAawDwgNYA4AGtAcADWgOAB7QGAA9oDQAe0BoAPKA1AHhAawDwgNYA4AGtAcADWgOAB7QGAA9oDQAe0BoAPKA1AHhAawDwgNYA4AGtAcADWgOAB7QGAA9oDQAe0BoAPKA1AHhAawDwgNYA4AGtAcADWgOAB7QGAA9oDQAe0BoAPKA1AHhAawDwgNYA4AGtAcADWgOAB7QGAA//HyS8YUc37nTVAAAAAElFTkSuQmCC)
For the second question, let's load some data for a phenotypic trait - mammal home range size. (I've chosen
this trait because I know it has an intermediate value of \(\lambda\).)
mammal.data<-read.csv("mammalHR.csv",
row.names=1)
homeRange<-setNames(log(mammal.data$homeRange),
rownames(mammal.data))
homeRange
## U_maritimus U_arctos U_americanus N_narica
## 4.75013596 4.41642806 4.03953633 0.04879016
## P_locor M_mephitis M_meles C_lupus
## 0.13102826 0.91629073 -0.13926207 5.31222027
## C_latrans L_pictus C_aureus U_cinereoargenteus
## 3.80666249 5.07517382 2.20827441 0.09531018
## V_fulva H_hyaena C_crocuta A_jubatus
## 1.35325451 5.02912988 3.21887583 4.12874599
## P_pardus P_tigris P_leo T_bairdii
## 3.14415228 4.24276457 5.46383181 0.69314718
## C_simum D_bicornis E_hemionus E_caballus
## 1.89461685 2.74727091 3.55534806 3.11351531
## E_burchelli L_guanicoe C_dromedarius G_camelopardalis
## 5.10594547 -0.69314718 4.60517019 4.43793427
## S_caffer B_bison T_oryx G_granti
## 4.92725368 4.89034913 4.47163879 2.99573227
## G_thomsonii A_cervicapra M_kirki O_americanus
## 1.66770682 1.87180218 -3.14655516 3.12456515
## O_canadensis H_equinus A_melampus C_taurinus
## 2.66235524 4.38202663 1.33500107 4.31748811
## D_lunatus A_buselaphus A_americana C_canadensis
## 0.78845736 1.60943791 2.30258509 2.55955019
## D_dama A_alces R_tarandus O_virginianus
## 0.26236426 2.77881927 3.40119738 0.67294447
## O_hemionus
## 1.04731899
We can estimate \(\lambda\) as follows:
hrangeLambda<-phylosig(mammal.tree,
homeRange,method="lambda",
test=TRUE)
hrangeLambda
##
## Phylogenetic signal lambda : 0.415998
## logL(lambda) : -100.813
## LR(lambda=0) : 0.783684
## P-value (based on LR test) : 0.376017
It's also now possible to plot the likelihood surface for \(\lambda\) - which by default computes the surface
between 0 & 1 or the MLE of \(\lambda\) (whichever is higher).
plot(hrangeLambda)
![plot of chunk unnamed-chunk-5](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAA2sAAANrCAMAAADva/CiAAAA7VBMVEUAAAAAADoAAGYAOjoAOmYAOpAAZrY6AAA6AGY6OgA6Ojo6OmY6OpA6Zjo6ZmY6ZpA6ZrY6kJA6kLY6kNtmAABmADpmOgBmOjpmOpBmZjpmZmZmZpBmZrZmkJBmkLZmkNtmtttmtv+QOgCQOmaQZgCQZjqQZmaQkGaQkLaQtraQttuQtv+Q29uQ2/+2ZgC2Zjq2Zma2kDq2kGa2tma2tpC2tra2ttu229u22/+2/9u2///bkDrbkGbbtmbbtpDbtrbbttvb25Db27bb29vb2//b/9vb////tmb/tpD/25D/27b/29v//7b//9v///9tYYugAAAACXBIWXMAABM5AAATOQGPwlYBAAAgAElEQVR4nO3dC39Tx52H8REYtIjdhK5d0zgNBXeBJllMQ9OUjda4u2ZxfDl6/y9ndXSxdbMu58z85zczz/fzaQiU+Iz+Zx7raskNAFhwsRcAFILWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYMOgNXIGBrQGWKE1wAatATZoDbBBa4ANWgNs0Bpgg9YAG7QG2KA1wAatATZoDbBBa4ANWgNs0Bpgg9YAG7QG2KA1wAatATZoDbBBa4ANWgNs0Bpgg9ZUuDVirw0+0FpE9/a0ODLCywGtmdry6mrtyEguUbRmxHsgJJcaWgsvaBMklwxaC6tZB7v/JwSnj9aCabH9m/5n5KaM1kKIdy1DbrpozbfoN+aiLwCr0ZpXnrZ56y/C3TdBtOaPv+3t5+sQnBZa80J3W8surDy01pZuZrfkF1gGWmsnzC72/0WpLT5aayPUDg7xZaktNlprLLlbZsktODO01oj+fbTVEl12Hmhtd+E3bMivT26x0NpubK7QAh+C2qJIf+MYyucqIZsLkhJa21pe+zOf7xvJoLUtmW5Nm2NRmy1a24rxtrQ6GrVZorWNUn18fys5XzY1tLZeAXsx/0sogtbWibUNjQ9LbSZo7X7xtqD5gQu4+o6P1u5R3O4r7fLao7WVitx4RV5oQ7S2QvRNF+v40S941mhtkcKNx3grELjw2aK1eew1JhAKrc1in9WYQhi0dkdnj8VeiM4kckJrU0r7K/5SlKaRC1qbSGOVhqjNN1obYWOtwEz8orWBYGki65GbS9poTXFHyaxIcDbpajHLq7dHvV7v6MXHcIcwwG5aj/l403iSl4duau9dmEMYYCdtxog8aTrIy67b+/689vOh66yNTfdcyZamtS7ZMSWm6RhP3OOV/+7xEKEJbyG1lQmPKiENh3hz8ODX1b/xdojA2D47YVztldqa5KKkMbG2Gk6wOnavbn9z4R598X+IkNS/S0suT31o8prO78J1Xo4Dq067M915PEQw+ptGdIGiy0pF4/H160f7e71e/ct+mEOEoV+aLmbXRvPhXX94Mnp27eF3G57Mljo/7JZ2GF9zZb1GK5XShJeZyggFldRaOttEeqHSi1NWTmvplKaOQTbjY24pPL9GaR4xzEYKaU1hDTuQX678AhV5Gdr15+CHaCW578P6601upAJKuL8WfQFZorZd5d8aeyIUJrub7FtLcj+ksuhU1qkh99bS3A3JrJqrth3k3RpbITgmvLWmP7/2vDfrqeZj/uwDAwx5W01/fu21myX5/FrCV2pJrTzhOdtq8TM1a99kxMchWkl6ByS29sSWG0vzMZ2s/wlRH4doLunS0sO4t9F8SNXx2luOPg7RGKfeGhPfrMWMfjv6MfQhmkn/u2yCFyDBJVvL8DH/DM56ihch/e9woeXXGqc8Fia/XqD5zD0jEOYQ9x3Y8miYw+zXajEexc+pyeRsp3oxUl23jcbTUfycmmyu1JK9HMku3ELT4Sh+Tg0nOj7Owf2azkbvc2qyuVJLG2fhXk1feyz32RlZneOUL0zKaw8rk9Yyu1JL+tJkdi78afo6f63PqeHsSuF0rNR0LFKfU8O5FcMJWaXFz9SofE5NhrdZkr9EGZ6T9prPROVzanI8qxlcpgwugm+p35niG6gqTsyixFvjhOri3CxIu7VcT2cel4vbHPNSbi3fc5nLBcvlcviRcGucSH2coxnpvj8kpzEFnKU7yb4/ZNYnMaMLl9FFaSvR94fM967aSE6XLqfL0k6a7w/J+UsIJ2siyfeH5OwlhdM1luL7Q+Z/7jK7hJldnKbSe8w/87tqI7ldxNwuTzPJtcZpSxFnbdB6CNXZL+t+TNTHIQJ+MZjhvLWewYa3P/BxiLkvVcgpy/BiZniRdpVUa8WcrxwvaCnfJ++XUmvFn6zElX7+Emqt9FOVvsLPYDqtlXSicr2suV6u7STzOGRRpynbC5vtBdtGKs+vFX2SMlLyeUyktZJPUV4KPpNptFbaCcr48mZ80TZJorXizk/OFzjny7ZeAq3xJGheij2d+q0Ve2qyVeoZlW+tyBOT+YXO/OLdR721Mk9L7pc698u3mnhrZZ6U/BV5XrVbK/KUFKHEM6vcWrkPQBZwwQu4iIuEWyvwbEyVcNFLuIzzdFsr71wUprgTLNtacWeiPKWdYtXWSjsP8wq59IVczCnR1go7C4tKufilXM4xzdbKOgcFK+pES7ZW1BkoW0mnWrG1kua/WkET4KLGPERB079PSSMo57LqtVbO7DFSzAmXa62YyWOqlFOu1lopc1+vsCkUcnHFWitk6psUNoZCLq5Wa4UMHQvKOO9SrZUxciwr4swrtVbEwLdS3iRKuMRCrZUw7i0VOIoCLrJOawUMG2vkf/5lWst/1Fgv+x1Aa4rKHEbul1qltdznvJtCp5H5xRZpLfMpYzt5bwON1vKeMbaW9UagNUXlziPnSy7RWs4DbqTggWR80RVay3i82FXGm0GgtYyni93lux3it5bvbJsreibZXvjorWU72TaKHkq2F57WoCbXLRG7tVznihYy3RSRW8t0qm0VPpZML37c1upPM+R/K/4X/qxIy/Pyx74NCayQ5Z6hNSjKcdPQmiJGluMIaE0RI8txBLQGTfltG1qDqOz2Da0pYmS13KZAa4oY2UhmY8i6ter4wa/3/X+X3cezv724/2/WPh06t/dD/W/99X8RPuUVG61Nfvdq+Nffu73Vf7/vRvZHX/Pxyr+CELKKjdZGTh59GRdV/7ri77qXg8Gp67wb/ubCvfK+0EVZbbE2shoErY1/U/fz14+Dq9+/W/F3+6NrtMkv1fHKHr3Kaou1ktMkSmjt9IlzD1+O/uTq0HX+NMpsrrWT2yYv9qf/qZt49KU6Hl2hDf+TUWV9gys2TGUUW/6tTaupOxneGBz6arG1m4Pbf794PP1PZ1sbp3hzMPp14VEVhJVPbPm31nedHwaDT93hDcBhQM++DO93uYXWLtzk2mya07ybg/Gtxklz098GlM/+ai+fWWTf2uwNwNsbgQut3d0q7I8f/5i30Nr0KwaUz/7yIJthZN/abCiT66+l+2sn03iGV2t/21/6Mgut3f11mMgltoJa67zb2Frf7d/8cfKf3nt/jdasZRJbQa1tvl6rY6r+/Ov4r9/7OKRBa5lsLm/ymEf2re1yf230a3/5Ef2559e4v2Yvj3lk39qKxyG79zwOOb4KvNj7ODj909zXmXvdiMXjkFiQRWz5t7b8/NrkMf/pjcRX0+fXTkZXb6O/v3DdNvN6SJ5fiyKH2PJvbfF1I27v5+5Ca+PXjUwbujx0zxa/0t3r/IfhhX9oJIed5VkGI8m6tXvcPXU9ddldfqT/PrweMooMRlJUa5fdBx9H121LD3+cbN+Pxev8sUxnGzVVVGvTu27L97fGL/Tf7mtwby0OnX3UUFGtDar33du7bvM2/Fz2HZOfyxYamZDUp1JWa6lgZKukPhVaQzIS30m0hnSkvZVoTREjWy3tudCaIkZ2j6QHQ2tIScqbidaQlIR3E60pYmT3Sng0tKaIkd0v3dnQGhKT7H6iNSQm2f1Ea4oY2TqpTofWFDGydVKdDq0hOYnuKFpDetLcUrSmiJGtl+Z8aE0RI9sgyQHRGhKU5J6iNaQoxU1Fa4oY2UYJjojWFDGyjRIcEa0hTeltK1pDopLbV7SmiJFtIbkh0ZoiRraN1KZEa0hWYjuL1pCutLYWrSliZNtJa060poiRbSmpQdEaUpbS5qI1pCylzUVrihjZ1hIaFa0pYmRbS2hUtIa0pbO9aA1pS2d70ZoiRraDZIZFa4oY2S5SmRatIXWpbDBaQ/IS2WG0poiR7SSRcdGaIka2mzTmRWvIQBJ7jNaQgST2GK0pYmS7SmFitKaIke0sgZHRGrKQwC6jNeRBf5vRmiJGtjv9mdGaIkbWgPzQaA25UN9otIZcqG+0luu7etN7+uJL0EOUiJE1Ij62hsu7OXjw6/CXvqt1XoU4RMkYWSPiY2vV2oXrvDw//8l13gU4BLAr7a3WprXq2I2u0frucYBDALvS3mptWpvckBxcdse/+j1EyRhZQ9KD89Ha9Fe/hygZI2tIenCt7q+djO+oXXYfrXsoUvryIy/Km61xa849/Optt76jVp1wfw0ilDdb07WdvX1eP9w/uiXp1t6ElL74ohhZY8Kja7O06uztH+rWvl7/ZLbwpVfFyBoTHh2v0UJedLcbrSEvutuN1hQxshZkh+djYTy/5hsja0F2eLSG3KhuOC/ruv4c/BDAtlQ3HPfXFDGyVkTHR2uKGFkrouOjNeRHc8vRGvKjueVoTREja0lygE1f5/+8N+spj/l7xchakhxgw0VVr90snl+DFsVN13hNG95kxMchgKYUN13zNZ249W9V5+EQxWJkrQmOsPmSquP1PyLq4RDFYmStCY6wxZJ+O/ox9CGApvS2HY/5I096247WFDEyD+SGGGhBc88IhDlEzhiZB3JDbLGgq7dHvV7v6MXHcIcAmlPbeI3Xc3l4e721t/ajM+QuMgqhtvGaruey6/a+P6/9fMjn1PjGyLwQG2PT5cy+1zHve+wbI/NCbIztPutwxW+8HQJoS2vr0RrypbX1mr7O/3jm1ZAXjs+p8YuReSI1yKaLqT++dxxYddpd/ypkqcubBkbmidQgW/xMzVCv16t/2Q9zCKAtpc3XfC3XH56Mnl17+N2GJ7OVLi4Ko7T5eD2kIkbmjdAoaU0RI/NGaJS0hrzpbD9aQ950th+tKWJkHskMk/eHVMTIPJIZJu8PidypbEDeHzIF1fn5+bqXwWEdlQ3I+0MqmhvZ1evu6MZD5+tNPwGPlVQ2IO8PqWhmZJ/qn3/vHR0d1a+G2/tHvDUlTGQH8v6Q0qr37uGL6a3H6w9d94ybkrsT2YE85q/s5njvh7k/+HT4aLsbE5ilsQVpTdF0ZNUvS//XGVdsu9PYgi1XUZ39svHca1zQpDAyzyQG2nIRG97+wMchgNYk9uDCInZ9IofWkASFTTi7hgZP5NBaEKtGVr3vbPmEJpYpbMK7NTR6IofWglg5sm1mjXsobMLbR7yaPZFDa3Yuu2vfrgxrCezCyRKaPpHD45BWrk+7e7xGqzmBXThZQsgncgQuZWoWRva/H77pfbf5mxrWib8NeS5b0czIrj/0urfvxYnm4m9DWtPWd/ujd7vtvKC2duJvw1Ur8PyAV/wLma7pqTh9Qm1tRd+HtKbodmTV8fSz7T4dug4v8m8j+j6kNUV3I6ve/+s0sE+HPL/WRvR9SGsJ4X0QWom9EWlN2c3vXy78yem/cN3WVOyNSGuKbkd22p19jUF1+qSzGB+2F3kn0pqiu5FdvXbu6V9+OT8/P/tw6Bxv79OGRmvVm6M733RpTcfV89u34XxKaa1otHZzsP17qzY8BBo7ezv8Dvjd8svosKO4W3F6vXb291leX3tHaztjZGFItJb4IXLDyAKJOlhaU1a9+fbuFsbN87UfUYItKLQW8okcWmts/Ijw5HFhfi67PYXWQj6RQ2s7m46M1nyLuRlvjx3uiRxa2xmthSLRWrgncmitMVrzTaS1QaAncmitMVrzLuJu5HFIRdyGDIbWMIfWghFp7fp86nOoQ2AntOZfvO04c+S510R6fC9CWmvs5qDzVf1a8Pqfvl8SXiqJ1qo33zj39Oi5c51/e+Km73Ph9RDYzt1tyHAvCS+VRGvTN7Gujt3+oO8ehzgEtnL7zu8BXxJeKonWbt+zqW7O45vH0xqURNuPc/fXJrdR6n/xeEec1qBEo7Xb6zVai4vPhApIobXBiRvfXzsZ3le7cNyGjIfWQoq1IWePe3PgOl8dfdN19dWa2w9xCDRHa54otDZ6rf/Q11+G5/Vrb4950ZoftOaJRGtD52fnoQ+BjbgNGVSkHcnrIRXRWlASrV296fV633p+F0Ja84PWfFForT95LZC3h0WWD4FdhHyL3JLF2ZKzR71w7t8/D65/cu5VqENgK7weMqzorY1eB1nz+FrIhUNgO7weMqzorc29RivMIQAJUfYkraFAsVsb3oYc30/z+PqshUNgO4wssNitDe+nPagf7r/y+PqsxUNgK7f31z4czfqW+2ueRG9teMXmHvaeOL9Xa7TWwgmPQ4YRY1POHbP60B2e0c4Lv98+aa25vtv7/jzMWy6VLXprNe8vh6S13d2OrPpUv9/7P2KuJU8SrSV5iNzM3bL/8ITc/IuwK/lc0QRc17ftyc2reK3xednirl4Pc3vG45DeRLxeexPuwWVa29mqkV39B49D+mS/Lbm/pmhpZKOrNd8PEJeN1rBsFNrDFzzi7xWtYcEotL2/cI3mG62hdjuyi2FoT38gtBDM9yWtKbod2Ynr8HrIQGgNM+pXqPJ6yEBoDTP4ueyQrDcmrSliZBZoDYzMBq0BNmgNMGK8M2lNESMzQWtgZDZoDbBBa4AR261Ja4oYmQ1aAyOzQWuADVoDjJjuTVpTxMiM0FrxGJkRWgNs0BpgxHJz0poiRmaF1krHyKzQGmCD1gAjhruT1hQxMjO0VjhGZobWABu0Bhix2560poiR2aG1sjEyO7QG2KA1wIjZ/qQ1RYzMEK0VjZEZojXABq0BRqw2KK0pYmSWaK1kjMwSrQE2aA2wQWslY2SmjMZNa4oYmSlaA2zQGmCD1grGyEzRWsEYmS2bedMaQGuADVorFyOzRWvlYmS2aA0wYrJHaQ2gtXIxMmO0VixGZozWABu0Bhix2KS0poiRWaO1UjEya7QG2KA1wAatlYqRmcsjBDbOzhiZuTxCYONAXx4h0Br05RECre2MkZnLIwQ2zs4Ymb3wM6c1oEZrgA1aKxMjsyfd2tXbo16vd/TiY7hDlIqR2RNu7fLQTe29C3MIwI5ua5ddt/f9ee3nQ9dZGxutIQXB92nTA5y4xyv/3eMhCsbIIlBt7ebgwa+rf+PtECVjZBHQGmBDtbXq2L26/c2Fe/TF/yEAU6qtDfvqvBwHVp12Z7rzeIiCMbIIZFsb9OtH+3u9Xv3LfphDlIuRxRB66s2//vWHJ6Nn1x5+t+HJbDYOkqDbmtIhgPZorUSMLAZaKxEjiyGF1nh+DTmgNcBGCq0Nrj8HP0RZGFkUgcfO/TVFjCwKWgNs6LZ29vfJqyCrN9/yekikT7W10279E9mjl4zw2IhvjCwK0dYunHt6NPmJbFrzjZFFIdrayegFx6ej2GgNeQi7VRv/rOj4PUYu3DAzWkMeRFub5NV3j77Qmm+MLA7J1qrj6XtnDW9M0ppvjCwOydYm99cGo3dD+B2tIQuarV123SSwm0PnaA050GxtcHU4Dax6T2ueMbJIgg7eyxev/snrRrxiZJHot7biy84KcwjAO9XW+Jwa5EazNT6nJiBGFolka3xOTUiMLBLJ1vicGmRIsTU+OwNZCrlZaU0RI4tFsDU+pyYoRhaLYGt8Tg2ypNgan1ODHEm2xufUBMTIYtFsTekQuWFk0QQcPa0BM2gNsEFrhWFk0dBaYRhZNHqt3TzvzXrK60aQB73WqtdzPw3Ka7SQCb3W6uey176438chysXIolFsbXCy/pVZPg5RLEYWjWRr1fHaW44+DgGYC7ddW3zl345+DH0IwJpka0KHyA0ji4fWysLI4pFtrTr7Zd2Pifo4BGBJtrUNb3/g4xCAJVorCyOLh9bKwsjioTXABq0BNmRb43HIIBhZPLKtiRwiN4wsomDDpzVgDq0BNmitKIwsIlorCiOLiNYAG7QG2KC1ojCyiGitKIwsIloDbNAaYCTUhqU1RYwsJlorCSOLidYAG7QG2KC1kjCymGitJIwsJloDbNAaYIPWSsLIYqK1kjCymGgNMBJox9IasIDWCsLIoqK1gjCyqGgNsEFrgA1aKwgji4rWCsLIoqI1wAatATZorSCMLCpaKwgjiyvM/GkNWERrgA1aKwcji4vWysHI4qI1wAatATZorRyMLC5aKwcji4vWABu0BhgJsmdpTREji4zWisHIIqM1wAatATZorRiMLDJaKwYji4zWABu0BtigtWIwsshorRiMLDJaA2zQGmAkxKalNUWMLDZaKwUji43WABu0BtigtVIwsthorRSMLDZaA2zQGmCD1krByGKjtVIwsthoDbBBa4ANWisFI4uN1krByKILcApoDViB1gAbtFYIRhYdrRWCkUVHa4ANWgNs0FohGFl0tFYIRhYdrQE2aA2wQWuFYGTR0VohGFl0tAbYoDXAiP9tS2uKGFl8tFYGRhYfrQE2aA2wQWtlYGTx0VoZGFl8tAbYoDXABq2VgZHFR2tlYGTx0Rpgg9YAG7RWBkYWH62VgZEJ8H4SaA1YidYAG7RWBEYmgNaKwMgE0Bpgg9YAG7RWBEYmgNaKwMgE0Bpgg9YAG7RWBEYmgNaKwMgE0Bpgg9YAI743Lq0pYmQKaK0EjEwBrQE2aA2wQWslYGQKaK0EjEwBrQE2aA2wQWslYGQKaK0EjEwBrQE2xFq7etN7+uJL0EMAcYi0dnPw4NfhL31X67wKcYiSMTIFUq1duM7L8/OfXOddgEOUjJFJ8Hwa2rRWHbvRNVrfPQ5wCCAyodYmNyQHl93xr34PAUQm2Nr0V7+HKBkjkyDU2uBkfEftsvto3UORbJydMTIJKq059/Crt936jlp1wv015EijtcHZ2+f1w/2jW5Ju7U1IWkOiRFqrVWdv/1C39vX6J7NpbWeMTIJQazqHyA0jk0BrgA1aA2wItsbza74xMgm0VgBGJkGwtcH15+CHAMwpthb9EEAIfrcurSliZBpkWjv7++Q57OrNt7we0itGpkGktdOuc27vY/2vPDaCPGm0duHc06PD8U9k0xrypNHaidsf/vN0FBut+cbINEi0dnMw/tG1CzfzM6N+D1EyRqZBpLVJXn336AutAZs1DKE6nr531vDGJK0Bm7W7vzaoq3O/ozXPGFmOmp7Vy+70p7FvDt36H8xm4+yMkeWo8Vm9OpwGVr2nNWAjLyFU/+R1I8AGgUJws8IcImeMLEctzurV26Ner3f04mO4Q5SKkeWo8Vm9PLy93tpb+9EZbByg1uJxyL3vz2s/H/I5NcBmzZ9fe7zy3z0eomCMLEdtX6O19BtvhygZI8sRrQE2Gr8e0t19cO+F43NqgE2a/6xo5+U4sOq069Z+YDat7YyR5ajxWR19LH2v16t/2Q9ziHIxshw1P6vXH56Mnl17+N2GJ7PZOMCA96wDrNCaIkaWI1pTxMhyRGuADVoDbDR93cjz3qynvG7EK0aWo6avG3k999Og61+jBeTHrLX6uey1L+4PfHSvNFYhsgyNVYgsQ+K9WAf1T9KsfWVW4KP7pLEKkWVorEJkGSqtVcdrbzkGPrpPGqsQWYbGKkSWodLa4LejHyMe3SONVYgsQ2MVIsuQaS35o09prEJkGRqrEFkGrfmmsQqRZWisQmQZUq1VZ7+s+zHRwEf3RGMVIsvQWIXIMqRa2/D2B4GP7onGKkSWobEKkWXQmm8aqxBZhsYqRJZBa75prEJkGRqrEFkGrfmmsQqRZWisQmQZtOabxipElqGxCpFlSLXW8nFIoBwa3z+A/NEaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2wYd/ap0Pn9n5Y/ycRVnH9vuvcw5e2q7jnop903kVexafnznVeGP9o4vIyTp9EOCeDQX/+B6B97U/z1vrjT/nYX/cnEVZxdTD+Ew8fCNJmGZM/tW3t3lPyyDS25WWcRNgZQ6fzn7vkbX9at3bZdS/rS3O3m5b/JMYqTtyzz8M/6Xr4QJA2y6hdONvWVp2SzsvRt5+4w7hwnR/qc2L7naf6af4zzvztT+vW+uPvD/27bxPLfxJhFTcH4+/hXj7pqvkyapfdvUPT3bW8islHEF1EHkaMnTG4PHR7h7Ot+VuFcWvV8XgbXXant0+W/yTGKqYuu4bba+UyquMH/3Vs2dryKqbfeEytGMbtLre8ej3pvPi/2U9g8rg/zVsbX4y7999a/pMYq5gy/R66chknnXeVcWuLqxh+w6k/N9b2sZEVwxjelB3ehvzUtdsYQ+df5j/tzOP+NG5t+i3z7uIs/0mMVUxcmp7XVcuov4nbtra8igv3uwPzx0ZWDWN4e25oz/ZB2YVt4XF/0tqMS9uHRlYsY3QXKX5rbnyFYnklv+qc/He3bq1j/aA/rQVdxcjp6BEnO8vLGN8piN+a/V3oVeek7/Y+jpK3/AY4yKY14ftr1XvXsT2py8uYPJXj5h91tl7FtDHLb39rdoZp8rMHvmddjfE45O2fjb6LWlpeRpzWFlcRqbX7Hg41XcbSAdN9HFLkWZTlYw5Te2b+QPc9F932NuTyKoazGF3B216hrHjOM8JN2dp83Mk+vyb7uhHbJ0zvXcaIcWsrX7Dx0vyO0qrX8uz952DwP4fWp2a+tXRfNzL76rLJN67Ir4ccreLmYHrjzfR76PIwasatrVjF+8gvDh0vY3j9GuFlmbN30fzuz5iv859ur7iv8x+t4sJFaW3FMAb2ra1YRZQX2C8to/pQL8P6xzC0o0MAAAGOSURBVA0WW0v3df5AoWgNsEFrgA1aA2zQGmCD1gAbtAbYoDXABq0BNmgNsEFrgA1aA2zQGmCD1gAbtAbYoDXABq0BNmgNsEFrgA1aA2zQGmCD1gAbtAbYoDXABq0BNmgNsEFrgA1aA2zQGmCD1gAbtAbYoDXABq0BNmgNsEFreajeu71fN/81RERreehH+Ax37IbW8vDXj4Or39t+oj12RGv5uNiPvQKsQ2v5uHgcewVYh9aycXPwgEdHlNFaNvquwx02ZbSWi+HV2t+4w6aM1nLRd/s3f4y9CKxBa5mo761Vf+YOmzBay0TfvRr+41XsZeB+tJaHm4P6VSMXex8Hp3+KvRasRmt5OKmv1gbVsXOO6zZRtJaFy+74eezLQ/cs8lJwH1oDbNAaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2wQWuADVoDbNAaYIPWABu0BtigNcAGrQE2aA2w8f8SsJhjdMSpeAAAAABJRU5ErkJggg==)
There is no way to specify a starting value for the optimizer in phylosig
since (absent sampling error
in our trait) R is doing a univariate optimization on the interval (0,maxLambda
); however, we can easily
compute the likelihood of any particular value of \(\lambda\) as follows:
hrangeLambda$lik(0)
## [1] -101.2049
hrangeLambda$lik(0.5)
## [1] -100.8297
hrangeLambda$lik(0.7)
## [1] -101.0179
If we do have sampling error, then optimization is multidimensional (because we have to simultaneously
optimize \(\lambda\) and \(\sigma^2\)) & it is consequently imaginable that the starting value can affect
convergence.
Here is a small demo of how to try different initial values for optimization using optim( )
. (Note that
I could have a different numerical optimization function, such as nlm( )
, and it would work more or
less the same.)
## imagine some standard errors for each species:
se<-homeRange
se[]<-0.5
se
## U_maritimus U_arctos U_americanus N_narica
## 0.5 0.5 0.5 0.5
## P_locor M_mephitis M_meles C_lupus
## 0.5 0.5 0.5 0.5
## C_latrans L_pictus C_aureus U_cinereoargenteus
## 0.5 0.5 0.5 0.5
## V_fulva H_hyaena C_crocuta A_jubatus
## 0.5 0.5 0.5 0.5
## P_pardus P_tigris P_leo T_bairdii
## 0.5 0.5 0.5 0.5
## C_simum D_bicornis E_hemionus E_caballus
## 0.5 0.5 0.5 0.5
## E_burchelli L_guanicoe C_dromedarius G_camelopardalis
## 0.5 0.5 0.5 0.5
## S_caffer B_bison T_oryx G_granti
## 0.5 0.5 0.5 0.5
## G_thomsonii A_cervicapra M_kirki O_americanus
## 0.5 0.5 0.5 0.5
## O_canadensis H_equinus A_melampus C_taurinus
## 0.5 0.5 0.5 0.5
## D_lunatus A_buselaphus A_americana C_canadensis
## 0.5 0.5 0.5 0.5
## D_dama A_alces R_tarandus O_virginianus
## 0.5 0.5 0.5 0.5
## O_hemionus
## 0.5
hrangeLambda.se<-phylosig(mammal.tree,
homeRange,se=se,method="lambda")
hrangeLambda.se
##
## Phylogenetic signal lambda : 0.441654
## logL(lambda) : -100.813
## MLE(sig2) : 0.0595599
## create temporary likelihood function
LIK<-function(par) hrangeLambda.se$lik(par[1],par[2])
## optimize using starting values 0.5 & 0.1
fit1<-optim(c(0.5,0.1),LIK,method="L-BFGS-B",
lower=c(0,0),upper=c(1,Inf),
control=list(fnscale=-1))
fit1
## $par
## [1] 0.44160308 0.05955602
##
## $value
## [1] -100.8131
##
## $counts
## function gradient
## 49 49
##
## $convergence
## [1] 0
##
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
## optimize using starting values of 1.0 & 10
fit2<-optim(c(1,10),LIK,method="L-BFGS-B",
lower=c(0,0),upper=c(1,Inf),
control=list(fnscale=-1))
fit2
## $par
## [1] 0.4416549 0.0595599
##
## $value
## [1] -100.8131
##
## $counts
## function gradient
## 33 33
##
## $convergence
## [1] 0
##
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
## optimize using starting values of 0.0 & 100
fit3<-optim(c(0,100),LIK,method="L-BFGS-B",
lower=c(0,0),upper=c(1,Inf),
control=list(fnscale=-1))
fit3
## $par
## [1] 0.44163973 0.05955861
##
## $value
## [1] -100.8131
##
## $counts
## function gradient
## 43 43
##
## $convergence
## [1] 0
##
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
In all these cases, we converge on the same solution - but for a more problematic optimization this might
not be the case.