Hi all,
I have a specific & more general question about estimating branch length
transformations.
1. Specific question:
I ran a gls using corPagel(..., fixed=FALSE) to estimate Pagel's lambda while
fitting the model. The model gives me a negative value of lambda which seems
very wrong. Does it suggest a problem with the data, simply no phylogenetic
signal, or something else? I am guessing it means no phylogenetic signal
because setting corPagel(value=0, ..., fixed=TRUE) gives similar but not
identical model. The code below & attached dataset should reproduce this.
I hope it's not too naive; I've searched, but cannot find anything about
negative lambda estimates.
2. General question:
When setting corPagel(value=x, ..., fixed=FALSE), the estimated lambda varies
depending on the value (x) specified. This also happens with other "corPhyl"
methods. Is there a preferred way of putting a sensible starting value? I'm
thinking an iterative process may do the trick, but again I'm guessing...
Thanks in advance for your input,
Mick
(PhD candidate)
McGill University
Montreal, Canada
### Here's a reproducible example with data in the attached file (dd.test.txt):
### read data file
dd <- read.csv("testdd.txt", header=T, row.names=1)
### read phylogeny
ttxt <- "((((((((Sit_avenae:1,Met_dirhodum:1):1,Mac_creelii:1):1,Acy_pisum:
1):1,Uro_jaceae:1):1,(Myz_persicae:1,(Bre_brassicae:1,Diu_noxia:1):1):1):
1,((Rho_maidis:1,Rho_padi:1):1,(Tox_citricida:1,((Aph_craccivora:1,Aph_fabae:
1):1,(Aph_glycines:1,Aph_gossypii:1):1):1):1):1):1,Dre_platanoides:1):
1,(Cin_piceicola:1,Cin_pilicornis:1):1);"
tree <- read.tree(text=ttxt)
### fit gls
corPa <- corPagel(value=1, phy=tree, fixed=FALSE)
form <- formula(th ~ aSize + pSize + instar)
summary( gls(form, correlation=corPa, data=dd) )
### output:
Generalized least squares fit by REML
Model: form
Data: dd
AIC BIC logLik
34.20101 38.03535 -11.10050
Correlation Structure: corPagel
Formula: ~1
Parameter estimate(s):
lambda
-0.6050088
Coefficients:
Value Std.Error t-value p-value
(Intercept) 1.8141584 0.6354861 2.8547571 0.0127
aSize 1.3394331 0.5901446 2.2696693 0.0396
pSize -1.0607982 0.5866774 -1.8081458 0.0921
instar 0.3872971 0.4276060 0.9057336 0.3804
Correlation:
(Intr) aSize pSize
aSize -0.613
pSize 0.058 -0.702
instar -0.745 0.383 -0.316
Standardized residuals:
Min Q1 Med Q3 Max
-1.7951511 -0.4593993 0.4495497 0.9947766 2.7149596
Residual standard error: 0.3977355
Degrees of freedom: 18 total; 14 residual
### session info
R version 2.9.0 (2009-04-17)
i386-pc-mingw32
locale:
LC_COLLATE=French_Canada.1252;LC_CTYPE=French_Canada.1252;LC_MONETARY=French_Canada.1252;LC_NUMERIC=C;LC_TIME=French_Canada.1252
attached base packages:
[1] grDevices utils datasets stats graphics methods
[7] base
other attached packages:
[1] geiger_1.2-14 ouch_2.5-4 subplex_1.1-3 msm_0.8.2
[5] mvtnorm_0.9-7 MASS_7.2-46 picante_0.7-0 nlme_3.1-92
[9] vegan_1.15-2 ape_2.3 geepack_1.0-16
loaded via a namespace (and not attached):
[1] gee_4.13-13 grid_2.9.0 lattice_0.17-22 splines_2.9.0
[5] survival_2.35-4
"","th","aSize","pSize","instar","rSize"
"Sit_avenae",2.34057803811701,0.941119058572024,0.850597568334247,0.943826746719236,1.03702283879018
"Met_dirhodum",2.34563724252011,0.973256826044084,0.871593656376963,1.03972077083992,0.837774244589395
"Mac_creelii",3.15102515796003,1.46851240142544,1.05239289777758,0.693147180559945,0.745315882
"Acy_pisum",2.80330188331362,1.20922931205038,1.04429028121510,0.946562873037042,0.878868303935916
"Uro_jaceae",3.69494564004702,1.06248029036020,0.936093359170335,0.95749834867258,0.850708924
"Myz_persicae",2.85650839099218,0.726036192173243,0.659734564359011,0.907262024479215,1.28951922680894
"Bre_brassicae",2.06595477383445,0.841323509008712,0.792615771267286,0.818868912327109,0.973140323
"Diu_noxia",2.03022150527321,0.655123972551085,0.664715363747263,0.80471895621705,1.01630487206375
"Rho_maidis",2.99573227355399,0.837425781805104,0.52669897272543,1.09861228866811,0.801577909
"Rho_padi",2.31253542384721,0.640936941696746,0.714812226888396,1.38629436111989,1.136077196
"Tox_citricida",2.53877299956916,0.65836106466922,0.644370010706087,0.549306144582964,0.93082912
"Aph_craccivora",1.85629799036563,0.648105167189013,0.693147180559945,1.6094379124341,1.078894134
"Aph_fabae",3.61361696961339,0.775801932182208,0.545227050483323,1.09861228866811,0.78125
"Aph_glycines",3.55915232109702,1.21191391735909,0.573701226418771,1.03972077092965,1.40784782247920
"Aph_gossypii",2.87495120513488,0.318754489522555,0.703979703724171,1.49786613677700,1.28995075405653
"Dre_platanoides",3.10863384891252,0.720903276377944,0.699092446888216,1.6094379124341,0.672285347
"Cin_piceicola",2.44560557187222,1.09861228866811,1.34807314829969,1.09861228866811,1.222222222
"Cin_pilicornis",2.7797505152674,1.46511117900518,1.34807314829969,1.09861228866811,0.977281381
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
