Hi Derek,
have a look at the following made-up example:
f1 - function(x){2*x}
f2 - function(x){-10*x+1}
x-rnorm(10)
x
(x0)*f1(x)
(x=0)*f2(x)
(x0)*f1(x) + (x=0)*f2(x)
Therefore I suggest you should specify the model as follows:
yourNLSmodel - nls(Y ~ (XZ) * f(X,a,b,c) + (X=Z) * g(X,a,d,e),
Hi Bob,
you fitted a log-logistic model specifying the model function LL.4().
If you want to fit a Boltzmann or logistic model then you can use the model
function L.4()
in place of LL.4():
drc.preTBS.2 - drm( Response ~ Dose, data = mydata, fct= L.4( names =
c( Slope, Lower Limit, Upper
Hi Antoon,
now that you mention trying out different methods, maybe you should consider
fitting a
sigmoidal curve to the entire dataset and not only the exponential part (which
constitutes
a very small dataset) as seems to have been the endeavour that initiated the
posting to
R-help.
One
Hi Doug,
you can add the fitted curve using the following general paradigm:
## Plotting the data
plot(p~z)
## Defining grid of z values
## (100 values ensures a smooth curve in your case)
zValues - seq(min(z), max(z), length.out = 100)
## Adding predicted values corresponding to the grid
Hi Hans,
I hope I can resolve your problems below (Marc, thank you very much for cc'ing
me on your
initial response!).
Have a look at the following R lines:
## Fitting the model using drm() (from the latest version)
m1- drm(response ~ dose, data = d, fct = LL.4())
summary(m1)
plot(m1)
##
Yes, you're right: taking logarithms is no longer needed!
Christian
__
R-help@r-project.org 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,
Hi Milton,
you're right that most of the functions in the package drc are suited for
sigmoidal/s-shaped curves defined on the positive axis.
However, there is one exception: the logistic model:
library(drc)
sp.coredep.attr.m1 - drm(preference~dst, data=sp.coredep.attr, fct=L.4())
Hi,
maybe the following line works for you:
beechworth.dt.2 - beechworth.dt[as.numeric(beechworth.dt$Year) %in% 1927:2007,
]
(using as.numeric() to make sure that Year is numeric, maybe not needed?).
Christian
__
R-help@r-project.org mailing list
Hi Thierry,
I think it should be possible to generate such correlated data using a mixed
model approach:
1) Generate pairs of correlated linear predictor values using an ordinary
linear mixed
model setup (for example using rnorm() repeatedly)
2) Back-transform these values using the inverse
Hi Dieter,
the following model assumes a linear relationship between the response
newbone and the
independent variable t with a common intercept equal to 0 and
treatment-dependent slopes:
grd.lme0 - lme(newbone~t:treat-1, data=grd, random=~1|subject)
summary(grd.lme0)
Christian
Hi Patrick,
there exist specialized functionality in R that offer both automated
calculation of
starting values and relatively robust optimization, which can be used with
success in many
common cases of nonlinear regression, also for your data:
library(drc) # on CRAN
## Fitting 3-parameter
Hi Greg,
you can use the extension package 'drc' from CRAN:
conc - c(10.3, 10.8, 11.6, 13.2, 15.8, 20.1) # Exposure concentrations
orign - c(76, 79, 77, 76, 78, 77) # Original number of subjects
ndead - c(16, 22, 40, 69, 78, 77) # Number dead after 96 h
d - data.frame(conc=conc, orign=orign,
Hi Saeed,
one approach is to try out several initial value combinations for a and b.
It often helps to find initial values of the same order of magnitude and of the
same sign
as the final estimates.
To get such initial values, you could linearize the model:
lm(log(q) ~ I(-depth))
and supply
Hi,
doing a search in R gives
help.search(loess)
?loess
Look out for the family argument in the help page.
Christian
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
Hi Mike,
the model you consider is a special case of the four-parameter logistic model
where the
lower and upper asymptotes are fixed at 0.5 and 1, respectively.
Therefore, this (dose-response) model can fitted using the R package 'drc':
library(drc)
xy.m - drm(y~x, fct = L.4(fixed=c(NA,0.5,
Dear Christoph,
using the package 'alr3' it's not difficult!
Have a look at the following example:
## Fitting a Michaelis-Menten model
Puromycin.m1-nls(rate~a*conc/(b+conc), data=Puromycin[1:12,],
start=list(a=200, b=1))
library(alr3)
## Predictions (with standard errors) at concentrations
Hi Göran,
the R package NADA is specifically designed for left-censored data:
http://www.statistik.uni-dortmund.de/useR-2008/slides/Helsel+Lee.pdf
The function cenreg() in NADA is a front end to survreg().
Christian
Göran Broström wrote:
On Fri, Oct 31, 2008 at 2:27 PM, Terry Therneau
Hi Christian,
I believe that the argument var has changed name to weights.
The following lines work for me:
...
sigma - rep(sqrt(.5), nrow(lumley1)) # not nrow=
lme1 - lme(Y1 ~ trt.B + trt.C + trt.D + trt.E, random = ~ 1 | trtpair,
data=lumley1, weights = varConstPower(form=~sigma,
Dear Agnes,
I think your model specification should look like this:
YourModel1 - lmerlmer(y ~ poptype*matingtype + (1|poptype:pop) +
(1|poptype:fam),
data = ...)
The 1 in front of | refers to models that are random intercepts models as
opposed to
general random coefficients models in which
Hi Benoit,
another way of making Petr's point is by looking at the profile log likelihood
function
for b; that is, only estimating the a parameter for a grid of b values:
## Defining mean function for fixed b
lgma - function(b){
function(C0, m, V, a){ (V + b * m * a + C0 * V * b - ((C0 *
Hi Mike,
try:
plot(1:10, main=expression(paste(Delta*i, , 0)))
Christian
__
R-help@r-project.org 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
Hi Peter,
I think one option for what anova could do in the nonlinear case is to report the analysis
of variance (or deviance) table obtained when doing a lack-of-fit test, that is comparing
the nonlinear regression model to an appropriate ANOVA model. This is for example the use
of anova in
Dear Irene,
if you have a vector of estimates of x1, x2, y1, y2 and the
corresponding estimated variance-covariance matrix (accessible through
the method vcov), then one possibility is to use the function
delta.method() in the package alr3 (on CRAN).
This function returns:
1) an estimate of
Hi Philip,
your data are event times because you're monitoring the same trees in
each plot over time, the event being death of a tree.
Therefore methods from survival analysis are more appropriate.
Start out having a look at the package survival, possibly considering
a Cox model with
Hi Steve,
I think you need to use apply() as in the following tiny example:
x - data.frame(response = c(-121,-131,-135))
apply(x, 1, function(response){rnorm(10, mean = response, sd =
rnorm(10, mean = 9.454398, sd = 1.980136))})
Christian
__
Hi,
have a look at ?try which leads to using a construct of the following
form:
myTtest - try(t.test(...))
if (inherits(myTtest, try-error))
{
cat(myTtest)
myTtest - NA
}
Christian
__
R-help@r-project.org mailing list
Hi Jarek,
an alternative approach is to provide more precise starting values!
It pays off to realise that it's possible to find quite good guesses
for some of the parameters in your model function:
t ~ tr+(ts-tr)/((1+(a*h)^n)^(1-(1/n)))
The parameters tr and ts correspond to the
Dear Simone,
you can use the package 'drc' to fit a logistic model, that is a non-linear
regression
model based on the equation c+(d-c)/(1+exp(b(x-e))), to your data (below named
'simone'):
## Fitting a 4-parameter logistic model (also called Boltzmann model)
simone.m1 - drm(size~x,
Dear Andreas,
try:
anova(incub.lme2, type=marginal)
Read more about marginal vs. sequential tests in the help page ?anova.lme.
Christian
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting
Dear Michael,
as to questions 1) and 2):
The following article by A. Zeileis provides more details on the
package 'sandwich':
http://www.jstatsoft.org/v16/i09/
(see also the references).
Use 'model.matrix' to construct a design matrix:
x - runif(10,0,10)
model.matrix(~x)
Christian
Hi again,
in the MLE setting the score function (no expectation taken) is the
estimating function. So for the OLS situation the basic estimating
function is: (in the terminology of Zeileis' paper)
psi(x,y,beta) = (y - x^t beta) x^t
Christian
__
Hi!
I would suggest trying out a few different starting values as a first
unsystematic approach.
For example changing hill=1 to hill=2 results in convergence:
foo.nls-nls(var~Emax*(Dose^hill)/((EC50^hill)+(Dose^hill)),
start=list(Emax=-4,EC50=269,hill=2),trace=T,data=foo)
Christian
Hi!
Use a list structure for all the components you want to have returned
by the function:
return(list(x=x, y=y, prob=prob))
Christian
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting
33 matches
Mail list logo
