[R] finding roots of multivariate equation
Hello, I want to find the roots of an equation in two variables. I am aware of the uniroot function, which can do this for a function with a single variable (as I understand it...) but cannot find a function that does this for an equation with more than one variable. I am looking for something implementing similar to a Newton-Raphson algorithm. Thanks. -- Bill Shipley North American Editor for Annals of Botany Subject Editor for Ecology Département de biologie Université de Sherbrooke Sherbrooke (Québec) J1K 2R9 Canada __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] installing new packages
Hello, I have just installed the newest version of R (2.4.1) for Windows XP. I can no longer install new packages. When trying to connect to a server (I have tried several) I get the following message: > chooseCRANmirror() Error in open.connection(file, "r") : unable to open connection In addition: Warning message: unable to connect to 'cran.r-project.org' on port 80. Have other people had the same problem with this version, or is it unique to my computer? Can someone suggest a solution? Thanks. Bill Shipley [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Duncan post-hoc tests in R?
Hello, I am looking for an R function that impliments Duncan's post-hoc test. I am aware of multcomp and its function "glht" but, unless I am missing something, this cannot impliment the Duncan test. Any help of pointers are welcome. Bill Shipley __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] dynamic loading error with Open Watcom object file
Hello. I am trying to use a FORTRAN subroutine from within R (Windows
version). This fortran subroutine is compiled using the Open Watcom Fortran
compiler and the compiled object file is called ritscale.obj. Following the
explanation on pages 193-194 of "The New S language" I use the dyn.load
command:
> dyn.load("f:/maxent/ritscale.obj")
Error in dyn.load(x, as.logical(local), as.logical(now)) :
unable to load shared library 'f:/maxent/ritscale.obj':
LoadLibrary failure: %1 n'est pas une application Win32 valide.
The error message says: LoadLibrary failure: %1 is not a valid Win32
application
I do not know what this means. Can someone help?
Bill Shipley
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] setting new working directories
Hello, and Happy New Year. My default working directory is getting very cluttered. I know that I should be using a different working directory for each project (I work in Windows), but do not know how to go about creating different ones and moving back and forth between them. I have read Venables & Ripley (Modern Applied Statistics with S-PLUS, 1994) but this seems out of date with respect to this topic and have searched through the documentation but cannot find a clear explanation for doing this. Can someone point me to the proper documentation for creating and using different working directories from within Windows (please, no comments about switching to UNIX...). Thanks. Bill Shipley __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] help with syntax of nlme call.
I am getting an error message in a call to nlme and cannot understand what is happening. I explain the steps below in the hope that someone can explain the error and how to correct it. STEP 1: Data set: name: marouane.data. This is a data frame whose first few lines are as follows: > marouane.data[1:13,] species plant leaf irradiance photosynthesis chlorophyll 1 Asclepias.incarnata 11 0 -2.091359 0.02619 2 Asclepias.incarnata 11 50 1.153241 0.02619 3 Asclepias.incarnata 11100 2.241963 0.02619 4 Asclepias.incarnata 11200 3.541258 0.02619 5 Asclepias.incarnata 11300 5.012547 0.02619 6 Asclepias.incarnata 11400 5.710689 0.02619 7 Asclepias.incarnata 11600 6.632571 0.02619 8 Asclepias.incarnata 11800 7.314284 0.02619 9 Asclepias.incarnata 11 1000 8.092402 0.02619 10 Asclepias.incarnata 11 1310 8.110368 0.02619 11 Asclepias.incarnata 12 0 -2.051934 0.02707 12 Asclepias.incarnata 12 50 1.213854 0.02707 13 Asclepias.incarnata 12100 2.271453 0.02707 absorbance nitrogen sla 1 0.8816 18.35913 77.53 2 0.8816 18.35913 77.53 3 0.8816 18.35913 77.53 4 0.8816 18.35913 77.53 5 0.8816 18.35913 77.53 6 0.8816 18.35913 77.53 7 0.8816 18.35913 77.53 8 0.8816 18.35913 77.53 9 0.8816 18.35913 77.53 10 0.8816 18.35913 77.53 11 0.8813 18.35913 77.53 12 0.8813 18.35913 77.53 13 0.8813 18.35913 77.53 The nesting structure is species (25), plants within species (4 each species) and leaves within plants (2 each plant). There is only 1 missing value in the data set. STEP 2: I constructed a self starting function (photo) that gives the nonlinear function and provides initial estimates of the three parameters. This self starting function works. I use this to call the nlsList function, which fits separate nonlinear functions to each species (I am only working on a 2 level model so far: between and within species): fit<-nlsList(photosynthesis~photo(irradiance,Am,Q,LCP)|species, + data=marouane.data,na.action=na.omit) This works and I get the three parameter estimates per species. STEP 3: use nlsList to call nlme. I use the nlsList object (fit) to fit a variance components model: fit.nlme<-nlme(fit) This works (at least it runs without an error message). To see what syntax is used, I extract the call: > fit.nlme$call nlme.formula(model = photosynthesis ~ photo(irradiance, Am, Q, LCP), data = marouane.data, fixed = list(Am ~ 1, Q ~ 1, LCP ~ 1), random = list(species = c(2.41887429868478, -1.53351210977838, 2.47282942435836, 0, 0, 0)), groups = ~species, start = list( fixed = c(11.7384877502532, 0.109091641284836, 9.91905527125462 )), na.action = na.omit) Question 1: the syntax to the random argument seems wrong! It should be something like: random=(list(Am~1,Q~1,LCP~1)). I have no idea where the 6 numbers (2.41887...) are coming from in the random argument. Can someone explain how the nlme function obtains such an argument for random? The values in fixed are the estimated fixed effects from nlsList. If I follow the instructions in Pinheiro & Bates, the call should be (with a self-starting function): nlme(model = photosynthesis ~ photo(irradiance, Am, Q, LCP), data = marouane.data, fixed = list(Am ~ 1, Q ~ 1, LCP ~ 1), random = list(Am~1,Q~1,LCP~1),groups=~species, na.action=na.omit). When I use this, I get an error message: > nlme(model=photosynthesis~photo(irradiance,Am,Q,LCP),data=marouane.data, + fixed=list(Am~1,Q~1,LCP~1),groups=~species, + random=list(Am~1,Q~1,LCP~1), + na.action = na.omit) Error in nlmeCall[[i]] <- NULL : subscript out of bounds Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] help with nlme function
Hello. I am trying to fit a nonlinear mixed model involving 3 parameters. I have successfully made a self-starting function. getInitial() correctly outputs the initial estimates. I can also use the nlsList with this function to get the separate nonlinear fits by group. However, I get an error message when using the nlme function. Here is the relevent code: fit<-nlsList(photosynthesis~photo(irradiance,Q,Am,LCP)|species/plant/leaf,da ta=marouane.data, + na.action=na.omit) This works, showing that the function "photo" works as a self-starting function. nlme(model=photosynthesis~photo(irradiance,Q,Am,LCP), + data=marouane.data,fixed=Q+Am+LCP~1, + random=Q+Am+LCP~1|species,na.action=na.omit) Error: subscript out of bounds This is what happens when I use the nlme function. I don't know what "subscript out of bounds" means but I assume that there is something wrong with the syntax of my code. The data frame (marouane.data) has dimensions of 1000 and 9. The first three columns give the group structure (species, plants within species, leaves within plants). "species" is a factor while "plants" is coded 1 or 2 and "leaves" is also coded 1 or 2. The other columns give values of measured variables. Could someone explain what I am doing wrong? Alternatively, is there another text besides Pinheiro & Bates that explains the basic syntax of the nlme function? Thanks. Bill Shipley __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] residual df in lmer and simulation results
Hello. Douglas Bates has explained in a previous posting to R why he does not output residual degrees of freedom, F values and probabilities in the mixed model (lmer) function: because the usual degrees of freedom (obs - fixed df -1) are not exact and are really only upper bounds. I am interpreting what he said but I am not a professional statistician, so I might be getting this wrong... Does anyone know of any more recent results, perhaps from simulations, that quantify the degree of bias that using such upper bounds for the demoninator degrees of freedom produces? Is it possible to calculate a lower bounds for such degrees of freedom? Thanks for any help. Bill Shipley North American Editor, Annals of Botany Editor, "Population and Community Biology" series, Springer Publishing Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] http://pages.usherbrooke.ca/jshipley/recherche/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] using split.screen?
Hello. I am having trouble understanding the use of split.screen. I want to divide the device surface first into 4 equal screens: split.screen(figs=c(2,2)) This works. I next want to subdivide each of these 4 screens into 10 subscreens. I do, for the first of these 4 screens: screen(1,new=T) and then: split.screen(figs=c(10,2)) My understanding is that this should split screen 1 (i.e. top right) into 10 screens within its area. However, this does not work and I get the following error message: Error in plot.new() : figure margins too large Does anyone know what I am doing wrong? Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] superimposing histograms con't
Earlier, I posted the following question: I want to superimpose histograms from three populations onto the same graph, changing the shading of the bars for each population. After consulting the help files and the archives I cannot find out how to do this (seemly) simple graph. To be clear, I want - a single x axis (from -3 to 18) - three groups of bars forming the histograms of each population (they will not overlap much, but this is a detail) - the bars from each histogram having different shadings or other visually distinguishing features. Gabor Grothendieck [EMAIL PROTECTED] pointed to some code to to this but I have found another way that works even easier. hist(x[sel1],xlim=c(a,b),ylim=c(A,B)) - this plots the histogram for the first group (indexed by sel1) but with an x axis and a y axis that spans the entire range. par(new=T) - to keep on the same graph hist(x[sel2],main=Null,xlab=NULL,ylab=NULL,axes=F) -superimposes the second histogram par(new=T) - to keep on the same graph hist(x[sel3],main=Null,xlab=NULL,ylab=NULL,axes=F) -superimposes the third histogram Bill Shipley North American Editor, Annals of Botany Editor, "Population and Community Biology" series, Springer Publishing Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] superimposing histograms
I want to superimpose histograms from three populations onto the same graph, changing the shading of the bars for each population. After consulting the help files and the archives I cannot find out how to do this (seemly) simple graph. To be clear, I want - a single x axis (from -3 to 18) - three groups of bars forming the histograms of each population (they will not overlap much, but this is a detail) - the bars from each histogram having different shadings or other visually distinguishing features. Can anyone explain how to do this? Thanks. Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] help with syntax of nlme function
I am having difficulty understanding the syntax of the nlme() function for nonlinear mixed models. The data frame is called Marouane.chlorophyll. The model involves a dependent variable (Absorb) and an independent variable (Ch.surf), which are both numeric variables in the data frame. The data are hierarchically grouped as Espece, Plante nested within Espece, and Feuille nested within Plante (i.e. leaves within plants within species). Espece is defined as a factor and Plante and Feuille are defined as numeric in the data frame. There are two parameters (k1 and b) and I want to fit a model in which these two parameters vary randomly across Espece and Plante (with Feuille being the residual level). Here is how I have specified the call to nlme: > fit<-nlme(model=Absorb~ k1 - (1/(b*Ch.surf)), + data=Marouane.chlorophyll,fixed=k1+b~1, + groups=~Espece/Plante, + start=list(k1=97.8, b=4.68)) However, there is something wrong with this syntax because I get the following error: Problem in names<-: Invalid length for names attribute: structu re(.Data = list(structure(.Data = numeric(0), class Can someone explain what I am doing wrong? Bill Shipley North American Editor, Annals of Botany Editor, "Population and Community Biology" series, Springer Publishing Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] using panel functions in pairs()
Hello. I want to construct a scatterplot matrix in which the above-diagonal
panels represent values from a control group while the below-diagonal values
represent values from a treatment group. I have managed to write primative
panel functions (below) but cannot figure out how to pass arguments to these
panel functions.
Here is the panel function for the upper.panel argument of pairs:
> my.upper.panel
function(x,y,treat.value=1,...)
{points(x[dat$treat==treat.value],y[dat$treat==treat.value],pch=16)}
Here is the panel function for the lower.panel agument of pairs:
> my.lower.panel
function(x,y,treat.value=2,...)
{points(x[dat$treat==treat.value],y[dat$treat==treat.value],pch=1)}
Now, if I simply call: >
pairs(~x+y+z,data=dat,lower.panel=my.lower.panel,upper.panel=my.upper.panel)
I get what I want. "dat" is a data frame with four columns called "treat",
x, y, z and "treat" has values of 1 or 2.
I want to be able to specify the value of "treat.value" in my panel
functions within the call to pairs in cases in which the values of "treat"
are different from 1 or 2. Something like:
pairs(~x+y+z,data=dat,lower.panel=my.lower.panel(treat.value=2),upper.panel=
my.upper.panel(treat.value=1))
However, when I do this, I get the error message:
>
pairs(~x+y+z,data=dat,lower.panel=my.lower.panel(treat.value=2),upper.panel=
my.upper.panel(treat.value=1))
Error in points(x[dat$treat == treat.value], y[dat$treat == treat.value], :
argument "x" is missing, with no default
Question: how can one pass along different values of arguments from within a
panel function such as I have described?
Bill Shipley
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] obtaining residuals from lmer
Hello. I cannot find out how to extract the residuals from a mixed model using the lmer function. Can someone help? Bill Shipley North American Editor, Annals of Botany Editor, "Population and Community Biology" series, Springer Publishing Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] error message explanation for lmer
I am getting the following error message using the lmer function for mixed models with method="Laplace": "nlminb returned message false convergence (8) in: LMEopt(x=mer,value=cv)" Could anyone explain what this means, and how I might overcome (or track down) the problem? Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] documentation for lmer?
Hello. Besides the short paper in the May 2005 edition of R News, I have not found any documentation concerning the lmer function in the lme4 package. Does anyone know of anything more substatial? Also, does anyone know if the nonlinear equivalent of lmer exists yet? Thanks. Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] more documentation on lmer?
Hello. Is there any documentation on the lmer function in the lme4 package beyond what was published in the May 2005 R News (vol.5/1)? As well, has the nonlinear version of lmer appeared yet? Bill Shipley North American Editor, Annals of Botany Editor, "Population and Community Biology" series, Springer Publishing Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] glm binomial with zero proportions
Hello. I must fit a logistic regression to data in the form of proportions, but in which some of the proportions are zero. I therefore cannot use the glm function with a binomial link since the link function is not defined for p=0 or 1. What other solutions are available? Any references to this specific problem (i.e. regression using proportions, of which some are zero) would be welcome. Thanks. Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] help downloading lme4 from CRAN
Hello. I am having trouble downloading the lme4 package from the CRAN
site. The error is:
> local({a <- CRAN.packages()
+ install.packages(select.list(a[,1],,TRUE), .libPaths()[1],
available=a, dependencies=TRUE)})
trying URL `http://cran.r-project.org/bin/windows/contrib/2.0/PACKAGES'
Content type `text/plain; charset=iso-8859-1' length 26129 bytes
opened URL
downloaded 25Kb
also installing the dependencies 'Matrix', 'latticeExtra', 'mlmRev'
trying URL
`http://cran.r-project.org/bin/windows/contrib/2.0/Matrix_0.95-5.zip'
Error in download.file(url, destfile, method, mode = "wb") :
cannot open URL
`http://cran.r-project.org/bin/windows/contrib/2.0/Matrix_0.95-5.zip'
In addition: Warning message:
cannot open: HTTP status was `404 Not Found'
Can someone please explain what I am doing wrong?
Bill Shipley
North American Editor, Annals of Botany
Editor, "Population and Community Biology" series, Springer Publishing
Département de biologie, Université de Sherbrooke,
Sherbrooke (Québec) J1K 2R1 CANADA
[EMAIL PROTECTED]
http://callisto.si.usherb.ca:8080/bshipley/
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] strategies to obtain convergence using nlme
Hello. I am working on an analysis involving the nonlinear mixed model function (nlme) in R. The data consist of measures of carbon fixation by leaves as a function of light intensity and the parametric function (standard in this area because it has a biological interpretation) is a non-rectangular hyperbola. I cannot get the nonlinear mixed model (nlme) function to converge cleanly. I am hoping that others have found strategies for getting starting values that might solve this problem. I have a request (below), some details of the problem, and then a question (at the end of this posting). Request 1: Any documents that you can point to would be welcome, especially if they involve a non-rectangular hyperbola. I have already searched the archives and have read the book by Pinheiro & Bates (2000) and, although they give many useful suggestions, none seem to work in my case. Details: In particular, I can obtain convergence for each group (plant) separately using either nls or nlsList, but cannot get convergence using nlme. My function is a non-rectangular hyperbola with 4 parameters (theta, Am, alpha, Rd): #Nonrectangular hyperbola for photosynthesis # myfunct (1/(2*theta))*( alpha*Irr + Am -sqrt((alpha*Irr+Am)^2-4*alpha*theta*Am*Irr))-Rd I have written a self-starting function (NRhyperbola) to provide starting values. This self-starting function works. When I use this self-starting function with nlsList using a data set with the grouping variable (and after removing a few groups for which there is an insufficient range of the independent variable to provide estimates), I get convergence and reasonable estimates. Similarly if I use nls separately for every group I obtain reasonable estimates of the 4 parameters. If I then use the fixed effect values from nlsList (i.e. the average value of each parameter over the groups) to do a nonlinear mixed model (nlme) using the fixed effects estimates from nlsList, I fail to converge after 50 iterations and get two types of warning: (1): "NaNs produced in: sqrt(.expr10)" - this is because the parameter values drift into regions in which the sqrt in the function is undefined. (2) "Singular precision matrix in level -1, block 1". - I do not know what this warning means. I have tried many different changes in the starting values and none work. Question: what does the second warning mean? Thanks for any help. Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] source of "susbcript out of bounds error" in nmle
A few days ago I posted a question to this discussion group concerning to origin of an error message < subscript out of bounds > while using the nonlinear mixed model (nlme) function in R with a self-starting function. Thanks for those who responded. This posting is to explain what (I think) it causing the error. Pinheiro & Bates (2000, pages 342-347) describe how to construct a self-starting function for use in nls, nlsList and nlme. A self-starting function provides, not only the function to be evaluated, but also estimates of initial values of the parameters that are to be estimated. I have written a self-starting function (NRhyperbola) for a non-rectangular hyperbola with 4 parameters (theta, Am, alpha, Rd). The function works: > getInitial(Photosynthese~NRhyperbola(Irr,theta,Am,alpha,Rd),data=lit.dat a2,na.action=na.omit) theta Am alpha Rd 0.69495045 35.32201003 0.03156335 1.26873011 If one calls the nlsList function with a grouping variable using NRhyperbola one gets the parameter estimates for each group: > nlsList(Photosynthese~NRhyperbola(Irr,theta,Am,alpha,Rd)|Groupe,data=lit .data2,na.action=na.omit) Coefficients: thetaAm alpha Rd 1 0.442150628 3.303208 0.00659723 -3.229906256 2 0.838342670 26.331938 0.05824155 0.239244579 3 0.513979519 8.650714 0.07326543 1.931373980 4 0.304795985 6.080392 0.06161286 1.194595984 5 0.870235949 5.245667 0.04243476 1.873685797 6 0.882610059 7.732813 0.04409574 0.039987285 7 0.901297969 12.655540 0.07025121 0.288532206 8 0.319134655 10.378213 0.01228871 -0.009558434 9 0.173420670 18.799837 0.02182253 0.025902109 10 NANA NA NA 11 NANA NA NA 12 NANA NA NA 13 0.898779206 28.232093 0.06589758 10.130898439 14 NANA NA NA 15 0.588541525 16.076403 0.05198140 4.102262827 16 -0.579787245 85.421300 0.06644187 4.879817577 17 0.538103607 34.046475 0.06338304 3.230418320 18 -0.149304124 46.348766 0.14660061 2.816890625 19 0.001826296 42.038952 0.09758133 0.175960147 20 NANA NA NA Degrees of freedom: 175 total; 115 residual Residual standard error: 1.257301 Parameter values for some groups cannot be estimated because there is insufficient variation in the independent value. Notice that the call to nlsList returns the 4 parameter values for each group. To get the average (fixed effect) values one must use fixef(obj): fixef(nlsList(Photosynthese~NRhyperbola(Irr,theta,Am,alpha,Rd)|Groupe,da ta=lit.data2,na.action=na.omit)) theta Am alpha Rd 0.43627516 23.42282078 0.05883306 1.84600701 When one calls the nonlinear mixed model function (nlme) using this selfstarting function one gets an error stating that the "subscript is out of bounds". However, this error does not occur if one explicitly gives the starting values as those obtained using the fixed effects. It seems as if the nlme function is receiving more than 1 starting value for each of the 4 parameters when it uses the self-starting function. I suspect that the nlme function calls nlsList to get its starting values but does not request the fixed effects (i.e. it uses the output from nlsList rather than the output from fixed(nlsList). I can't get inside the nlme function to verify this. Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] coding nesting in data for nlme example of Wafer data set.
I am trying to understand the proper way to encode the nesting structure for data in the context of nlme, in the specific case of individuals nested within species for which each individual is unique. I have searched through Pinheiro & Bates and also past postings, but without success. Take the Wafer data set which has 2 levels: Wafer (8 values) and Site nested within Wafer (10 values for each value of Wafer) (Pinheiro & Bates book). The data set has Sites coded as values from 1 to 8 for Wafer 1, values 1 to 8 for Wafer 2 etc. Does this mean that the SAME sites are used for each Wafer (i.e. site 1 of water 1 is the same as site 1 of wafer 2)? If different sites were chosen for each of the wafers would one have to code the sites = 1 to 8 for Wafer 1, sites = 9 to 16 for Wafer 2, and so on? In the case of individuals nested within species, each individual is unique - analogous to different sites being chosen for each wafter. Thanks. Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] syntax of nlme with nesting
This may appear too elementary to some on this list, but not to me. My apologies if this is the case. I have mastered the lme function but the nlme function has me stumped. I am attempting to fit a nonlinear mixed model with 4 levels of nesting. I am getting a cryptic error message and do not know what is wrong with the syntax of the call. This is the call: > nlme(Photosynthese~NRhyperbola(Irr,theta,Am,alpha,Rd), + fixed=theta+Am+alpha+Rd~1, + random=theta~1|Reference/Espece/Plante/Groupe, + data=lit.data) NRhyperbola is a self-starting function with one variable (Irr) and four parameters (theta,Am,alpha,Rd). The data set (lit.data) contains Photosynthese (dependent variable) and Irr, as well as the grouping structure, which is Reference, Espece nested in Reference, Plante nested in Espece and Groupe nested in Plante. I want to allow only the parameter theta to vary randomly. I get the following error message: "Error: subscript out of bounds". What does this mean? There are some "Plante" for which there is only one "Groupe" , some "Espece" for which there is only one "Plante" etc. Is this the source of the error? If so, how can one solve this? Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] self starting function for nonlinear least squares.
Following on my posting of this morning, concerning a problem that I am
having constructing a self-starting function for use with nls (and
eventually with nlsList and nlme), the following is the self-starting
function called NRhyperbola:
> NRhyperbola
function (Irr,theta,Am,alpha,Rd)
{
# Am is the maximum gross photosynthetic rate
# Rd is the dark resiration rate (positive value)
# theta is the shape parameter
# alpha is the initial quantum yeild
#
(1/(2*theta))*
(alpha*Irr + Am - sqrt((alpha*Irr + Am)^2 -4*alpha*theta*Am*Irr))- Rd
}
attr(,"initial")
function (mCall, LHS, data)
{
xy <- sortedXyData(mCall[["Irr"]], LHS, data)
if (nrow(xy) < 3)
stop("Too few unique irradiance values")
fit <- smooth.spline(xy[, "x"], xy[, "y"])
alpha <- predict(fit, x = 0, deriv = 1)$y
if(min(xy[,"x"],na.rm=T)>0)delta.for.zero<-alpha*(min(xy[,"x"],na.rm=T)-
0)
if(min(xy[,"x"],na.rm=T)<=0)delta.for.zero<-0
# to correct for minimal values greater that zero
Rd <- (predict(fit, 0)$y - delta.for.zero)
# to correct for maximal values
alpha <- predict(fit, x = max(xy[,"x"],na.rm=T), deriv = 1)$y
if(alpha>0){
Am <- predict(fit,
max(xy[,"x"],na.rm=T))$y+alpha*(3000-max(xy[,"x"],na.rm=T))
Am <- Am - Rd
}
theta1 <- rep(NA, 4)
light <- c(100, 200, 500, 800)
for (i in 1:4) {
y.pred <- predict(fit, light[i])$y
theta1[i] <- ((y.pred + Rd) * (alpha * light[i] + Am) -
alpha * Am * light[i])/(y.pred + Rd)^2
}
theta <- mean(theta1)
Rd<- -1*Rd
value <- c(theta, Am, alpha, Rd)
names(value) <- mCall[c("theta", "Am", "alpha", "Rd")]
value
}
attr(,"class")
[1] "selfStart"
This function was created, following pp. 345-346 of Pinheiro & Bates
(Mixed-effects models in S and S-PLUS) by first creating NRhyperbola:
NRhyperbola<- function (Irr,theta,Am,alpha,Rd)
{
# Am is the maximum gross photosynthetic rate
# Rd is the dark resiration rate (positive value)
# theta is the shape parameter
# alpha is the initial quantum yeild
#
(1/(2*theta))*
(alpha*Irr + Am - sqrt((alpha*Irr + Am)^2 -4*alpha*theta*Am*Irr))- Rd
}
then creating a function
NRhyperbolaInit<-
function (mCall, LHS, data)
{
xy <- sortedXyData(mCall[["Irr"]], LHS, data)
if (nrow(xy) < 3)
stop("Too few unique irradiance values")
fit <- smooth.spline(xy[, "x"], xy[, "y"])
alpha <- predict(fit, x = 0, deriv = 1)$y
if(min(xy[,"x"],na.rm=T)>0)delta.for.zero<-alpha*(min(xy[,"x"],na.rm=T)-
0)
if(min(xy[,"x"],na.rm=T)<=0)delta.for.zero<-0
# to correct for minimal values greater that zero
Rd <- (predict(fit, 0)$y - delta.for.zero)
# to correct for maximal values
alpha <- predict(fit, x = max(xy[,"x"],na.rm=T), deriv = 1)$y
if(alpha>0){
Am <- predict(fit,
max(xy[,"x"],na.rm=T))$y+alpha*(3000-max(xy[,"x"],na.rm=T))
Am <- Am - Rd
}
theta1 <- rep(NA, 4)
light <- c(100, 200, 500, 800)
for (i in 1:4) {
y.pred <- predict(fit, light[i])$y
theta1[i] <- ((y.pred + Rd) * (alpha * light[i] + Am) -
alpha * Am * light[i])/(y.pred + Rd)^2
}
theta <- mean(theta1)
Rd<- -1*Rd
value <- c(theta, Am, alpha, Rd)
names(value) <- mCall[c("theta", "Am", "alpha", "Rd")]
value
}
and then using "selfStart"
NRhyperbola<-selfStart(NRhyperbola,initial=NRhyperbolaInit)
Bill Shipley
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] help with a self-starting function in nonlinear least squares regression.
Hello. I am having a problem setting up a self-starting function for use in nonlinear regression (and eventually in the mixed model version). The function is a non-rectangular hyperbola - called "NRhyperbola" - which is used for fitting leaf photosynthetic rate to light intensity. It has one independent variable (Irr) and four parameters (theta, Am, alpha and Rd). I have created this to act as a self-starting function. The self-starting function seems to work (i.e. when I call "getInitial" it provides the initial values), but I can't get it to work when used within "nls". Here is an example: 1) > getInitial(Photosynthese~NRhyperbola(Irr,theta,Am,alpha,Rd),data=lit.dat a[1:11,]) thetaAm alphaRd 0.5021546914 3.7466359015 0.0005743723 -3.0685671752 So, "getInitial" succeeds in extracting the initial values from the function. 2) > nls(Photosynthese~NRhyperbola(Irr,theta,Am,alpha,Rd),data=lit.data[1:11, ]) Error in eval(expr, envir, enclos) : Object "theta" not found But the call to "nls" does not find the parameter theta, even though the call "getInitial" did find it and returned its initial value. I am working from the Pinheiro & Bates book on mixed-effects models in S and S-PLUS. According to that book (p. 346): "When nls is called without initial values for the parameters and a self-start model function is provided, nls calls getInitial to provide the initial values." Two questions: 1) what am I doing wrong? 2) When a self-Starting model is called from within nlsList or nlme, is getInitial only called one (to get the values ignoring any hierarchical structure in the data) or is it called for each group? Thanks. Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] summary of problem with mCall function.
Last week I posted a question concerning the mCall function, which is used to create self-starting functions and is described in the book by Pinheiro, J.C. and Bates, D.M. (Mixed-effects models in S and S-PLUS). On page 345 one finds the following call: xy<-sortedXyData(mCall[["x"]], LHS,data) It is necessary to replace the "x" in the call to mCall by the actual variable name for the dependent variable. The error message that I was getting was due to the fact that I had called by dependent variable in the data frame (data) by another name than x without changing this in the call to mCall. Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] making self-starting function for nls
Hello. Following pages 342-347 of Pinheiro & Bates, I am trying to
write a self-starting nonlinear function (a non-rectagular hyperbola) to
be used in nonlinear least squares regression (and eventually for a
mixed model). When I use the getInitial function for my self-starting
function I get the following error message:
> getInitial(photo~NRhyperbola(Irr,theta,Am,alpha,Rd),dat)
Error in tapply(y, x, mean, na.rm = TRUE) :
arguments must have same length
Since I do not explicitly call tapply in my function that makes
NRhyperbola a self-starting function (called NRhyperbolaInit, see
below), I assume that the error is coming from within the mCall function
but I can't figure out where or how.
Would someone who has successfully done this be willing to look at my
code and see where the problem arises?
> NRhyperbolaInit
function(mCall,LHS,data)
{
xy<-sortedXyData(mCall[["x"]],LHS,data)
if(nrow(xy)<3){
stop("Too few unique irradiance values")
}
theta<-0.75
Rd<-min(xy[,"y"])
Am<-max(xy[,"y"]) + abs(Rd)
if(sum(xy[,"x"]<50)>3)alpha<-coef(lm(y~x,data=xy,subset=x<50))[2]
if(sum(xy[,"x"]<50)<=3)alpha<-0.07
value<-c(theta,Am,alpha,Rd)
names(value)<-mCall[c("theta","Am","alpha","Rd")]
value
}
Bill Shipley
[EMAIL PROTECTED]
<http://callisto.si.usherb.ca:8080/bshipley/>
http://callisto.si.usherb.ca:8080/bshipley/
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] function for maximum entropy distributions
Hello. I have written a function (Maxent), based on the improved iterative scaling algorithm of Della Pietra et al., that calculates the maximum entropy distribution given moment constraints. That is, it obtains the probability distribution (p) with maximum entropy given constraints linear in p. If anyone is interested in this function, or is willing to test it, I will send it to them. Bill Shipley North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] 1st derivatives using gam.
Hello. Using the smooth.spline function one can obtain the 1st derivatives by specifiying < deriv=1 > in the predict function (i.e. predict(fit,fit$x,deriv=1) where fit is the object created by smooth.spline). Can one obtain the first derivatives using gam() and predict.gam? If so, how? Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] planned means and loess?
Hello. I am trying to figure out how (if?) I can adapt a test of planned means (Dunnett's test?) to a loess fit. Context: The experiment consists of growing plants in hydroponic solution of a period of time (t=0 to crit), after which the nutrient solution concentration is decreased to a new (lower) level. Plants are harvested daily during the course of the experiment and the foliar nitrogen concentration is measured each day from these destructive harvests. The statistical goal is to determine the shortest time following the change in nutrient solution concentration at which one can detect a significant (p=0.05) decrease in the foliar nitrogen concentration. Therefore, I wish to compare the predicted foliar N concentration at t=crit to the foliar N concentrations at time t=crit+1, +2, .. Because the time course of N does not follow any simple parametric function, I am using smoothing splines via loess. Is there any way to adapt a test of planned means (say, Dunnett's test) to this problem? If not, can anyone suggest how this problem can be tackled via loess? Any suggestions are gratefully welcomed. Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] problem with intervals in mixed model
Hello. I am analysing data from a mixed model perspective using the lme() function. The fixed effects model is a quadratic (Y~X+X2) where X2 is the square of X and the data have a 3-level structure. I fitted a series of three models with the same fixed effects but differing in the random effects (only intercept, intercept + X, intercept +X +X2). The anova shows that all three parameters vary significantly (p<0.001) across groups. I have therefore chosen the third model, in which all three parameters vary. When I attempted to obtain the confidence intervals for the correlations between the random components, using: intervals(fit3,which="var-cov") I get the following error message: Problem in intervals.lme(fit3, which = "var..: Cannot get confi dence intervals on var-cov components: Non-positive definite ap proximate variance-covariance I assume that this arises because the correlation between two of the parameters at the 2nd lowest level is -0.998. Can anyone tell me how to deal with this problem? Specifically, 1) how should I interpret such a strong correlation? 2) how can I obtain confidence intervals for these correlations between the random components? Any help is appreciated. Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] anova with gam?
Hello. In SPLUS I am used to comparing nested models in gam using the anova function. When I tried this in R this doesn't work (the error message says that anova() doesn't recognise the gam fit). What must I do to use anova with gam? To be clear, I want to do the following: fFit1<-gam(y~x,.) fit1<-fam(y~s(x),.) anova(fit2,fit1,test="F") Thanks. Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] constraining initial slope in smoother.spline
Hello. I want to fit a smoother spline (or an equivalent local regression method) to a series of data in which the initial value of the 1st derivative (slope) is constrained to a specific value. Is it possible to do this? If so, how? Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] shrinkage estimates in lme
Hello. Slope estimates in lme are shrinkage estimates which pull the OLS slope estimates towards the population estimates, the degree of which depends on the group sample size and the distance between the group-based estimate and the overall population estimate. Although these shrinkage estimates as said to be more precise with respect to the true values, they are also biased. So there is a tradeoff between precision and bias. Are there rules of thumb to help determine when it is better to use the OLS slope estimates and when to use the mixed model (lme) shrinkage estimates? I have 35 groups but the numbers per group vary from over 50 to as low as 4. Thanks for any help. Bill Shipley Subject Matter Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] testing equality of variances across groups in lme?
Hello. I am fitting a two-level mixed model which assumes equality of variance in the lowest-level residuals across groups. The call is: fit3<-lme(CLnNAR~CLnRGR,data=meta.analysis, + na.action="na.omit",random=~1+CLnRGR|study.code) I want to test the assumption of equality of variances across groups at the lowest level. Can someone tell me how to do this? I know that one can test a nested model in which the residuals are given weights, such as: fit3<-lme(CLnNAR~CLnRGR,data=meta.analysis, + na.action="na.omit",random=~1+CLnRGR|study.code, weights=varPower(form~CLnRGR)) This allows the residual variances to increase with the independent variable. However, I have no expectation of how the variances might change, only that they might be different across groups. Bill Shipley Subject Matter Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] cubic spline smoother with heterogeneous variance.
Hello. I want to estimate the predicted values and standard errors of Y=f(t) and its first derivative at each unique value of t using the smooth.spline function. However, the data (plant growth as a function of time) show substantial heterogeneity of variance since the variance of plant mass increases over time. What is the consequence of such heterogeneity of variance in terms of bias in the estimate of the predicted value of Y and its first derivative? I could Ln-transform the data to achieve homogeneity of variance, but this would give me the mean of Ln(Y) at each time (i.e. the mode of Y when back-transformed) and the derivative of Ln(Y) with time (i.e. d(Ln(Y))/dt = dY/YDt), not dY/dt. Can anyone suggest the best strategy for solving this problem? Bill Shipley Subject Matter Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] variance of combinations of means - off topic
Hello, and please excuse this off-topic question, but I have not been able to find an answer elsewhere. Consider a value Z that is calculated using the product (or ratio) of two means X_mean and Y_mean: Z=X_mean*Y_mean. More generally, Z=f(X_mean, Y_mean). The standard error of Z will be a function of the standard errors of the means of X and Y. I want to calculate this se of Z. Can someone direct me to a reference (text book or other) that gives the solution to this *general* problem? Thanks. Bill Shipley [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] inferential test for smoother df?
Hello. In loess regression (or gam with cubic spline smoothers, I think) it is possible to fit models with different numbers of equivalent parameters thus model df and then conduct an inferential test via anova. Is this a valid way of choosing the smoother df? Specifically, I fix a significance level of alpha and then fit a sequence of models with increasing numbers of model df (say 2,3,4 ). I conduct an anova to compare this sequence of models and choose the smoother df as the one at which models fit with further increases do not result in a significant improvement. If this is not an acceptable strategy, what would people recommend beyond using the built in cross-validation criterion? Thanks for any leads. Bill Shipley Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] question about df in linear mixed model function lme
I have looked through Pinheiro and Bates (2000) but have not found an explanation, so I post my question to the group. Pages 90-92 explain how the denominator df are calculated, but not why. Consider a mixed model with two predictor variables (Y, Z) in which Y only varies between groups while Z varies at both levels. The denominator df for the test of the marginal effect of Y is based on the number of groups. The denominator df for the marginal effect of Z is based on the total number of observations. I understand the logic of both. However, why is the denominator df for the interaction term (Y:Z) also based on the total number of observations? Bill Shipley Subject Matter Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] (off topic) article on advantages/disadvantages of types of SS?
Hello. Please excuse this off-topic request, but I know that the question has been debated in summary form on this list a number of times. I would find a paper that lays out the advantages and disadvantages of using different types of SS in the context of unbalanced data in ANOVA, regression and ANCOVA, especially including the use of different types of contrasts and the meaning of the hypotheses that are tested in such cases. Thanks for any leads. Bill Shipley Subject Matter Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] lme function with marginal terms for ANOVA?
Hello. The linear mixed model (lme) function has an anova (anova()) function. This anova function has a type argument that can be marginal. I understood that this was equivalent to the type III (type=ssType3) sum of squares of the ordinary linear model. Is this correct? Specifically, if fit is a lme object, and I specify anova(fit,type=marginal), so I get the type III sum of squares? Thanks Bill Shipley Subject Matter Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] controlling colour in Trellis histogram
Hello. I am sorry for posting a (seemingly) simple question, but I have
just spent 2 hours trying to find the answer, without success. I want
to make a histogram with conditioning on a factor, using Trellis
graphics. However, I do not want any colours (only black and white)
either in the histograms or in the strip. There must be some simple
argument but I cant find it. Here is my code so far:
> histogram(~GRC.SLA|as.factor(species.type),data=group.means,
+ strip=function(...)
+
strip.default(...,style=1,factor.levels=c("Conifers","Trees","Herbs","Mo
nocots")))
As a default, this produces the bars in a pale blue and the strip in
orange-yellow.
Thanks.
Bill Shipley
Subject Matter Editor, Ecology
North American Editor, Annals of Botany
Département de biologie, Université de Sherbrooke,
Sherbrooke (Québec) J1K 2R1 CANADA
[EMAIL PROTECTED]
<http://callisto.si.usherb.ca:8080/bshipley/>
http://callisto.si.usherb.ca:8080/bshipley/
[[alternative HTML version deleted]]
__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] histograms with more than one variable
Hello. I want to plot the distribution of a continuous variable (y) in each of two groups on the same graph as histograms. I suppose one could call this a 2-d histogram? Can this be done in R? Here is a typical data.set: y group 1.2 1 3.3 1 2.4 2 5.7 1 0.2 2 etc. Bill Shipley Subject Matter Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] specifying a complex function in nls
Earlier, I posted a question concerning the syntax of nls() when the
nonlinear function within the call is specified outside of nls().
Sundar showed me how to do this. Now I have a slightly more complex
problem that is analogous to a broken-stick regression but in a
non-linear context.
I want my data to be fit to a non-rectangular hyperbola if the value of
the independent variable (I) is less than some value (Icut, to be
estimated) and to be fit to a constant value (to be estimated) if the
value of the independent variable is greater than or equal to Icut.
That is, if the value of the independent variable is less than Icut,
then fit the non-rectangular hyperbola, otherwise fit a constant.
Here is the function that I wrote. a is the non-rectangular hyperbola
and b is the constant. When I give this to nls() I get an error
message stating:
Error in nlsModel(formula, mf, start) : singular gradient matrix at
initial parameter estimates
> NRHyperbolic.cut
function (Am,alpha,theta,Rd,I,Icut,const)
{
b<-rep(1e6,length(I))
yes<- (I>=Icut)
vec<-as.numeric(yes)
a<-(1/(2*theta))*(alpha*I+Am-sqrt((alpha*I+Am)^2-4*theta*alpha*I*Am))-Rd
b[yes]<-vec[yes]*rep(const,length(I[yes]))
apply(cbind(a,b),1,min)
}
This is a non-rectangular hyperbola (a)
Bill Shipley
Subject Matter Editor, Ecology
North American Editor, Annals of Botany
Département de biologie, Université de Sherbrooke,
Sherbrooke (Québec) J1K 2R1 CANADA
[EMAIL PROTECTED]
<http://callisto.si.usherb.ca:8080/bshipley/>
http://callisto.si.usherb.ca:8080/bshipley/
[[alternative HTML version deleted]]
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] specifying a function in nls
Hello. I am trying to understand the syntax of the nonlinear least
squares function (nls) when the function definition is made outside of
the call. Here is the context.
1. If I specify the following command, it works fine:
> fit2<-nls(
+ A~Am*(1-exp(-alpha*(I-LCP))),data=dat1,
+ start=list(Am=3.6,alpha=0.01,LCP=20))
2. Now, I want to be able to specify the function definition
outside of nls. I do the following:
> Mitscherlich<-function(Am,alpha,I,LCP,...){
Am*(1-exp(-alpha*(I-LCP)))
}
and then:
> fit3<-nls(
+ A~Mitscherlich,data=dat1,
+ start=list(Am=2.7,alpha=0.006,LCP=45))
and I get the error message: Error in lhs - rhs : non-numeric argument
to binary operator
What am I doing wrong?
Thanks.
Bill Shipley
Subject Matter Editor, Ecology
North American Editor, Annals of Botany
Département de biologie, Université de Sherbrooke,
Sherbrooke (Québec) J1K 2R1 CANADA
[EMAIL PROTECTED]
<http://callisto.si.usherb.ca:8080/bshipley/>
http://callisto.si.usherb.ca:8080/bshipley/
[[alternative HTML version deleted]]
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] error message in mle function
I am getting an error message concerning the estimation of confidence intervals when fitting a mixed model and don't know what the problem is, or its solution. Just to provide context: the model is describing the effects of age, exp(age), harvest age, and climate variables on bighorn horn annular length. The data structure is repeated measures (between individuals, within individuals over time). Id is a random effect (there are between 3-11 horn measurements per ram, one horn measurement per age, over the 25 year period in the dataset). The mixed effect results is unable to provide confidence intervals for the fixed and random effects because: of an error in the variance-covariance structure. The error says that the intervals are non-positive definitive. Bill Shipley [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] nls function
Hello. I am trying to fit a non-rectangular hyperbola function to data of photosynthetic rate vs. light intensity. There are 4 parameters that have to be estimated. I find the nls function very difficult to use because it often fails to converge and then gives out cryptic error messages. I have tried playing with the control parameters but this does not always help. Is there another non-linear regression function in R that I might try (other than regression smoothers, which wont give the parameter estimates of the specified function)? Bill Shipley Subject Matter Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] controlling nls errors
Hello. I am using the nonlinear least squares function (nls). The function that I am trying to fit seems to be very sensitive to the starting values and, if these are not chosen properly, the nls function stops and gives an error message: Error in numericDeriv(form[[3]], names(ind), env) : Missing value or an Infinity produced when evaluating the model In addition: Warning message: NaNs produced in: sqrt((quantum.yeild * I + Amax)^2 - 4 * theta * quantum.yeild * I want to embed the nls function within a loop so that if it gives an error message I can automatically adjust the starting value and retry. To do this I have to be able to test if the function has converged within the loop. I need something like an if.error(fit) function that will return true if there has been an error. Does such a thing exist? Bill Shipley Associate Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] specifying partial nesting in lme
Hello. I am fitting a repeated measures model using the mixed model function of R (lme). However, the hierarchical structure is complicated. Each individual sheep is measured a number of times (once per year) over its life (ranging from once to 12 consecutive years). However, this longitudinal study involves many different cohorts and so many different individuals (of different ages) are measured each year. Therefore, individuals measured in the same year are not independent. A typical data set might look like this: Ind age year mass 1 1 1984 2.3 1 2 1985 4.3 1 3 1986 5.5 ... 2 1 1985 3.3 2 2 1986 4.1 ... Can someone explain how to specify this in the random part of the model specification in the lme function of R? If I were to ignore the "year" effect, I would specify: lme(mass~age, random=age|ind) Bill Shipley Associate Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] typeIII SS for lme?
To avoid angry replies, let me first say that I know that the use of Type III sums of squares is controversial, and that some statisticians recommend instead that significance be judged using the non-marginal terms in the ANOVA. However, given that type III SS is also demanded by some is there a function (equivalent to drop1 for lm) to obtain type III sums of squares for mixed models using the lme function? Bill Shipley Associate Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] using pdMAT in the lme function?
Hello. I want to specify a diagonal structure for the covariance matrix of random effects in the lme() function. Here is the call before I specify a diagonal structure: > fit2<-lme(Ln.rgr~I(Ln.nar-log(0.0011)),data=meta.analysis, + random=~1+I(Ln.nar-log(0.0011)|STUDY.CODE,na.action=na.omit) and this works fine. Now, I want to fix the covariance between the between-groups slopes and intercepts to zero. I try do do this using the pdDiag command as follows, but it does not work: > fit2<-lme(Ln.rgr~I(Ln.nar-log(0.0011)),data=meta.analysis, + random=pdDiag(diag(2),~1+I(Ln.nar-log(0.0011))|STUDY.CODE),na.action=na. omit) I get back an error saying that I have zero degrees of freedom. Clearly, the syntax of the command is wrong but I cant figure out why. The data set (meta.analysis) is not defined as a groupedData object. Any help is appreciated. Bill Shipley Associate Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] help specifying error structure in lme
I am having difficulty in properly specifying the error structure for an experiment that is being analysed using the lme() function. Physically, I have 15 plots. Each plot contains 8 subplots. Each subplot is measured at 3 different times. One treatment variable (L), containing 3 levels, is randomly assigned to the whole plots; thus there are 5 replicate plots for each of the 3 levels of L. Two other treatment variables (K, F), each having 2 levels, are randomly assigned to the subplots within each plot; thus there are 2 replicates of each combination of K and F in each of the 15 plots. The fixed formula is y~L*K*F*time. What is the proper formula for the random effects? I had thought that it was random=~1|L/plots/subplots (i.e. subplots nested within plots nested within the L treatments) but this is wrong because, using this formula, I have no degrees of freedom for the L term. Bill Shipley Associate Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] interpreting output of lme
Hello. This is not a technical question but rather an interpretational one, so I apologise if it is not appropriate for this discussion group. I am fitting a mixed model involving two predictor variables: x1 and x2. x1 varies between groups but is constant within groups while x2 varies both between and within groups. I fit two nested models. Model 1: lme(y~x1*x2, random=~1|groups). In this model, in which only the intercepts randomly vary, both main fixed terms and their interaction are highly significant. Model 2: lme(y~x1*x2, random=~1+x1|groups). In this model, in which both the intercepts and the partial slope of x1 randomly vary, only x2 remains significant as a fixed term; both x1 and x1:x2 loose significance (p>0.05). Here is my interpretation, and I would like to know if it is correct: The partial slope of x2 changes as a function of x1 (thus, a real interaction). However, because x1 only varies between groups, the variation in the partial slopes of x2 only occurs between groups. Thus, when I allow the partial slope of x2 to be random, these between-group changes in its partial slope remove the effect of x1 (whose only effect is to change these partial slopes between groups). Does this make sense? If not, how should one explain the differences in the significance of the fixed terms in these two nested models? Bill Shipley Associate Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] help with lme()
Hello. I am trying to determine whether I should be using ML or REML methods to estimate a linear mixed model. In the book by Pinheiro & Bates (Mixed-effects models in S and S-PLUS, page 76) they state that one difference between REML and ML is that « LME models with different fixed-effects structures fit using REML cannot be compared on the basis of their restricted likelihoods. In particular, likelihood ratio tests are not valid under these circumstances. I am not sure what fixed-effects structures means. Does it mean that, as long as the types of contrasts are the same between two models, they ARE comparable, but are NOT comparable if the types of contrasts are changes? Or rather, does it simply mean that one should use t or F tests for the fixed effects, and restrict the likelihood ratio tests to the random effects only if using REML? Bill Shipley Associate Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] help with constrOptim function
Hello. I had previously posted a question concerning the optimization of a nonlinear function conditional on equality constraints. I was pointed towards the contrOptim function. However, I do not understand the syntax of this function with respect to specifying the constraints and so I dont know if it is what I need. The command is: constrOptim(theta, f, grad,ui,ci, ). theta is the initial value of the variables in the function, f and grad are the function to be optimized and its gradient functions (partial deriviatives). Now, ui and ci are stated to be constraints in the help documentation but it is not explained how these constraints are to be specified. To be concrete, here is my objective function: f<- -1*pi*log(pi) i.e. the entropy of a vector of probabilities (pi) The constraints are the first moments of the distribution: C1=sum(k1*pi); C2<-sum(k2*pi); etc. Any help is appreciated. Bill Shipley Associate Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] constrained nonlinear optimisation in R?
Hello. I have searched the archives but have not found anything. I need to solve a constrained optimisation problem for a nonlinear function (maximum entropy formalism). Specifically, Optimise: -1*SUM(p_ilog(p_i)) for a vector p_i of probabilities, conditional on a series of constraints of the form: SUM(T_i*p_i)=k_i for given values of T_i and k_i (these are constraints on expectations). Can this be done in R? Bill Shipley Associate Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] setting up complicated ANOVA in R
Hello. I am about to do a rather complicated analysis and am not sure how to do it. The experiment has a split-plot design and also repeated measures. Both of these complications require one to define an error term and it seems that one cannot specify two such terms. The split-plot command is: aov(y~covariates +A*B+Error(C), data=) where A and B are the fixed effects and C is the plot-level source of error. The repeated measures command is: aov(y~covariates+A*B*time + Error(subject), data=) where subject is the error source for the repeated measures over time. Can these be somehow combined to include a split-plot & repeated-measures design? If not, can I perhaps use a mixed-model analysis with random subjects nested within the whole-plot? Any suggestions or leads are appreciated. Bill Shipley Associate Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] summary of "explaining curious results of aov"
Earlier, I had posted the following question to the group : > Hello. I have come across a curious result that I cannot explain. > Hopefully, someone can explain this. I am doing a 1-way ANOVA with 6 > groups (example: summary(aov(y~A)) with A having 6 levels). I get an > F of 0.899 with 5 and 15 df (p=0.51). I then do the same analysis but > using data only corresponding to groups 5 and 6. This is, of course, > equivalent to a t-test. I now get an F of 142.3 with 1 and 3 degrees > of freedom and a null probability of 0.001. I know that multiple > comparisons changes the model-wise error rate, but even if I did all > 15 comparisons of the 6 groups, the Bonferroni correction to a 5% > alpha is 0.003, yet the Bonferroni correction gives conservative > rejection levels. > > How can such a result occur? Any clues would be helpful. Brian Ripley, Robert Balshaw, Peter Dalgaard and Ted Harding all responded. The answer was basically the same from all: If there is heterogeneity of variances between the groups, and the variances of groups 5 and 6 are smaller than the others, then my result could occur because the average within-group variance over all groups in the general ANOVA is higher than the within-group variance when looking only at groups 5 and 6. Combine this with the very small sample size and unequal group membership. A number of reference books state that ANOVA is fairly robust to moderate degrees of heterogeneity of variance but not what constitutes moderate! Bill Shipley Associate Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] explaining curious result of aov
Hello. I have come across a curious result that I cannot explain. Hopefully, someone can explain this. I am doing a 1-way ANOVA with 6 groups (example: summary(aov(y~A)) with A having 6 levels). I get an F of 0.899 with 5 and 15 df (p=0.51). I then do the same analysis but using data only corresponding to groups 5 and 6. This is, of course, equivalent to a t-test. I now get an F of 142.3 with 1 and 3 degrees of freedom and a null probability of 0.001. I know that multiple comparisons changes the model-wise error rate, but even if I did all 15 comparisons of the 6 groups, the Bonferroni correction to a 5% alpha is 0.003, yet the Bonferroni correction gives conservative rejection levels. How can such a result occur? Any clues would be helpful. Thanks. Bill Shipley Associate Editor, Ecology North American Editor, Annals of Botany Département de biologie, Université de Sherbrooke, Sherbrooke (Québec) J1K 2R1 CANADA [EMAIL PROTECTED] <http://callisto.si.usherb.ca:8080/bshipley/> http://callisto.si.usherb.ca:8080/bshipley/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
