Re: [R-sig-phylo] fitting cuadratic model using gnls
Thanks Ted and Ben for your useful comments. I used poly() to generate orthogonal polynomials of degree 2 and it seems to work well. Cheers, Jorge -Mensaje original- De: Ben Bolker [mailto:bbol...@gmail.com] Enviado el: martes, 10 de diciembre de 2013 19:44 Para: r-sig-phylo@r-project.org; Jorge Sánchez Gutiérrez Asunto: Re: [R-sig-phylo] fitting cuadratic model using gnls -BEGIN PGP SIGNED MESSAGE- Hash: SHA1 On 13-12-10 12:57 PM, Theodore Garland Jr wrote: > I don't know, but one thing to watch out for is a high correlation > between X and Xsquared. That gives multicollinearity and can lead to > numerical problems. One thing that is commonly done is: 1. > standardize the X variable (subtract mean and divide by standard > deviation) 2. square that value and use it for the quadratic This is > referred to as an "orthogonol polynomial." > > Cheers, Ted > Yes. Therefore, you would be (much) better off wit m3 <- gls(duration ~ poly(latitude_used,2), data=data, correlation = dataCor) (A quadratic model is still a *linear* fit in the statistical sense.) poly() will generate orthogonal polynomials by default, which will be more stable and are more appropriate if you're planning to test the significance of terms. If it is important to you to have the parameters interpretable as slope and curvature at the equator (i.e. where latitude=0), you may be able to get away with m3 <- gls(duration ~ poly(latitude_used,2,raw=TRUE), data=data, correlation = dataCor) > Theodore Garland, Jr., Professor Department of Biology University of > California, Riverside Riverside, CA 92521 Office Phone: (951) > 827-3524 Wet Lab Phone: (951) 827-5724 Dry Lab Phone: (951) > 827-4026 Home Phone: (951) 328-0820 Skype: theodoregarland > Facsimile: (951) 827-4286 = Dept. office (not confidential) > Email: tgarl...@ucr.edu > http://www.biology.ucr.edu/people/faculty/Garland.html > http://scholar.google.com/citations?hl=en&user=iSSbrhwJ > > Inquiry-based Middle School Lesson Plan: "Born to Run: Artificial > Selection Lab" > http://www.indiana.edu/~ensiweb/lessons/BornToRun.html > > Fail Lab: Episode One > http://testtube.com/faillab/zoochosis-episode-one-evolution > http://www.youtube.com/watch?v=c0msBWyTzU0 > > From: > r-sig-phylo-boun...@r-project.org > [r-sig-phylo-boun...@r-project.org] on behalf of Jorge Sánchez > Gutiérrez [jorgesgutier...@unex.es] Sent: Tuesday, December 10, > 2013 9:48 AM To: r-sig-phylo@r-project.org Subject: [R-sig-phylo] > fitting cuadratic model using gnls > > Hello everyone, > > > > Im trying to fit nonlinear models using the function gnls in ape. > For example, I run a simple quadratic regression of moult duration > (duration) and latitude (latitude_used) as follows: > > > > mod <- duration ~ latitude_used +I(latitude_used^2) > > init<-list(alpha=1,beta=1) > > m3 <- gnls(mod, data=data, start=init, correlation = dataCor) > > > > However, I get the following error message: > > > > Error in gnls(mod, data=data, start = init, correlation = > dataCor) : > > step halving factor reduced below minimum in NLS step > > > > Ive also tried with several correlation structures and different > starting values for estimation, but without success. I would be very > grateful if you can provide any suggestions. > > > > Thanks in advance, > > Jorge S. Gutiérrez > > > [[alternative HTML version deleted]] > > > ___ R-sig-phylo mailing > list - R-sig-phylo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive > at http://www.mail-archive.com/r-sig-phylo@r-project.org/ > -BEGIN PGP SIGNATURE- Version: GnuPG v1.4.11 (GNU/Linux) Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/ iQEcBAEBAgAGBQJSp2DoAAoJEOCV5YRblxUHn4gH/jYtxFo6aj2DlHAPw2FXufA4 CwtJmIHpWlnD5mhg5YqBgJNldn2XPxqx/qEf0mMzRgd2gAyIQiXd91WI6sqkaVhF FmYSPzClqx14jZSlF64KE54l/8W4mp81rDyYnh+/SSebltCrX8Cv7NqYqvNDhYQC M7i/avh0NNv/PZuDuUwbj95doQmdeUyeSlig0sxHaP14C7FxzvxODUBhcTB6yQzN 0wuK+iIaUKqSGm3jmcpY08oy0RJ34QfyH175UezRmSaE/RI+XCohhBbdtU5+WgRX lBRYmwtCIGAJ+DNmhlnhE3u6ZFdnzhw95gYd1LGL3Sr+JZl1U+VIKfIjhlrTgLQ= =lLum -END PGP SIGNATURE- ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] Mismatch distribution (pegas): expected versus empirical and other questions
Hi Alistair, I wrote a function to make the "standard" expectation under constant population size quite a while ago. It's up here https://github.com/dwinter/Pegas/blob/new_functions/R/MMD.stable.R It's been a little while since I've had to think about mismatch distributions, so, though I'm sure I would checked everything out at the time, you should probably confirm this function behaves as expected. (BTW, Emmanuel, if you want to include any of that function in PEGAS you are most welcome to) David Winter On Tue, Dec 10, 2013 at 6:27 AM, Alastair Potts wrote: > Hi all, > > The mismatch distribution function in pegas (MMD) produces the familiar > histogram and then also a line that I always took as the theoretical > "expected" distribution under a scenario of constant population size. > > I had a quick look at the MMD function and noticed that this line is not > calculated using the theoretical expectations based on theta, but rather is > some sort of moving-average based on the density() function. > > Does anyone know how to implement the correct theoretical expectation? I > have dipped into the literature and am now completely lost. > > Also, there are a range of statistics (such as the raggedness index) and > Arlequin has somehow managed to calculate a significance value for these > statistics. Has anyone implemented (or can easily implement) these > statistics and significance tests in R? > > Thanks in advance for any comments, thoughts and solutions. > > Cheers, > Alastair Potts > > [[alternative HTML version deleted]] > > ___ > R-sig-phylo mailing list - R-sig-phylo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-phylo > Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/ ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] fitting cuadratic model using gnls
-BEGIN PGP SIGNED MESSAGE- Hash: SHA1 On 13-12-10 12:57 PM, Theodore Garland Jr wrote: > I don't know, but one thing to watch out for is a high correlation > between X and Xsquared. That gives multicollinearity and can lead > to numerical problems. One thing that is commonly done is: 1. > standardize the X variable (subtract mean and divide by standard > deviation) 2. square that value and use it for the quadratic This > is referred to as an "orthogonol polynomial." > > Cheers, Ted > Yes. Therefore, you would be (much) better off wit m3 <- gls(duration ~ poly(latitude_used,2), data=data, correlation = dataCor) (A quadratic model is still a *linear* fit in the statistical sense.) poly() will generate orthogonal polynomials by default, which will be more stable and are more appropriate if you're planning to test the significance of terms. If it is important to you to have the parameters interpretable as slope and curvature at the equator (i.e. where latitude=0), you may be able to get away with m3 <- gls(duration ~ poly(latitude_used,2,raw=TRUE), data=data, correlation = dataCor) > Theodore Garland, Jr., Professor Department of Biology University > of California, Riverside Riverside, CA 92521 Office Phone: (951) > 827-3524 Wet Lab Phone: (951) 827-5724 Dry Lab Phone: (951) > 827-4026 Home Phone: (951) 328-0820 Skype: theodoregarland > Facsimile: (951) 827-4286 = Dept. office (not confidential) > Email: tgarl...@ucr.edu > http://www.biology.ucr.edu/people/faculty/Garland.html > http://scholar.google.com/citations?hl=en&user=iSSbrhwJ > > Inquiry-based Middle School Lesson Plan: "Born to Run: Artificial > Selection Lab" > http://www.indiana.edu/~ensiweb/lessons/BornToRun.html > > Fail Lab: Episode One > http://testtube.com/faillab/zoochosis-episode-one-evolution > http://www.youtube.com/watch?v=c0msBWyTzU0 > > From: > r-sig-phylo-boun...@r-project.org > [r-sig-phylo-boun...@r-project.org] on behalf of Jorge Sánchez > Gutiérrez [jorgesgutier...@unex.es] Sent: Tuesday, December 10, > 2013 9:48 AM To: r-sig-phylo@r-project.org Subject: [R-sig-phylo] > fitting cuadratic model using gnls > > Hello everyone, > > > > I’m trying to fit nonlinear models using the function “gnls” in > ape. For example, I run a simple quadratic regression of moult > duration (duration) and latitude (latitude_used) as follows: > > > > mod <- duration ~ latitude_used +I(latitude_used^2) > > init<-list(alpha=1,beta=1) > > m3 <- gnls(mod, data=data, start=init, correlation = dataCor) > > > > However, I get the following error message: > > > > “Error in gnls(mod, data=data, start = init, correlation = > dataCor) : > > step halving factor reduced below minimum in NLS step” > > > > I’ve also tried with several correlation structures and different > starting values for estimation, but without success. I would be > very grateful if you can provide any suggestions. > > > > Thanks in advance, > > Jorge S. Gutiérrez > > > [[alternative HTML version deleted]] > > > ___ R-sig-phylo > mailing list - R-sig-phylo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable > archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/ > -BEGIN PGP SIGNATURE- Version: GnuPG v1.4.11 (GNU/Linux) Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/ iQEcBAEBAgAGBQJSp2DoAAoJEOCV5YRblxUHn4gH/jYtxFo6aj2DlHAPw2FXufA4 CwtJmIHpWlnD5mhg5YqBgJNldn2XPxqx/qEf0mMzRgd2gAyIQiXd91WI6sqkaVhF FmYSPzClqx14jZSlF64KE54l/8W4mp81rDyYnh+/SSebltCrX8Cv7NqYqvNDhYQC M7i/avh0NNv/PZuDuUwbj95doQmdeUyeSlig0sxHaP14C7FxzvxODUBhcTB6yQzN 0wuK+iIaUKqSGm3jmcpY08oy0RJ34QfyH175UezRmSaE/RI+XCohhBbdtU5+WgRX lBRYmwtCIGAJ+DNmhlnhE3u6ZFdnzhw95gYd1LGL3Sr+JZl1U+VIKfIjhlrTgLQ= =lLum -END PGP SIGNATURE- ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] fitting cuadratic model using gnls
I don't know, but one thing to watch out for is a high correlation between X and Xsquared. That gives multicollinearity and can lead to numerical problems. One thing that is commonly done is: 1. standardize the X variable (subtract mean and divide by standard deviation) 2. square that value and use it for the quadratic This is referred to as an "orthogonol polynomial." Cheers, Ted Theodore Garland, Jr., Professor Department of Biology University of California, Riverside Riverside, CA 92521 Office Phone: (951) 827-3524 Wet Lab Phone: (951) 827-5724 Dry Lab Phone: (951) 827-4026 Home Phone: (951) 328-0820 Skype: theodoregarland Facsimile: (951) 827-4286 = Dept. office (not confidential) Email: tgarl...@ucr.edu http://www.biology.ucr.edu/people/faculty/Garland.html http://scholar.google.com/citations?hl=en&user=iSSbrhwJ Inquiry-based Middle School Lesson Plan: "Born to Run: Artificial Selection Lab" http://www.indiana.edu/~ensiweb/lessons/BornToRun.html Fail Lab: Episode One http://testtube.com/faillab/zoochosis-episode-one-evolution http://www.youtube.com/watch?v=c0msBWyTzU0 From: r-sig-phylo-boun...@r-project.org [r-sig-phylo-boun...@r-project.org] on behalf of Jorge Sánchez Gutiérrez [jorgesgutier...@unex.es] Sent: Tuesday, December 10, 2013 9:48 AM To: r-sig-phylo@r-project.org Subject: [R-sig-phylo] fitting cuadratic model using gnls Hello everyone, I’m trying to fit nonlinear models using the function “gnls” in ape. For example, I run a simple quadratic regression of moult duration (duration) and latitude (latitude_used) as follows: mod <- duration ~ latitude_used +I(latitude_used^2) init<-list(alpha=1,beta=1) m3 <- gnls(mod, data=data, start=init, correlation = dataCor) However, I get the following error message: “Error in gnls(mod, data=data, start = init, correlation = dataCor) : step halving factor reduced below minimum in NLS step” I’ve also tried with several correlation structures and different starting values for estimation, but without success. I would be very grateful if you can provide any suggestions. Thanks in advance, Jorge S. Gutiérrez [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] fitting cuadratic model using gnls
Hello everyone, Im trying to fit nonlinear models using the function gnls in ape. For example, I run a simple quadratic regression of moult duration (duration) and latitude (latitude_used) as follows: mod <- duration ~ latitude_used +I(latitude_used^2) init<-list(alpha=1,beta=1) m3 <- gnls(mod, data=data, start=init, correlation = dataCor) However, I get the following error message: Error in gnls(mod, data=data, start = init, correlation = dataCor) : step halving factor reduced below minimum in NLS step Ive also tried with several correlation structures and different starting values for estimation, but without success. I would be very grateful if you can provide any suggestions. Thanks in advance, Jorge S. Gutiérrez [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] Mismatch distribution (pegas): expected versus empirical and other questions
Hi all, The mismatch distribution function in pegas (MMD) produces the familiar histogram and then also a line that I always took as the theoretical "expected" distribution under a scenario of constant population size. I had a quick look at the MMD function and noticed that this line is not calculated using the theoretical expectations based on theta, but rather is some sort of moving-average based on the density() function. Does anyone know how to implement the correct theoretical expectation? I have dipped into the literature and am now completely lost. Also, there are a range of statistics (such as the raggedness index) and Arlequin has somehow managed to calculate a significance value for these statistics. Has anyone implemented (or can easily implement) these statistics and significance tests in R? Thanks in advance for any comments, thoughts and solutions. Cheers, Alastair Potts [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
