Re: [R-sig-phylo] ape-package/ possibility to get the transformed residuals (phylogentic corrected residuals)?

2010-10-12 Thread Simon Blomberg
Actually, the function residuals.gls calculates 3 different types of 
residuals: response or raw residuals, pearson or standardized 
residuals, and normalized residuals. For diagnostic purposes, you 
probably want the normalized residuals. The default is the response  
residuals (although the ?residuals.gls says that the pearson residuals 
are the default. This seems to be a bug in the help file). See 
?residuals.gls


Also, it's better programming practice to use the extractor function 
(residuals.gls in this case), rather than using result$residuals 
directly, since the contents of slots can change on the whim of the 
developer, but the extractor function should, in general, return what 
you ask for.


Cheers,

Simon.

On 13/10/10 03:40, Liam J. Revell wrote:

Hi Jan:

Yes.  gls() returns an object of class gls which contains a vector 
containing the residuals as one of its components.

So for instance you might compute:
 
result-gls(y~x,data=test.data,correlation=corPagel(phy=tree,value=1,fixed=FALSE)) 


in which case:
 result$residuals
contains the model residuals.
You can also get the residuals using the function resid(), i.e.:
 resid(result)

Residual standard error is:
 result$sigma
so residual sum of squares should just be:
 resid.SS-n*result$sigma^2
for n species.

Hope this is helpful. - Liam

Liam J. Revell
NESCent, Duke University
web: http://anolis.oeb.harvard.edu/~liam/
NEW email: lrev...@nescent.org



Werner, Jan wrote:

Dear R-sig-phylo members,
for some analyses I used the gls and gnls functions (package nlme) 
and controlled for phylogenetic effects by the phylogenetic 
correlation structure corPagel provided by the ape-package. My 
question is, is there a possibility to get the transformed residuals 
(phylogentic corrected

residuals) or the phylogenetic corrected residual sum of squares?

Many thanks for your efforts.

Best regards

Jan Werner

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo


___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo


--
Simon Blomberg, BSc (Hons), PhD, MAppStat, AStat
Lecturer and Consultant Statistician
School of Biological Sciences
The University of Queensland
St. Lucia Queensland 4072
Australia
T: +61 7 3365 2506
email: S.Blomberg1_at_uq.edu.au
http://www.uq.edu.au/~uqsblomb/

Policies:
1.  I will NOT analyse your data for you.
2.  Your deadline is your problem

Statistics is the grammar of science - Karl Pearson.

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo


[R-sig-phylo] geiger: applying fitDiscrete to multiple variables over a tree sample

2010-10-12 Thread Lara Poplarski
Dear List,

I am new to R, and hope someone can kindly help with the following
task. I have a Bayesian sample of trees in nexus format and discrete
data; example trees and data are at the bottom of the email. I would
like to use fitDiscrete in geiger to estimate parameter lambda for all
variables. The idea is to then compare the distributions across trees
of lambda values for the different variables.

I have been using the following:

for(i in 1:5) {
cat(t,i,\n);
fd-fitDiscrete(trees[[i]],data,treeTransform=lambda)
print(fd)
}

A few specific questions:

1. Since trees in my sample are non-ultrametric, I get the following warning:
Warning: some tree transformations in GEIGER might not be sensible for
nonultrametric trees.

Is the lambda transformation sensible for non-ultrametric trees? I
could not find further information on this in the manual or in the
Pagel references.

2. In some cases, I get the following warning:
[1] Warning: may not have converged to a proper solution.

Does it make sense to get fitDiscrete to repeat the ML estimation for
these cases, until convergence? Can someone suggest how to modify the
loop above to do this?

3. Finally, I would be grateful for pointers on how to tackle the
output. For example, can someone suggest how to go about calculating
the mean lambda value, across trees, for each of the 4 variables?

I also have a more general question. Both the tree sample and dataset
are relative large (1000 trees * 80 variables), so based on
preliminary runs I expect the analysis to take rather long. Any
suggestions on how to speed things up and/or on a different approach
to the tree sample/multiple variable set-up would be most welcome!

Many thanks in advance,
Lara Poplarski

=
trees
=
#NEXUS
[R-package APE, Tue Oct 12 15:40:42 2010]

BEGIN TAXA;
DIMENSIONS NTAX = 5;
TAXLABELS
taxon4
taxon5
taxon1
taxon2
taxon3
;
END;
BEGIN TREES;
TRANSLATE
1   taxon4,
2   taxon5,
3   taxon1,
4   taxon2,
5   taxon3
;
TREE * UNTITLED = [R]
(((1:0.03870620631,2:0.03870620631):0.01327593656,(3:0.009785519975,4:0.009785519975):0.04219662289):0.9758290676,5:1.02781121);
TREE * UNTITLED = [R]
(5:1.546301171,((1:0.1587038879,(2:0.1024085511,3:0.1024085511):0.05629533677):1.183393702,4:1.34209759):0.204203581);
TREE * UNTITLED = [R]
(2:0.09245329862,(3:0.07043250347,(1:0.7179110514,(4:0.8067701138,5:0.03902596165):0.9586401181):0.6771807231):0.4711165503);
TREE * UNTITLED = [R]
((4:0.3565618952,3:0.6832435813):0.6823437791,((5:0.6852910472,1:0.5428841521):0.9435752912,2:0.7694561274):0.5071240654);
TREE * UNTITLED = [R]
((4:0.2373158785,3:0.3414158912):0.9065216579,(2:0.995543058,(5:0.8844613975,1:0.7886170356):0.9908593148):0.627468311);
END;

=
data
=

V01 V02 V03 V04
taxon1  h   h   a   NA
taxon2  g   h   d   NA
taxon3  h   j   h   h
taxon4  g   h   g   NA
taxon5  i   j   NAa

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo