Re: [R] regular expressions : extracting numbers
Dear David, does the following work for you? sVec - c(lema, rb 2%, rb 2%, rb 3%, rb 4%, rb 3%, rb 2%,mineuse, rb, rb, rb 12, rb, rj 30%, rb, rb, rb 25%, rb, rb, rb, rj, rb) reVec - regexpr([[:digit:]]+, sVec) # see ?regex for details on '[:digit:]' and '+' substr(sVec ,start = reVec, stop=reVec + attr(reVec, match.length) - 1) # see ?substr for details Christian __ 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.
Re: [R] dose-response on a grid
Hi Bill, have a look at the following artificial example: ## Loading the package 'drc' (on CRAN) library(drc) ## Generating dataset with four dose-response curves finneyx4 - rbind(finney71, finney71, finney71, finney71) ## Generating artificial points (x,y) ## different pairs for each of the 4 curves in the above dataset finneyx4$x - rep(1, 24) finneyx4$y - rep(1:4, c(6, 6, 6, 6)) ## Fitting the two-parameter log-logistic model (logistic regression) m1 - drm(affected/total ~ dose, as.factor(x):as.factor(y), weights = total, data = finneyx4, fct = LL.2(), type = binomial) ## Calculating ED50/LD50 for each location (they are all the same for this dataset) ED(m1, 50) You could try the same approach for your data! Christian __ 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] Research assistant in biostatistics in Copenhagen
The Statistics Group at the Department of Natural Sciences, Faculty of Life Sciences, University of Copenhagen, is looking for a research assistant with R skills. For details see: http://www.matfys.kvl.dk/~torbenm/eu Please note that the deadline for applications is June 22 2007. Christian __ 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.
Re: [R] ED50 from logistic model with interactions
Hi Kate, try looking at the package 'drc' on CRAN and in particular look at the example in the help page for the dataset 'daphnids' (?daphnids). You can obtain arbitrary ED values with approximate standard errors using the function 'ED'. Christian __ 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.
Re: [R] glm() problem
Hi, try using the function 'glm.control' in the first place: glm(n~., data = mDat, family = poisson, control = glm.control(trace = TRUE)) Christian __ 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] Error in nlme with factors in R 2.4.1
Hi, the following R lines work fine in R 2.4.0, but not in R 2.4.1 or any devel versions of R 2.5.0 (see below for details). library(drc) # to load the dataset 'PestSci' library(nlme) ## Setting starting values sv - c(0.43355869, 2.49963220, 0.05861799, 1.73290589, 0.38153146, 0.24316978) ## No error m1 - nlme(SLOPE ~ c + (d-c)/(1+exp(b*(log(DOSE)-log(e, fixed = list(b+c+d+e~1), random = d~1|CURVE, start = sv[c(2:5)], data = PestSci) ## No error m2 - nlme(SLOPE ~ c + (d-c)/(1+exp(b*(log(DOSE)-log(e, fixed = list(b~HERBICIDE, c+d+e~1), random = d~1|CURVE, start = sv[c(1:5)], data = PestSci) ## Error in R 2.4.1! m3 - nlme(SLOPE ~ c + (d-c)/(1+exp(b*(log(DOSE)-log(e, fixed = list(b~factor(HERBICIDE)-1, c~1, d~1, e~factor(HERBICIDE)-1), random = d~1|CURVE, start = sv, data = PestSci) Error in dimnames(data) - dimnames : length of 'dimnames' [1] not equal to array extent An ensuing call to traceback() yields: 7: array(0, c(n, n), list(levs, levs)) 6: contr.treatment(n = 0) 5: do.call(contr[[nm]], list(n = length(levs))) 4: FUN(factor(HERBICIDE)[[1]], ...) 3: lapply(nms, contrMat, contr = contr, data = dataMix) 2: nlme.formula(SLOPE ~ c + (d - c)/(1 + exp(b * (log(DOSE) - log(e, fixed = list(b ~ factor(HERBICIDE) - 1, c ~ 1, d ~ 1, e ~ factor(HERBICIDE) - 1), random = d ~ 1 | CURVE, start = sv, data = PestSci) 1: nlme(SLOPE ~ c + (d - c)/(1 + exp(b * (log(DOSE) - log(e, fixed = list(b ~ factor(HERBICIDE) - 1, c ~ 1, d ~ 1, e ~ factor(HERBICIDE) - 1), random = d ~ 1 | CURVE, start = sv, data = PestSci) Output from sessionInfo() for R 2.4.0 R version 2.4.0 (2006-10-03) i386-pc-mingw32 locale: LC_COLLATE=Danish_Denmark.1252;LC_CTYPE=Danish_Denmark.1252;LC_MONETARY=Danish_Denmark.1252;LC_NUMERIC=C;LC_TIME=Danish_Denmark.1252 attached base packages: [1] methods stats graphics grDevices utils datasets [7] base other attached packages: drc plotrix nlme MASS lattice 1.0-7 2.1-1 3.1-77 7.2-29 0.14-9 Output from sessionInfo() for R 2.4.1 R version 2.4.1 (2006-12-18) i386-pc-mingw32 locale: LC_COLLATE=Danish_Denmark.1252;LC_CTYPE=Danish_Denmark.1252;LC_MONETARY=Danish_Denmark.1252;LC_NUMERIC=C;LC_TIME=Danish_Denmark.1252 attached base packages: [1] stats graphics grDevices utils [5] datasets methods base other attached packages: drc plotrix nlme MASS lattice 1.0-7 2.1-6 3.1-78 7.2-30 0.14-16 Christian __ 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.
Re: [R] Select the last two rows by id group
Hi Lauri, here is a little modification of the solution for retrieving the last row only : score[as.vector(unlist(tapply(rownames(score), score$id, tail, 2))),] giving the last two rows. Replacing 2 by 6 or 10 gives you the last 6 or 10 rows (if they exist). Christian __ 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] nlme with a factor in R 2.4.0beta
Hi, the following R lines work fine in R 2.4.0 alpha (and older R versions), but not in R 2.4.0 beta (details below): library(drc) # to load the dataset 'PestSci' library(nlme) ## Starting values sv - c(0.328919, 1.956121, 0.097547, 1.642436, 0.208924) ## No error m1 - nlme(SLOPE ~ c + (d-c)/(1+exp(b*(log(DOSE)-log(e, fixed = list(b+c+d+e~1), random = d~1|CURVE, start = sv[c(2,3,4,5)], data = PestSci) ## Error: attempt to select more than one element m2 - nlme(SLOPE ~ c + (d-c)/(1+exp(b*(log(DOSE)-log(e, fixed = list(b~HERBICIDE, c+d+e~1), random = d~1|CURVE, start = sv, data = PestSci) Output from sessionInfo() for R 2.4.0 alpha R version 2.4.0 alpha (2006-09-16 r39365) i386-pc-mingw32 locale: LC_COLLATE=Danish_Denmark.1252;LC_CTYPE=Danish_Denmark.1252;LC_MONETARY=Danish_Denmark.1252;LC_NUMERIC=C; LC_TIME=Danish_Denmark.1252 attached base packages: [1] methods stats graphics grDevices utils datasets [7] base other attached packages: nlme drc 3.1-75 1.0-1 Output from sessionInfo() for R 2.4.0 beta R version 2.4.0 beta (2006-09-24 r39497) i386-pc-mingw32 locale: LC_COLLATE=Danish_Denmark.1252;LC_CTYPE=Danish_Denmark.1252;LC_MONETARY=Danish_Denmark.1252;LC_NUMERIC=C; LC_TIME=Danish_Denmark.1252 attached base packages: [1] methods stats graphics grDevices utils datasets [7] base other attached packages: nlme drc 3.1-76 1.0-1 Christian __ 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.
Re: [R] linear terms within a nonlinear model
Hi, the contributed package 'drc' allows specification of non-linear regression models with individual parameter models that include covariates. For an example see section 8 the accompanying paper in J. Statist. Software (http://www.jstatsoft.org/v12/i05/v12i05.pdf). Christian __ 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.
Re: [R] How to request AIC information from lm object?
Hi Michael, use: extractAIC to get AIC from an lm object: y - rnorm(10) extractAIC(lm(y~1)) Christian __ 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
Re: [R] problems while correlating values
Hi. Try compute the correlation for the transposed data frame: cor( t(person.data), use=pairwise.complete.obs ) Assuming that your data frame person.data contains NAs to indicate where no value is available, I think the following computation yields N: (!is.na(person.data)) %*% t(!is.na(person.data)) using the functions is.na, ! (negation), %*% (matrix multiplication) and t (matrix transpose). For details see their help pages. Christian __ 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
Re: [R] How to reference R in papers?
Hi JJ, try the following function in R: citation() Christian __ 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
Re: [R] how to vectorize a matrix
Hi Lorenzo, maybe the following example is of use? a - matrix(1:25,5,5) stack(as.data.frame(a[, c(1,3,5,2,4)])) Note that 'stack' takes a data frame or list as first argument (not a matrix). Therefore the matrix is first converted to a data frame using 'as.data.frame'. Christian __ 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
Re: [R] probit analysis
Hi Jinsong, you can use the package drc on CRAN to fit logit and Weibull models (but not the probit model) with natural mortality/response and/or natural immunity to binomial data (maximum likelihood estimation). To get an idea try: library(drc) ?earthworms Christian __ 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
Re: [R] nls and factor
Hi Manuel,
an alternative to the approach pointed out by Prof. Ripley is to use
the package 'drc' which allows one or more parameters in a non-linear
regression model to depend on a factor.
You will need the latest version available at www.bioassay.dk (an older
version is available on CRAN).
Bconc - runif(30,0.1,10)
Aconc - runif(30,0.1,10)
At - runif(30,1,30)
Bt - runif(30,1,30)
BKm - 1
AKm - 0
EBoth - -0.41
yB - 100*exp(EBoth*Bt)*Bconc/(BKm+Bconc)+rnorm(30,0,1)
yA - 75*exp(EBoth*At)*Aconc/(AKm+Aconc)+rnorm(30,0,1)
## Constructing the dataset
tvalues - c(At, Bt)
conc - c(Aconc, Bconc)
yfactor - factor(rep(c(A, B), c(30, 30)))
y - c(yA, yB)
## Defining the non-linear function
twodimMM - function()
{
fct - function(x, parm)
{
parm[, 1]*exp(parm[, 2]*x[, 1])*x[, 2]/(parm[, 3] + x[, 2])
}
ssfct - function(dataset) # a very simple self starter function
{
return(c(90, -0.5, 0.8))
}
names - c(lev, Ev, km) # parameter names
return(list(fct=fct, ssfct=ssfct, names=names))
}
library(drc)
model1 - multdrc(y~cbind(tvalues,conc), yfactor, fct=twodimMM())
summary(model1)
## Collapsing the Ev parameters into a single parameter
model2 - multdrc(y~cbind(tvalues,conc), yfactor,
collapse=data.frame(yfactor,1,yfactor), fct=twodimMM())
anova(model2, model1) # model reduction is insignificant
summary(model2)
Christian
__
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
Re: [R] SE estimates for treatment groups from nlme
Hi Katie, maybe the easiest solution is to create a new factor that corresponds to the combinations of the three factors A, B and C. A quick and dirty way to create such a factor is: ABC - factor(paste(A, x, B, x, C, sep = )) ABC and then fit the model using the variable ABC instead of A*B*C. Christian __ 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
Re: [R] How to get around heteroscedasticity with non-linear least squares in R?
Hi Quin, the package 'drc' on CRAN deals with modelling dose-response curves. Moreover it allows adjustment for heterogeneity by means of transformation (Box-Cox transformation) modelling the variance as a power of the mean. See the package documentation for more features. Christian __ 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
Re: [R] calculating IC50
Hi! Yes, the 'drc' package can be used to obtain IC50 or any other ICx value for several, commonly used dose-response models. The vignette is more up-to-date than the article in JSS (which dates back to the start of 2005). Christian Liaw, Andy wrote: Perhaps also of interest is the `drc' package on CRAN. There's an article in JSS describing it (probably the same as the package vignette--- haven't checked). Andy From: Prof Brian Ripley On Thu, 2 Feb 2006 [EMAIL PROTECTED] wrote: I was wondering if there is an R-package to automatically calculate the IC50 value (concentration of a substrance that inhibits cell growth to 50%) for some measurements. Function dose.p in recommended package MASS. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ 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@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@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
Re: [R] suggestions for nls error: false convergence
Hi Spencer.
When using 'optim' and the first try fails you could:
1) try some other methods: Nelder-Mead, BFGS, ...
2) increase the maximum number of iterations (argument maxit in the control
list)
3) specify the argument parscale in the control list, in order to have all
parameters of same magnitude during
optimisation (this is useful if the parameters are suspected to be of different
magnitudes).
Using the default method (Nelder-Mead) with maxit=1000 results in convergence,
and essentially the same estimates are
obtained if you use the method BFGS and set maxit=1000 and parscale=c(277, 100,
101, 10) (the initial starting values):
x - 1:100
y - c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,1,1,1,2,2,2,2,2,3,4,4,4,5,
5,5,5,6,6,6,6,6,8,8,9,9,10,13,14,16,19,21,
24,28,33,40,42,44,50,54,69,70,93,96,110,127,127,141,157,169,
178,187,206,216,227,236,238,244,246,250,255,255,257,260,261,262,266,268,
268,270,272,272,272,273,275,275,275,276)
func2 - function( par,y, x, rescale ) {
par - rescale*par
a = par[1]
m = par[2]
n = par[3]
tau = par[4]
y. - a * (1+m*exp(-x/tau)) / (1+n*exp(-x/tau))
sum((y-y.)2)
}
est.no2 - optim(c(277, 100, 101, 10), func2, hessian=TRUE, y=y, x=x,
rescale=1, control=list(maxit=1000))
est.no3 - optim(c(277, 100, 101, 10), func2, hessian=TRUE, method=BFGS,
y=y, x=x, rescale=1,
control=list(maxit=1000, parscale=c(277, 100, 101, 10)))
The optimisation in the package 'drc' uses BFGS with the maxit and parscale
arguments specified.
Best wishes
Christian
__
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
Re: [R] suggestions for nls error: false convergence
Hi.
An alternative is to use the package 'drc' on CRAN to fit your data!
x - 1:100
y - c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,1,1,1,2,2,2,2,2,3,4,4,4,5,
5,5,5,6,6,6,6,6,8,8,9,9,10,13,14,16,19,21,
24,28,33,40,42,44,50,54,69,70,93,96,110,127,127,141,157,169,
178,187,206,216,227,236,238,244,246,250,255,255,257,260,261,262,266,268,
268,270,272,272,272,273,275,275,275,276)
## Defining the function (in a bit different notation)
logi - function(dose, parm){parm[, 1] * (1+parm[, 2]*exp(-dose/parm[, 3])) /
(1+parm[, 4]*exp(-dose/parm[, 3]))}
## Fitting the function to the data (see ?multdrc for details)
library(drc)
model - multdrc(y~x, fct = list(logi, NULL, c(a, m, tau, n)),
startVal=c(277, 100, 10, 101))
summary(model)
plot(model, log=)
Hope this helps?
Christian
__
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
Re: [R] retrieving p-values in lm
Hi Patrick, try: lm.res.2$coefficients which I found by looking at the content of the function 'summary.lm'. Christian __ 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
Re: [R] language settings
Hi Ariel, ok, you want to change the language? Right click on the R icon on the desktop, choose Properties. In the field Destination you add at the end of the line: language=it (Italian) orlanguage=fr (French) Then the line should look somewhat like: (for Italian) C:\Programs\r\R-2.2.0\bin\Rgui.exe language=it Now you can start R. Christian __ 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
Re: [R] Non linear modeling
Hi Angelo, have a look at the following example which uses 'gls' in the nlme package. library(nlme) x - runif(100, 0, 1) y - x + exp(4*x)*rnorm(100, 0, 2) gls(y~x, correlation = varExp(form=~x)) For details see ?gls and ?varExp. Christian __ 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
Re: [R] fitting weibull distribution
Hi. I think there may be one or more zeros in your data set, causing the problem: x - rgamma(100) fitdistr(x, weibull) fitdistr(c(x,0), weibull) Maybe you should omit the zeros. Christian __ [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
Re: [R] lme parameterization question
Hi, try something like: lme(y~w,random=list(~1|year,~1+w|site)) Christian - Original Message - From: John Fieberg [EMAIL PROTECTED] To: [EMAIL PROTECTED] Sent: Wednesday, April 02, 2003 10:36 PM Subject: [R] lme parameterization question Hi, I am trying to parameterize the following mixed model (following Piepho and Ogutu 2002), to test for a trend over time, using multiple sites: y[ij]=mu+b[j]+a[i]+w[j]*(beta +t[i])+c[ij] where: y[ij]= a response variable at site i and year j mu = fixed intercept Beta=fixed slope w[j]=constant representing the jth year (covariate) b[j]=random effect of jth year, iid N(0,sigma2[b]) a[i]=random effect of the ith site, iid N(0, sigma2[a]) t[i]=random effect of ith site, iid N(0, sigma2[t]) c[ij]=random error associated with ith site and jth year __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
