Re: [R] how do i use the get function to obtain an element from alist...
To Mark Leeds and the others, thank you for solving my problem, and so quickly, JM __ 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] how do i use the get function to obtain an element from a list...
my problem can be explained with the following example: x - 1:12 y - 13:24 a - data.frame(x = x, y = y) ## if i write a$x ## it returns [1] 1 2 3 4 5 6 7 8 9 10 11 12 ## but the function get doesn't recognize a$x. Instead it produces the following error: get(a$x) Error in get(x, envir, mode, inherits) : variable a$x was not found i intend to do it inside a loop, using a new object (and hence, a new name) for each iteration (i.e., instead of a$x, it would be a$1, a$2, a$3, and so on, for a million times). i would greatly appreciate it if someone could help me on this issue, thanks in advance, Juan Manuel Barreneche, Zoología de Vertebrados, Facultad de Ciencias, UDELAR, Uruguay. __ 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] Convert string to list?
Let's say I have the following string: str - P = 0.0, T = 0.0, Q = 0.0 I'd like to find a function that generates the following object from 'str'. list(P = 0.0, T = 0.0, Q = 0.0) Thanks! -- http://mutualism.williams.edu __ 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] Alternative to xyplot()?
Hello list, Is anyone aware of a non-lattice-based alternative to xyplot()? Thanks! Manuel -- http://mutualism.williams.edu __ 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] Alternative to xyplot()?
On Tue, 2007-07-17 at 22:18 +0100, Gavin Simpson wrote: On Tue, 2007-07-17 at 16:39 -0400, Manuel Morales wrote: Hello list, Is anyone aware of a non-lattice-based alternative to xyplot()? x - rnorm(20) y - rnorm(20) plot(x, y) ? If you mean some specific aspect of xyplot(), you'll have to tell us what this is. HTH G Sorry. I was thinking of the groups functionality, as illustrated below: grps-rep(c(1:3),10) x-rep(c(1:10),3) y-x+grps+rnorm(30) library(lattice) xyplot(y~x,group=grps, type=c(r,p)) __ 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] running Rcmdr
Yes, of course, i know it, and i have installed all packages i need, but this is not what i am searching. I want to do that only with a command line without have to type only first R and after library(Rcmdr); i only want to type something like R library(Rcmdr) but i know it´s doesn´t work and i´m looking for something like that. I hope you have already understood my question. Thank you Felipe Carrillo escribió: If you have already installed Rcmdr from the internet (Cran site) then when you open R window if you type library(Rcmdr) you should see the R commander window and use it just like you would use R command window. */Manuel [EMAIL PROTECTED]/* wrote: Thanks for your answer, but i think i dont do correctly my question. I need the command line to run Rmdr, like that: R Rcmdr or R loadRcmdr.R where loadRcmdr has library(Rcmdr). or something like that. I tried the last example, but when Rcmdr is executed, later it is closed. About RProfile.site, i dont know what i have to change. If you think its useful to me, please explain me. Thanks. Felipe Carrillo escribió: First, you need to install the Rcmdr packages and then in the R Command window or Tinn-R window type library(Rcmdr) without the quotation marks. In addition, if you want Rcmdr to start automatically everytime you start R, go to the following path C:\Program Files\R\R-2.5.0\etc\RProfile.site if you installed R in program files, otherwise follow your own path and type the same command library(Rcmdr) and the R commander window should pop up everytime you start R */Manuel /* wrote: Hi to all, i want to know how can run Rcmdr automatically , or how to load a library in the call of R greetings __ 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. Choose the right car based on your needs. Check out Yahoo! Autos new Car Finder tool. __ 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. Park yourself in front of a world of choices in alternative vehicles. Visit the Yahoo! Auto Green Center. http://us.rd.yahoo.com/evt=48246/*http://autos.yahoo.com/green_center/;_ylc=X3oDMTE5cDF2bXZzBF9TAzk3MTA3MDc2BHNlYwNtYWlsdGFncwRzbGsDZ3JlZW4tY2VudGVy __ 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] degrees of freedom in lme
On Mon, 2007-06-25 at 13:15 -0400, Doran, Harold wrote: This is such a common question that it has a an FAQ-like response from Doug Bates. Google lmer p-values and all that to find the response. Isn't this a different question, though, since Jean-Baptiste is using nlme. Details on the calculation of DF in nlme can be found in chapter 4 of the book by Pinheiro and Bates Mixed Effects Models in S and S-PLUS. Using the formula provided, I get denDF of 10 for level 1 and 32 for level 2. I'm not sure why lme is using the denDF estimated at level 2 in this example ... -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Jean-Baptiste Ferdy Sent: Monday, June 25, 2007 12:26 PM To: r-help@stat.math.ethz.ch Subject: [R] degrees of freedom in lme Dear all, I am starting to use the lme package (and plan to teach a course based on it next semester...). To understand what lme is doing precisely, I used balanced datasets described in Pinheiro and Bates and tried to compare the lme outputs to that of aov. Here is what I obtained: data(Machines) summary(aov(score~Machine+Error(Worker/Machine),data=Machines)) Error: Worker Df Sum Sq Mean Sq F value Pr(F) Residuals 5 1241.89 248.38 Error: Worker:Machine Df Sum Sq Mean Sq F valuePr(F) Machine2 1755.26 877.63 20.576 0.0002855 *** Residuals 10 426.53 42.65 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Error: Within Df Sum Sq Mean Sq F value Pr(F) Residuals 36 33.287 0.925 anova(lme(fixed=score~Machine,random=~1|Worker/Machine,data=Machines)) numDF denDF F-value p-value (Intercept) 136 773.5709 .0001 Machine 210 20.5762 3e-04 No problem here: the results are essentially the same, which is expected. Now I turn to an ANCOVA with a random grouping factor. data(Orthodont) OrthoFem - Orthodont[Orthodont$Sex==Female,]; summary(aov(distance~age+Error(Subject/age),data=OrthoFem)) Error: Subject Df Sum Sq Mean Sq F value Pr(F) Residuals 10 177.227 17.723 Error: Subject:age Df Sum Sq Mean Sq F valuePr(F) age1 50.592 50.592 52.452 2.783e-05 *** Residuals 10 9.645 0.965 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Error: Within Df Sum Sq Mean Sq F value Pr(F) Residuals 22 9.8250 0.4466 anova(lme(fixed=distance~age,random=~1+age|Subject,data=OrthoFem)) numDF denDF F-value p-value (Intercept) 132 1269.7764 .0001 age 132 52.4517 .0001 This time the F values are (almost) identical, the numerator degrees of freedom are the same, but the denominator degrees of freedom are very different (10 for aov vs. 32 for lme). I understand that there is an issue with the estimation of that number, but I would naively expect the number given by lme to be close to that provided by aov is the case of a balanced dataset. That's obviously not true in the case of an ANCOVA... But why?? And how should I interpret the F-test given by anova.lme? Thanks in advance for your help ! -- Jean-Baptiste Ferdy Institut des Sciences de l'Évolution de Montpellier CNRS UMR 5554 Université Montpellier 2 34 095 Montpellier cedex 05 tel. +33 (0)4 67 14 42 27 fax +33 (0)4 67 14 36 22 __ 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@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. -- Manuel A. Morales http://mutualism.williams.edu __ 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] running Rcmdr
Thanks for your answer, but i think i dont do correctly my question. I need the command line to run Rmdr, like that: R Rcmdr or R loadRcmdr.R where loadRcmdr has library(Rcmdr). or something like that. I tried the last example, but when Rcmdr is executed, later it is closed. About RProfile.site, i dont know what i have to change. If you think its useful to me, please explain me. Thanks. Felipe Carrillo escribió: First, you need to install the Rcmdr packages and then in the R Command window or Tinn-R window type library(Rcmdr) without the quotation marks. In addition, if you want Rcmdr to start automatically everytime you start R, go to the following path C:\Program Files\R\R-2.5.0\etc\RProfile.site if you installed R in program files, otherwise follow your own path and type the same command library(Rcmdr) and the R commander window should pop up everytime you start R */Manuel [EMAIL PROTECTED]/* wrote: Hi to all, i want to know how can run Rcmdr automatically , or how to load a library in the call of R greetings __ 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. Choose the right car based on your needs. Check out Yahoo! Autos new Car Finder tool. http://us.rd.yahoo.com/evt=48518/*http://autos.yahoo.com/carfinder/;_ylc=X3oDMTE3NWsyMDd2BF9TAzk3MTA3MDc2BHNlYwNtYWlsdGFncwRzbGsDY2FyLWZpbmRlcg--%20 __ 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] running Rcmdr
Hi to all, i want to know how can run Rcmdr automatically , or how to load a library in the call of R greetings __ 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] Grid-drawing in plot.its OR how to influence axTicks
Dear List, I am ploting an its object with user-specified ticks on the x axis, but cannot figure out how to adjust the plotted grid to these ticks. temp - its(1:3,dates=as.POSIXct(strptime(c(1999-12-31 01:00:00,2000-01-01 02:00:00,2000-01-10 02:00:00),format=%Y-%m-%d %X))) plot(temp) Gives the grid alined with the ticks on the y-axis but not on the x-axis d-as.POSIXct(round(range(dates(temp)),days)) atTicks-seq(d[1], d[2], by=day) plot(temp, at=atTicks) Gives daily ticks but does not affect the grid axTicks(1) [1] 94660 94680 94700 94720 94740 unclass(atTicks) [1] 946594800 946681200 946767600 946854000 946940400 947026800 947113200 [8] 947199600 947286000 947372400 947458800 attr(,tzone) [1] axTicks(1) gives where grid is drawn, unclass(atTicks) gives where the Ticks are (and where I would like the grid is drawn). plot(temp, at=atTicks, xaxp=c(unclass(range(atTicks)),length(atTicks)-1)) Does not help and my question simply is how to make axTicks(1) to be unclass(atTicks). Many thanks for any hints. Manuel __ 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] Problem adjusting x-labels with bargraphCI
On Tue, 2007-03-20 at 13:15 -0400, Vera, Pedro L. wrote: Hello: I'm having quite a bit of difficulty adjusting the x-labels using bargraphCI. I've tried using text and srt=45 to rotate the labels or mtext for 2 lines to break up the labels. However, using either method, I cannot line up the labels with the midpoints of the bars (they line up with some sort of tick mark that is off the midpoint of the bars). Any suggestions would be greatly appreciated. Regards Pedro L. Vera, Ph.D. University of South Florida, Dept of Surgery Assuming that you are asking about bargraph.CI from the sciplot package, and that you are trying to get x-axis tick-labels that are perpendicular to the x-axis, specifying las=2 in the call to bargraph.CI should work. If you want more control, you can assign the output of bargraph.CI to an object (a list that contains the x-values of the bars, summary stats, and CIs). This will allow you to position labels at the x-values of the plotted bars. For example: test - bargraph.CI(x.factor = dose, response = len, data = ToothGrowth, xaxt=n) axis(side=1,at=test$xvals,labels=c(bar1,bar2,bar3),las=2) -- Manuel A. Morales http://mutualism.williams.edu __ 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] lme4/lmer: P-Values from mcmc samples or chi2-tests?
On Tue, 2007-02-13 at 14:45 +0100, Christoph Scherber wrote: Dear R users, I have now tried out several options of obtaining p-values for (quasi)poisson lmer models, including Markov-chain Monte Carlo sampling and single-term deletions with subsequent chi-square tests (although I am aware that the latter may be problematic). However, I encountered several problems that can be classified as (1) the quasipoisson lmer model does not give p-values when summary() is called (see below) (2) dependence on the size of the mcmc sample (3) lack of correspondence between different p-value estimates. How can I proceed, left with these uncertainties in the estimations of the p-values? Below is the corresponding R code with some output so that you can see it all for your own: ## m1-lmer(number_pollinators~logpatch+loghab+landscape_diversity+(1|site),primula,poisson,method=ML) m2-lmer(number_pollinators~logpatch+loghab+landscape_diversity+(1|site),primula,quasipoisson,method=ML) summary(m1);summary(m2) #m1: [...] Fixed effects: Estimate Std. Error z value Pr(|z|) (Intercept) -0.403020.57403 -0.7021 0.48262 logpatch 0.109150.04111 2.6552 0.00793 ** loghab 0.087500.06128 1.4279 0.15331 landscape_diversity 1.023380.40604 2.5204 0.01172 * snip The p-values from mcmc are: ## markov1=mcmcsamp(m2,5000) HPDinterval(markov1) lower upper (Intercept) -1.394287660 0.6023229 logpatch 0.031154910 0.1906861 loghab0.002961281 0.2165285 landscape_diversity 0.245623183 1.6442544 log(site.(In)) -41.156007604 -1.6993996 attr(,Probability) [1] 0.95 ## mcmcpvalue(as.matrix(markov1[,1])) #i.e. the p value for the intercept [1] 0.3668 mcmcpvalue(as.matrix(markov1[,2])) #i.e. the p-value for logpatch [1] 0.004 mcmcpvalue(as.matrix(markov1[,3])) #i.e. the p-value for loghab [1] 0.0598 mcmcpvalue(as.matrix(markov1[,4])) #i.e. the p-value for landscape.div [1] 0.0074 If one runs the mcmcsamp function for, say, 50,000 runs, the p-values are slightly different (not shown here). ##here are the p-values summarized in tabular form: snip [MCMC] logpatch 0.004 loghab0.0598 landscape_diversity 0.0074 snip [single-term deletions] logpatch 0.007106 loghab0.1704 landscape_diversity 0.01276 snip To summarize, at least for quasipoisson models, the p-values obtained from mcmcpvalue() are quite different from those obtained using single-term deletions followed by a chisquare test. Especially in the case of loghab, the difference is so huge that one could tend to interpret one of the p-values as marginally significant. The P-values from the MCMC chain look suspiciously like 1/2 the others. Are you effectively doing a 1-tailed test in your MCMC comparison? -- Manuel A. Morales http://mutualism.williams.edu signature.asc Description: This is a digitally signed message part __ 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] R vs Matlab {Re: R in Industry}
On Thu, 2007-02-08 at 16:53 +0100, Martin Maechler wrote:
Albr == Albrecht, Dr Stefan (AZ Private Equity Partner) [EMAIL
PROTECTED]
on Thu, 8 Feb 2007 16:38:18 +0100 writes:
snip
Albr And, I was very astonished to realise, Matlab is very, very much
faster
Albr with simple for loops, which would speed up simulations
considerably.
Can you give some evidence for this statement, please?
At the moment, I'd bet that you use forgot to pre-allocate a
result array in R and do something like the notorious horrible (:-)
1-dimensional
r - NULL
for(i in 1:1) {
r[i] - verycomplicatedsimulation(i)
}
instead of the correct
r - numeric(1)
for(i in 1:1) {
r[i] - verycomplicatedsimulation(i)
}
Would a similar speed issue arise for the construction:
r - vector()
...
--
Manuel A. Morales
http://mutualism.williams.edu
__
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] getting a new factor
On Sun, 2006-12-31 at 00:25 +0100, Ricardo Rodríguez wrote: Hi all, Given a data frame as... head(veg) genus species trophia type geo zone importance 1 Sphagnum subsecundum MA En100 2 Sphagnum denticulatum MA En200 3 Molinia caerulea MA En300 4 Sphagnumflexuosum MA En400 5 Juncus squarrosus MA En500 6Carex echinata MA En600 I do need a new one with a new factor constructed from a concatenation of the two first columns, genus and species. I am guessing this must be a simple matter while working with R, but I am stuck at this point. veg$genus.species - paste(veg$genus,veg$species) will add a new column combining genus and species. -- Manuel A. Morales http://mutualism.williams.edu signature.asc Description: This is a digitally signed message part __ 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] Wombling
Hello, Does anyone know if there is any package in R to perform wombling (sensu Barbujani et al. 1990) or a related method? RSiteSearch only gives 2 results with wombling and I am not aware of any other name for the method. Some weeks ago there was a related question on the list but there was no response... I acknowledge that programming it wouldn't be so difficult (only needs to apply solve() over a grid), but I am not a good programmer and I am a little bit concerned with efficiency. Thank you for any information, -- José-Manuel Blanco-Moreno --- Department of Biological Sciences University of Alberta CW-405 Biological Sciences Bldg. Edmonton, Alberta, Canada T6G 2E9 --- Phone: (1) 780 492 3289 Fax: (1) 780 492 9457 __ 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] Force square crosstabulation
Hello list members, I'm looking for a way to force the results of a crosstabulation to be square - that is, to include 0 values. For example: table(letters[1:4],letters[c(1:3,3)]) yields: a b c a 1 0 0 b 0 1 0 c 0 0 1 d 0 0 1 I would like to return: a b c d a 1 0 0 0 b 0 1 0 0 c 0 0 1 0 d 0 0 1 0 Any suggestions? Thanks! -- Manuel A. Morales http://mutualism.williams.edu signature.asc Description: This is a digitally signed message part __ 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] s.e. on interaction plots
On Fri, 2006-11-17 at 05:21 -0500, Chuck Cleland wrote: Murray Pung wrote: Is it possible to add standard error bars to the means on interaction plots? Not with interaction.plot(), as far as I know. You might consider the effects package by John Fox. For example: Another possibility is the function lineplot.CI from the package sciplot (a wrapper that adds SE bars to the output from interaction plots). For example: library(sciplot) data(ToothGrowth) lineplot.CI(dose, len, group = supp, data = ToothGrowth) -- Manuel A. Morales http://mutualism.williams.edu __ 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] Announcement: sciplot (includes functions for graphs w/ error bars)
The package sciplot is now available for download from CRAN. This package includes a collection of functions that create graphs with error bars for data collected from one-way or higher factorial designs, as well as a function to plot bifurcation diagrams resulting from analysis with XPPAUTO. The functions in this package replicate some of the functionality of plotmeans() from the package gplots, with differences in the treatment of two-way and higher designs. Example graphs can be seen at: http://mutualism.williams.edu/sciplot -- Manuel A. Morales http://mutualism.williams.edu signature.asc Description: This is a digitally signed message part __ 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 data and subset arguments in user-defined functions
Dear list,
A while ago, I posted a question asking how to use data or subset
arguments in a user-defined function. Duncan Murdoch suggested the
following solution in the context of a data argument:
data - data.frame(a=c(1:10),b=c(1:10))
eg.fn - function(expr, data) {
x - eval(substitute(expr), envir=data)
return(mean(x))
}
eg.fn(a,data)
I've tried various approaches to add a subset argument to the example
above, but no luck. I'm looking for something like the following (but
that works!)
eg.fn2 - function(expr, data, subset) {
data - subset(data,subset)
x - eval(substitute(expr), envir=data)
return(mean(x))
}
eg.fn2(a, data, subset=a3)
This returns the error:
Error in eg.fn2(a, data, subset = a 3) :
object a not found
Any suggestions?
Thanks!
--
Manuel A. Morales
http://mutualism.williams.edu
signature.asc
Description: This is a digitally signed message part
__
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] Fwd: rarefy a matrix of counts
On Wed, 2006-10-11 at 14:25 -0400, Brian Frappier wrote: I tried all of the approaches below. the problem with: x - data.frame(matrix(NA,100,3)) for (i in 2:ncol(DF)) x[,i-1] - sample(rep(DF[,1], DF[,i]),100) if you want result in data frame or x-vector(list, 3) for (i in 2:ncol(DF)) x[[,i-1]] - sample(rep(DF[,1], DF[,i]),100) is that this code still samples the rows, not the elements, i.e. returns 100 or 300 in the matrix cells instead of red or a matrix of counts by color (object type) like: x1x2 x3 red 32 560 gr6895 40 sum 100 100 100 It looks like Tony is right: sampling without replacement requires listing of all elements to be sampled. snip How about the following approach which generates a new sample using the rMultinom function from Hmisc. library(Hmisc) data - matrix(c(400, 300, 2500, 100, 25, 200, 300, 1000, 500), nrow=3, byrow=TRUE) col.sums - apply(data,2,sum) probs - t(data)/col.sums w - rMultinom(probs,100) apply(w, 1, table) Note that I replaced the zero in your example data set with 25 because the table function doesn't seem to output the results nicely when there are zero values. HTH, Manuel On 10/11/06, Tony Plate [EMAIL PROTECTED] wrote: Here's a way using apply(), and the prob= argument of sample(): df - data.frame(sample1=c(red=400,green=100,black=300), sample2=c(300,0,1000), sample3=c(2500,200,500)) df sample1 sample2 sample3 red 400 3002500 green 100 0 200 black 3001000 500 set.seed(1) apply(df, 2, function(counts) sample(seq(along=counts), rep=T, size=7, prob=counts)) sample1 sample2 sample3 [1,] 1 3 1 [2,] 1 3 1 [3,] 3 3 1 [4,] 2 3 2 [5,] 1 3 1 [6,] 2 3 1 [7,] 2 3 3 Note that this does sampling WITH replacement. AFAIK, sampling without replacement requires enumerating the entire population to be sampled from. I.e., you cannot do sample(1:3, prob=1:3, rep=F, size=4) instead of sample(c(1,2,2,3,3,3), rep=F, size=4) -- Tony Plate From reading ?sample, I was a little unclear on whether sampling without replacement could work Petr Pikal wrote: Hi a litle bit different story. But x1 - sample(c(rep(red,400),rep(green, 100), rep(black,300)),100) is maybe close. With data frame (if it is not big) DF color sample1 sample2 sample3 1 red 400 3002500 2 green 100 0 200 3 black 3001000 500 x - data.frame(matrix(NA,100,3)) for (i in 2:ncol(DF)) x[,i-1] - sample(rep(DF[,1], DF[,i]),100) if you want result in data frame or x-vector(list, 3) for (i in 2:ncol(DF)) x[[,i-1]] - sample(rep(DF[,1], DF[,i]),100) if you want it in list. Maybe somebody is clever enough to discard for loop but you said you have 80 columns which shall be no problem. HTH Petr On 11 Oct 2006 at 10:11, Brian Frappier wrote: Date sent:Wed, 11 Oct 2006 10:11:33 -0400 From: Brian Frappier [EMAIL PROTECTED] To: Petr Pikal [EMAIL PROTECTED] Subject: Fwd: [R] rarefy a matrix of counts -- Forwarded message -- From: Brian Frappier [EMAIL PROTECTED] Date: Oct 11, 2006 10:10 AM Subject: Re: [R] rarefy a matrix of counts To: r-help@stat.math.ethz.ch Hi Petr, Thanks for your response. I have data that looks like the following: sample 1 sample 2 sample 3 red candy400 300 2500 green candy1000 200 black candy 3001000500 I don't want to randomly select either the samples (columns) or the candy types (rows), which sample as you state would allow me. Instead, I want to randomly sample 100 candies from each sample and retain info on their associated type. I could make a list of all the candies in each sample: sample 1 red red red red green green black red black ... and then randomly sample those rows. Repeat for each sample. But, I am not sure how to do that without alot of loops, and am wondering if there is an easier way in R. Thanks! I should have laid this out in the first email...sorry. On 10/11/06, Petr Pikal [EMAIL PROTECTED] wrote: Hi I am not experienced in Matlab and from your explanation I do not understand what exactly do you want. It seems that you want randomly choose a sample of 100 rows from your martix, what can be achived by sample. DF-data.frame(rnorm(100), 1:100, 101:200, 201:300) DF[sample(1:100, 10),] If you want to do this several times, you need to save your
Re: [R] Simulate p-value in lme4
On Sun, 2006-10-08 at 07:34 -0400, Michael Kubovy wrote:
Dear r-helpers,
Spencer Graves and Manual Morales proposed the following methods to
simulate p-values in lme4:
preliminary
require(lme4)
require(MASS)
summary(glm(y ~ lbase*trt + lage + V4, family = poisson, data =
epil), cor = FALSE)
epil2 - epil[epil$period == 1, ]
epil2[period] - rep(0, 59); epil2[y] - epil2[base]
epil[time] - 1; epil2[time] - 4
epil2 - rbind(epil, epil2)
epil2$pred - unclass(epil2$trt) * (epil2$period 0)
epil2$subject - factor(epil2$subject)
epil3 - aggregate(epil2, list(epil2$subject, epil2$period 0),
function(x) if(is.numeric(x)) sum(x) else x[1])
simulation (SG)
o - with(epil3, order(subject, y))
epil3. - epil3[o,]
norep - with(epil3., subject[-1]!=subject[-dim(epil3)[1]])
subj1 - which(c(T, norep))
subj.pred - epil3.[subj1, c(subject, pred)]
subj. - as.character(subj.pred$subject)
pred. - subj.pred$pred
names(pred.) - subj.
iter - 10
chisq.sim - rep(NA, iter)
set.seed(1)
for(i in 1:iter){
s.i - sample(subj.)
# Randomize subject assignments to 'pred' groups
epil3.$pred - pred.[s.i][epil3.$subject]
fit1 - lmer(y ~ pred+(1 | subject),
family = poisson, data = epil3.)
fit0 - lmer(y ~ 1+(1 | subject),
family = poisson, data = epil3.)
chisq.sim[i] - anova(fit0, fit1)[2, Chisq]
}
simulation (MM)
iter - 10
chisq.sim - rep(NA, iter)
Zt - slot(model1,Zt) # see ?lmer-class
n.grps - dim(ranef(model1)[[1]])[1]
sd.ran.effs - sqrt(VarCorr(model1)[[1]][1])
X - slot(model1,X) # see ?lmer-class
fix.effs - matrix(rep(fixef(model1),dim(X)[1]), nrow=dim(X)[1],
byrow=T)
model.parms - X*fix.effs # This gives parameters for each case
# Generate predicted values
pred.vals - as.vector(apply(model.parms, 1, sum))
for(i in 1:iter) {
rand.new - as.vector(rnorm(grps,0, sd.ran.effs)) #grps should be
#n.grps?
Yes. Change grps to n.grps.
rand.vals - as.vector(rand.new%*%Zt) # Assign random effects
mu - pred.vals+rand.vals # Expected mean
resp - rpois(length(mu), exp(mu))
sim.data - data.frame(slot(model2,frame), resp) # Make data frame
sim.model1 - lmer(resp~1+(1|subject), data=sim.data,
family=poisson)
sim.model2 - lmer(resp~pred+(1|subject), data=sim.data,
family=poisson)
chisq.sim[i] - anova(sim.model1,sim.model2)$Chisq[[2]]
}
end
Unfortunately the latter fails (even if I replace grps with n.grps):
Error in slot(model2, frame) : object model2 not found
You need to specify the models fit0 and fit1 from SG's example as model1
and model2. E.g., model1 - fit0, etc.
In any event, I would be eager to hear more discussion of the pros
and cons of these approaches.
One practical problem with my approach (MM) is that the fitting
algorithms for lmer would often choke - in particular for those
randomly generated data sets that by chance included variance components
close to zero.
In any case, the MCMC approach advocated by Douglas Bates is *much*
faster. That's the approach I've been using since he suggested a
function for estimating P-values from MCMC samples for factors (or model
comparisons) with greater than 2 levels.
See:http://tolstoy.newcastle.edu.au/R/e2/help/06/09/1003.html
and the very long thread that accompanies it.
HTH,
Manuel
--
Manuel A. Morales
http://mutualism.williams.edu
signature.asc
Description: This is a digitally signed message part
__
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] Conservative ANOVA tables in lmer
On Wed, 2006-09-13 at 08:04 +1000, Andrew Robinson wrote: On Tue, September 12, 2006 7:34 am, Manuel Morales wrote: On Mon, 2006-09-11 at 11:43 -0500, Douglas Bates wrote: Having made that offer I think I will now withdraw it. Peter's example has convinced me that this is the wrong thing to do. I am encouraged by the fact that the results from mcmcsamp correspond closely to the correct theoretical results in the case that Peter described. I appreciate that some users will find it difficult to work with a MCMC sample (or to convince editors to accept results based on such a sample) but I think that these results indicate that it is better to go after the marginal distribution of the fixed effects estimates (which is what is being approximated by the MCMC sample - up to Bayesian/frequentist philosophical differences) than to use the conditional distribution and somehow try to adjust the reference distribution. Am I right that the MCMC sample can not be used, however, to evaluate the significance of parameter groups. For example, to assess the significance of a three-level factor? Are there better alternatives than simply adjusting the CI for the number of factor levels (1-alpha/levels). I wonder whether the likelihood ratio test would be suitable here? That seems to be supported. It just takes a little longer. require(lme4) data(sleepstudy) fm1 - lmer(Reaction ~ Days + (Days|Subject), sleepstudy) fm2 - lmer(Reaction ~ Days + I(Days^2) + (Days|Subject), sleepstudy) anova(fm1, fm2) So, a brief overview of the popular inferential needs and solutions would then be: 1) Test the statistical significance of one or more fixed or random effects - fit a model with and a model without the terms, and use the LRT. I believe that the LRT is anti-conservative for fixed effects, as described in Pinheiro and Bates companion book to NLME. 2) Obtain confidence intervals for one or more fixed or random effects - use mcmcsamp Did I miss anything important? - What else would people like to do? Cheers Andrew Andrew Robinson Senior Lecturer in Statistics Tel: +61-3-8344-9763 Department of Mathematics and StatisticsFax: +61-3-8344 4599 University of Melbourne, VIC 3010 Australia Email: [EMAIL PROTECTED]Website: http://www.ms.unimelb.edu.au __ 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@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] Conservative ANOVA tables in lmer
On Mon, 2006-09-11 at 11:43 -0500, Douglas Bates wrote: On 9/10/06, Andrew Robinson [EMAIL PROTECTED] wrote: On Thu, Sep 07, 2006 at 07:59:58AM -0500, Douglas Bates wrote: I would be happy to re-institute p-values for fixed effects in the summary and anova methods for lmer objects using a denominator degrees of freedom based on the trace of the hat matrix or the rank of Z:X if others will volunteer to respond to the these answers are obviously wrong because they don't agree with whatever and the idiot who wrote this software should be thrashed to within an inch of his life messages. I don't have the patience. This seems to be more than fair to me. I'll volunteer to help explain why the anova.lmer() output doesn't match SAS, etc. Is it worth putting a caveat in the output and the help files? Is it even worth writing a FAQ about this? Having made that offer I think I will now withdraw it. Peter's example has convinced me that this is the wrong thing to do. I am encouraged by the fact that the results from mcmcsamp correspond closely to the correct theoretical results in the case that Peter described. I appreciate that some users will find it difficult to work with a MCMC sample (or to convince editors to accept results based on such a sample) but I think that these results indicate that it is better to go after the marginal distribution of the fixed effects estimates (which is what is being approximated by the MCMC sample - up to Bayesian/frequentist philosophical differences) than to use the conditional distribution and somehow try to adjust the reference distribution. Am I right that the MCMC sample can not be used, however, to evaluate the significance of parameter groups. For example, to assess the significance of a three-level factor? Are there better alternatives than simply adjusting the CI for the number of factor levels (1-alpha/levels). Thanks! Manuel = Manuel A. Morales Asst. Prof., Biology Williams College http://mutualism.williams.edu __ 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] ess-remote question
Dear list, Is there any way to load a local data file when connected to a remote machine via ESS? Thanks! Manuel __ 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] Simulate p-value in lme4
Spencer,
Thanks for the reply. I concluded the same wrt between group variation
soon after posting. However, the approach I ended up with was fully
parametric as opposed to the resampling approach that you use in your
reply. Interestingly, the two approaches yield different P-values, I
think because your approach retains overdispersion in the data (?). In
any case, my parametric stab at this is below.
iter - 10
chisq.sim - rep(NA, iter)
Zt - slot(model1,Zt) # see ?lmer-class
n.grps - dim(ranef(model1)[[1]])[1]
sd.ran.effs - sqrt(VarCorr(model1)[[1]][1])
X - slot(model1,X) # see ?lmer-class
fix.effs - matrix(rep(fixef(model1),dim(X)[1]), nrow=dim(X)[1],
byrow=T)
model.parms - X*fix.effs # This gives parameters for each case
# Generate predicted values
pred.vals - as.vector(apply(model.parms, 1, sum))
for(i in 1:iter) {
rand.new - as.vector(rnorm(grps,0, sd.ran.effs))
rand.vals - as.vector(rand.new%*%Zt) # Assign random effects
mu - pred.vals+rand.vals # Expected mean
resp - rpois(length(mu), exp(mu))
sim.data - data.frame(slot(model2,frame), resp) # Make data frame
sim.model1 - lmer(resp~1+(1|subject), data=sim.data,
family=poisson)
sim.model2 - lmer(resp~pred+(1|subject), data=sim.data,
family=poisson)
chisq.sim[i] - anova(sim.model1,sim.model2)$Chisq[[2]]
}
Manuel
On Sun, 2006-08-20 at 11:22 -0700, Spencer Graves wrote:
You've raised a very interesting question about testing a
fixed-effect factor with more than 2 levels using Monte Carlo. Like
you, I don't know how to use 'mcmcsamp' to refine the naive
approximation. If we are lucky, someone else might comment on this for us.
Beyond this, you are to be commended for providing such a simple,
self-contained example for such a sophisticated question. I think you
simulation misses one important point: It assumes the between-subject
variance is zero. To overcome this, I think I might try either the
bootstrap or permutation distribution scrambling the assignment of
subjects to treatment groups but otherwise keeping the pairs of
observations together.
To this end, consider the following:
# Build a table to translate subject into 'pred'
o - with(epil3, order(subject, y))
epil3. - epil3[o,]
norep - with(epil3., subject[-1]!=subject[-dim(epil3)[1]])
subj1 - which(c(T, norep))
subj.pred - epil3.[subj1, c(subject, pred)]
subj. - as.character(subj.pred$subject)
pred. - subj.pred$pred
names(pred.) - subj.
iter - 10
chisq.sim - rep(NA, iter)
set.seed(1)
for(i in 1:iter){
## Parameteric version
s.i - sample(subj.)
# Randomize subject assignments to 'pred' groups
epil3.$pred - pred.[s.i][epil3.$subject]
fit1 - lmer(y ~ pred+(1 | subject),
family = poisson, data = epil3.)
fit0 - lmer(y ~ 1+(1 | subject),
family = poisson, data = epil3.)
chisq.sim[i] - anova(fit0, fit1)[2, Chisq]
}
Hope this helps.
Spencer Graves
Manuel Morales wrote:
Dear list,
This is more of a stats question than an R question per se. First, I
realize there has been a lot of discussion about the problems with
estimating P-values from F-ratios for mixed-effects models in lme4.
Using mcmcsamp() seems like a great alternative for evaluating the
significance of individual coefficients, but not for groups of
coefficients as might occur in an experimental design with 3 treatment
levels. I'm wondering if the simulation approach I use below to estimate
the P-value for a 3-level factor is appropriate, or if there are any
suggestions on how else to approach this problem. The model and data in
the example are from section 10.4 of MASS.
Thanks!
Manuel
# Load req. package (see functions to generate data at end of script)
library(lme4)
library(MASS)
# Full and reduced models - pred is a factor with 3 levels
result.full - lmer(y~pred+(1|subject), data=epil3, family=poisson)
result.base - lmer(y~1+(1|subject), data=epil3, family=poisson)
# Naive P-value from LR for significance of pred factor
anova(result.base,result.full)$Pr(Chisq)[[2]] # P-value
(test.stat - anova(result.base,result.full)$Chisq[[2]]) # Chisq-stat
snip Wrong approach here/snip
# Script to generate data, from section 10.4 of MASS
epil2 - epil[epil$period == 1, ]
epil2[period] - rep(0, 59); epil2[y] - epil2[base]
epil[time] - 1; epil2[time] - 4
epil2 - rbind(epil, epil2)
epil2$pred - unclass(epil2$trt) * (epil2$period 0)
epil2$subject - factor(epil2$subject)
epil3 - aggregate(epil2, list(epil2$subject, epil2$period 0),
function(x) if(is.numeric(x)) sum(x) else x[1])
epil3$pred - factor(epil3$pred, labels = c(base, placebo, drug))
__
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
[R] Simulate p-value in lme4
Dear list,
This is more of a stats question than an R question per se. First, I
realize there has been a lot of discussion about the problems with
estimating P-values from F-ratios for mixed-effects models in lme4.
Using mcmcsamp() seems like a great alternative for evaluating the
significance of individual coefficients, but not for groups of
coefficients as might occur in an experimental design with 3 treatment
levels. I'm wondering if the simulation approach I use below to estimate
the P-value for a 3-level factor is appropriate, or if there are any
suggestions on how else to approach this problem. The model and data in
the example are from section 10.4 of MASS.
Thanks!
Manuel
# Load req. package (see functions to generate data at end of script)
library(lme4)
library(MASS)
# Full and reduced models - pred is a factor with 3 levels
result.full - lmer(y~pred+(1|subject), data=epil3, family=poisson)
result.base - lmer(y~1+(1|subject), data=epil3, family=poisson)
# Naive P-value from LR for significance of pred factor
anova(result.base,result.full)$Pr(Chisq)[[2]] # P-value
(test.stat - anova(result.base,result.full)$Chisq[[2]]) # Chisq-stat
# P-value from simulation. Note that in the simulation, I use the
# estimated random effects for each subject rather than generating a new
# distribution of means. I'm not sure if this is appropriate or not ...
intercept - fixef(result.base)
rand.effs - ranef(result.base)[[1]]
mu - exp(rep(intercept+rand.effs[[1]],2))
p.value - function(iter, stat) {
chi.stat - vector()
for(i in 1:iter) {
resp - rpois(length(mu), mu) # simulate values
sim.data - data.frame(y=resp,subject=epil3$subject,pred=epil3$pred)
result.f - lmer(y~pred+(1|subject), data=sim.data,
family=poisson)
result.b - lmer(y~1+(1|subject), data=sim.data, family=poisson)
chi.stat[i] - anova(result.b,result.f)$Chisq[[2]]
}
val - sum(unlist(lapply(chi.stat, function(x) if(xstat) 1 else
0)))/iter
hist(chi.stat)
return(val)
}
p.value(10,test.stat) # Increase to =1000 to get a reasonable P-value!
# Script to generate data, from section 10.4 of MASS
epil2 - epil[epil$period == 1, ]
epil2[period] - rep(0, 59); epil2[y] - epil2[base]
epil[time] - 1; epil2[time] - 4
epil2 - rbind(epil, epil2)
epil2$pred - unclass(epil2$trt) * (epil2$period 0)
epil2$subject - factor(epil2$subject)
epil3 - aggregate(epil2, list(epil2$subject, epil2$period 0),
function(x) if(is.numeric(x)) sum(x) else x[1])
epil3$pred - factor(epil3$pred, labels = c(base, placebo, drug))
__
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] Using data=x or subset=y in user-defined functions
On Wed, 2006-06-07 at 21:36 +0100, Prof Brian Ripley wrote:
I suggest you investigate with().
Thanks! Just to be sure, the following will do what I want?
eg.function - function(x, data=NULL, subset=NULL, ...) {
with(if(is.null(subset)) data else subset(data,subset), {
...
})
}
On Wed, 7 Jun 2006, Manuel Morales wrote:
Dear list members,
In some of my functions, I attach the data internally to allow subset
commands or to specify a data frame. This works well except for cases
where there is a masking conflict (which returns a warning). I see
some alternative listed in ?attach, but I'm not sure which of them do
what I'd like. Any suggestions?
Below is how I've been setting up my functions:
# Set up environment
on.exit(detach(data))
attach(data)
if(!is.null(subset)) {
data-subset(data,subset)
detach(data)
attach(data)
}
subset = NULL
# Function body here
output - x
return(output)
}
Thanks!
Manuel
__
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] Using data=x or subset=y in user-defined functions
On Thu, 2006-06-08 at 11:54 -0400, Manuel Morales wrote:
On Wed, 2006-06-07 at 21:36 +0100, Prof Brian Ripley wrote:
I suggest you investigate with().
Thanks for the suggestion! Unfortunately, it's not clear to me how to
call with() from a function. The example below fails with the error
message: 'Error in mean(x) : object a not found'
data.test - data.frame(a=c(1:10),b=rnorm(10))
eg.function - function(x, data)
with(data, return(mean(x)))
eg.function(a,data.test)
Manuel
Thanks! Just to be sure, the following will do what I want?
eg.function - function(x, data=NULL, subset=NULL, ...) {
with(if(is.null(subset)) data else subset(data,subset), {
...
})
}
On Wed, 7 Jun 2006, Manuel Morales wrote:
Dear list members,
In some of my functions, I attach the data internally to allow subset
commands or to specify a data frame. This works well except for cases
where there is a masking conflict (which returns a warning). I see
some alternative listed in ?attach, but I'm not sure which of them do
what I'd like. Any suggestions?
Below is how I've been setting up my functions:
# Set up environment
on.exit(detach(data))
attach(data)
if(!is.null(subset)) {
data-subset(data,subset)
detach(data)
attach(data)
}
subset = NULL
# Function body here
output - x
return(output)
}
Thanks!
Manuel
__
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] error bars in lattice xyplot *with groups*
Hi Mike,
If you're not committed to using a panel function, another option is to
use the function lineplot.CI, available in the package sciplot at
http://mutualism.williams.edu/sciplot
E.g.
# Define and generate variables in long format
range - vector()
voice - vector()
string - strsplit(as.character(singer$voice.part), )
for(i in 1:dim(singer)[1]) {
range[i] - string[[i]][1]
voice[i] - string[[i]][2]
}
# Define function for CI
conf.int - function(x) {
st - boxplot.stats(x)
c((st$conf[2]-st$conf[1])/2)
}
# Plot
library(sciplot)
lineplot.CI(response=height, x.factor=voice, trace.factor=range,
data=singer, fun=median, ci.fun=conf.int)
lineplot.CI(response=height, x.factor=voice.part, data=singer,
fun=median, ci.fun=conf.int)
Manuel
On Tue, 2006-06-06 at 00:20 -0300, Mike Lawrence wrote:
Hi all,
I'm trying to plot error bars in a lattice plot generated with xyplot.
Deepayan
Sarkar has provided a very useful solution for simple circumstances
(https://stat.ethz.ch/pipermail/r-help/2005-October/081571.html), yet I am
having trouble getting it to work when the groups setting is enabled in
xyplot (i.e. multiple lines). To illustrate this, consider the singer data
generated by the above linked solution previously submitted:
#
library(lattice)
singer.split -
with(singer,
split(height, voice.part))
singer.ucl -
sapply(singer.split,
function(x) {
st - boxplot.stats(x)
c(st$stats[3], st$conf)
})
singer.ucl - as.data.frame(t(singer.ucl))
names(singer.ucl) - c(median, lower, upper)
singer.ucl$voice.part -
factor(rownames(singer.ucl),
levels = rownames(singer.ucl))
#now let's split up the voice.part factor into two factors,
singer.ucl$voice=factor(rep(c(1,2),4))
singer.ucl$range=factor(rep(c(Bass,Tenor,Alto,Soprano),each=2))
#here's Deepayan's previous solution, slightly modified to depict
# the dependent variable (median) and the error bars on the y-axis
# and the independent variable (voice.part) on the x-axis
prepanel.ci - function(x, y, ly, uy, subscripts, ...)
{
x - as.numeric(x)
ly - as.numeric(ly[subscripts])
uy - as.numeric(uy[subscripts])
list(ylim = range(y, uy, ly, finite = TRUE))
}
panel.ci - function(x, y, ly, uy, subscripts, pch = 16, ...)
{
x - as.numeric(x)
y - as.numeric(y)
ly - as.numeric(ly[subscripts])
uy - as.numeric(uy[subscripts])
panel.arrows(x, ly, x, uy, col = black,
length = 0.25, unit = native,
angle = 90, code = 3)
panel.xyplot(x, y, pch = pch, ...)
}
#this graph works
xyplot(median ~ voice.part,
data=singer.ucl,
ly = singer.ucl$lower,
uy = singer.ucl$upper,
prepanel = prepanel.ci,
panel = panel.ci,
type=b
)
#this one does not (it will plot, but will not seperate the groups)
xyplot(median ~ voice,
groups=range,
data=singer.ucl,
ly = singer.ucl$lower,
uy = singer.ucl$upper,
prepanel = prepanel.ci,
panel = panel.ci,
type=b
)
Any suggestions?
__
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 data=x or subset=y in user-defined functions
Dear list members,
In some of my functions, I attach the data internally to allow subset
commands or to specify a data frame. This works well except for cases
where there is a masking conflict (which returns a warning). I see
some alternative listed in ?attach, but I'm not sure which of them do
what I'd like. Any suggestions?
Below is how I've been setting up my functions:
eg.function - function(x, data=NULL, subset=NULL, ...) {
# Set up environment
on.exit(detach(data))
attach(data)
if(!is.null(subset)) {
data-subset(data,subset)
detach(data)
attach(data)
}
subset = NULL
# Function body here
output - x
return(output)
}
Thanks!
Manuel
__
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] MART(tm) vs. gbm
Hello, I tried two different implementations of the stochastic gradient boosting (Friedman 2002) : the MART(tm) with R tool (http://www-stat.stanford.edu/~jhf/R-MART.html) and the gbm R package. To me, both seemed fairly comparable, except maybe regarding the different loss criterion proposed and the fact that the gbm tool is sightly more convenient to use. However, it seemed that the MART with R tool systematically outperforms the gbm tool in terms of goodness of fit (whatever the way of choosing the best iteration for the gbm package). I tried to find out if there were specific options that could have explained it but nothing came out. See below for an example of how I compare both implementations. Did any one had the same experience, and can anyone give me hints about such performance differences or tell me if I am missing something obvious? Thank you in advance,Manuel Here are the arguments and options I used for comparison purposes, working on a 1600 records * 15 variables dataset : # the MART with R tool lx - mart( as.matrix(x), y, c(1,1,1,1,1,1,1,1,1,1,1,1,1,2,2) niter=1000, tree.size=6, learn.rate=0.01, loss.cri=2 # gaussian ) # for gbm gbm1 - gbm(y ~ v1 + v2 + v3 + v4 + v5+ v6+ v7+ v8 + v9 + v10 + v11 + v12 + v13 + v14 + v15, data=data, var.monotone=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), distribution=gaussian, n.trees=1000, shrinkage=0.01, interaction.depth=6, bag.fraction = 0.5, train.fraction = 0.5, n.minobsinnode = 10, cv.folds = 1, keep.data=TRUE) # I then do predictions on the same dataset, and further perform goodness of fit comparisons #... __ 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] Can't there be a cd command?
Hi, I am afraid this discussion is going into the wrong direction. I do agree with some comments of Joerg van den Hoff, however, I cannot adhere to some things said in the last paragraphs. I was just expressing a personal observation and I was expecting many people to have disagree and perhaps very few to agree. However, it was just an opinion. Unfortunately, during the next months I am going to have little free time, so I cannot: (1) gather examples of usability issues in R, (2) support them by user cases from R-help, (3) propose solutions and put forth reasons and argue for the changes one by one.[*] In a software libre project, complaining without contributing is futile, pointless and even insulting to the people that do contribute. That is the reason why I did not send any further comments on this. Of course, I would be very happy if someone (or some group) decided to start to work on this. However, as said before, I don't consider I have the right to tell anybody that they must do it. Thus, a discussion on whether R is hard/easy or R has to be hard/easy or R is harder/easier than program/language X or R is like deciphering hieroglyphics / R is like piloting an Apache helicopter is pointless in my opinion. So please, don't quote me or mention me in such kind of discussion. Please, don't even reply to such messages: there is already enough traffic in R-help. Now, Joerg van de Hoff points out particular cases (like subsetting/indexing issues) where new users seem to always have problems. That is a better approach: focusing on actual cases. Still, more work needs to be done to identify where the problems are and how to solve them. That would imply to examine the reaction of users (from R-help or your own students), since your (my) personal experience is almost useless once you know the answer (magic syntax or workaround). Therefore, the value of I think this is hard/easy because it is hard/easy for me is close to zero. Such discussion may end up on people calling names and I don't want to be involved on that. I think working mainly off-list on particular proposals (perhaps sharing information in the R-Wiki) would be the ideal way to work on this. I am sorry I cannot invest more time on this at the present moment. Regards, Manuel López-Ibáñez. PS: by the way, I do use Perl, Emacs and LaTeX (almost everyday) and I think they are great, yet they could be improved in terms of usability. [*] A clear candidate is the cd command. cd means change directory in Windows, Unix and dozens of applications. It can be argued that ls doesn't mean list files. For me, the fact that ls has another meaning is just unfortunate and I understand that it may be problematic to change it. However, cd doesn't mean anything in R yet. Actually, there is a dir command. Could you guess what dir() does? If still you are not convinced, let it be, I cannot discuss this further (perhaps after summer). Joerg van den Hoff wrote: Duncan Murdoch wrote: On 5/16/2006 5:46 AM, Joerg van den Hoff wrote: Manuel López-Ibáñez wrote: Jan T. Kim wrote: That's an idea I like very much too -- much better than the currently popular idea of protecting users from the unfriendliness of programming, anyway... It is just my opinion that the amount of mail in R-help speaks volumes about the current friendliness [1], or lack thereof, of R. Perhaps I am just the only one who thinks this way... [1] http://en.wikipedia.org/wiki/Usability I think you are 100% right: the r-help list says it all. needless to say, R is a great achievment without any doubt, but claiming that it's easy to use (beyond the most basic arithmetics) is really wishful thinking. This is sloppy thinking. The volume of mail here shows that there are a lot of questions, perhaps because there are a lot of users. well, as far as my english goes, 'sloppy' is a strong word (and apart from mathematicians physicists (my background) probably are the people who are most allergic to being accused of it :-)) and it's an overhasty conclusion on your side, I'd say. I want to make clear beforehand, that I do _not_ think this a very important discussion, but rather an informal exchange of opinions, so maybe this takes us all a bit to far, but anyway: for one, I myself (and I think manuel, too) was not talking of the shear volume of mails (this obviously would have to be 'calibrated' against the total number of R users and the resulting quantity had to be compared to other help-lists). at least my impression is, that there are a number of reoccuring difficulties in the mail, which are rather specific to R's design (whether this situation could or should be altered: that would be a different topic). certainly, among these are the subsetting/indexing issues, certainly lazy evaluation, certainly anything related to environments, namespaces, computing on the language (substitute, eval, ...). You're also
Re: [R] Can't there be a cd command?
Jan T. Kim wrote: That's an idea I like very much too -- much better than the currently popular idea of protecting users from the unfriendliness of programming, anyway... It is just my opinion that the amount of mail in R-help speaks volumes about the current friendliness [1], or lack thereof, of R. Perhaps I am just the only one who thinks this way... [1] http://en.wikipedia.org/wiki/Usability __ LLama Gratis a cualquier PC del Mundo. Llamadas a fijos y móviles desde 1 céntimo por minuto. http://es.voice.yahoo.com __ 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] Making R talk to Win/OpenBUGS in Linux (again)
As far as I can figure out, the problem with running LinBUGS on FC5 is
that support for linuxthreads was removed after being deprecated in FC4.
http://fedora.redhat.com/docs/release-notes/fc5/#id2887615
As a result, the LD_ASSUME_KERNEL workaround, which presumably forced
the use of linuxthreads (see
http://www.math.helsinki.fi/openbugs/phpBB2/viewtopic.php?t=4), doesn't
work anymore.
Ideally, bugs.so / brugs.so would be compiled with support for NPTL
instead of linuxthreads, but I don't think it's possible to do this as
an end user. Unless changes are made by the OpenBUGS developers,
OpenBUGS will not run natively in FC5.
Manuel
On Mon, 2006-05-01 at 13:38 -0500, Paul Johnson wrote:
Dear Jun:
How about telling us which version of Linux you use and how you make
linbugs run? As far as I can tell, the OpenBUGS people have
intentionally removed linbugs from their version 2.2. That leaves us
with various scripts that people have posted tried, none work for
me. Fedora Core 5, I cannot get any of the linbugs scripts that are
floating around to work.
I tried re-building the program thus:
# gcc -o cbugs CBugs.c -ldl
That creates an executable cbugs, but there is no joy for me.
$ ./cbugs
failed to install signal
32
failed to install signal
33
* BlackBox
* illegal memory read [ad = ]
- HostFiles.NewFileRef (pc=0D9B, fp=BFC73BE4)
- HostFiles.OpenFile (pc=0E0E, fp=BFC73BFC)
- HostFiles.Directory.Old (pc=29B2, fp=BFC73EA4)
- StdLoader.ThisObjFile (pc=031D, fp=BFC744C0)
- StdLoader.ReadHeader (pc=075F, fp=BFC746E0)
- StdLoader.LoadMod (pc=107B, fp=BFC747F8)
- StdLoader.LoadMod (pc=11E0, fp=BFC74910)
- StdLoader.LoadMod (pc=11E0, fp=BFC74A28)
- StdLoader.Hook.ThisMod (pc=1355, fp=BFC74A3C)
- Kernel.ThisMod (pc=1224, fp=BFC74B58)
- Init.Init (pc=004B, fp=BFC74B70)
- Init.$$ (pc=000A, fp=BFC74B80)
# export LD_ASSUME_KERNEL=2.4.1
# ./cbugs
./cbugs: error while loading shared libraries: libdl.so.2: cannot open
shared object file: No such file or directory
$ ldd ./cbugs
linux-gate.so.1 = (0x00a2b000)
libdl.so.2 = /lib/libdl.so.2 (0x00ba9000)
libc.so.6 = /lib/libc.so.6 (0x00a4d000)
/lib/ld-linux.so.2 (0x00a2c000)
And after doing that LD_ASSUME_KERNEL thingie, nothing in the shell
works anymore...
$ which cbugs
/usr/bin/which: error while loading shared libraries: libc.so.6:
cannot open shared object file: No such file or directory
On 5/1/06, jun yan [EMAIL PROTECTED] wrote:
I have used linbugs with the rbugs package for a recent work. It might be
worthwhile trying.
Jun
On 5/1/06, Uwe Ligges [EMAIL PROTECTED]
wrote:
Paul Johnson wrote:
Thank you very much. With the insertion of WINEPATH declaration, then
the following example program does run. And really fast, too!
library(R2WinBUGS)
WINEPATH - /usr/bin/winepath
Will be fixed in the package real soon now.
Gregor, many thanks for the patches.
Best,
Uwe Ligges
# An example
model file is given in:
model.file - system.file(package = R2WinBUGS, model, schools.txt)
# Let's take a look:
file.show(model.file)
# Some example data (see ?schools for details):
data(schools)
schools
J - nrow(schools)
y - schools$estimate
sigma.y - schools$sd
data - list (J, y, sigma.y)
inits - function(){
list(theta = rnorm(J, 0, 100), mu.theta = rnorm(1, 0, 100),
sigma.theta = runif(1, 0, 100))
}
parameters - c(theta, mu.theta, sigma.theta)
schools.sim - bugs(data, inits, parameters, model.file, n.chains = 3,
n.iter = 1000,bugs.directory = /usr/local/share/WinBUGS14/,
working.directory = /tmp, clearWD = FALSE, useWINE=T, newWINE=T,
WINE=/usr/bin/wine,WINEPATH=WINEPATH)
On 4/30/06, Gregor Gorjanc [EMAIL PROTECTED] wrote:
Hello Paul,
thank you very much for this report. You caught a bug in R2WinBUGS that
was introduced by me. I added support for winepath in 1.1-1 version.
Since I switch between Windows and Linux I always set WINEPATH and then
use it in bugs(). That's why I forgot to add it in further calls in
bugs(). How silly :( This actually means that nobody else beside you
tried to run R2WinBUGS under Linux. To bad. Please report summary to
BUGS list as you have also asked there for help and I did not had time
to answer you. I have been successfully running R2WinBUGS under Linux
for quite some time now.
Now you have the following options:
A Set WINEPATH before you call bugs()
i.e.
WINEPATH - /usr/bin/winepath
bugs(..., WINEPATH=WINEPATH)
This is the fastest approach for you.
B Apply the following diffs, build and install R2WinBUGS by yourself
I am also
Re: [R] nls and factor
Thanks, it was actually p.249, at least in my MASS3. but that solved my doubt. I've have another doubt, can this factor interact with one of the parameters in the model? My problem is basically a Michaelis Menten term, where this factor determines a different Km. The rest of the parameters in the model are the same. But I don't know how to write the nls formula, or if it is possible. This is a toy example, B and A are the two factors, conc is the concentration and t is the temperature: ## Generate independent variables Bconc-runif(30,0.1,10) Aconc-runif(30,0.1,10) At-runif(30,1,30) Bt-runif(30,1,30) ## These are the parameters I want to calculate from ## my real data BKm-1 AKm-0 EBoth--0.41 # These are my simulated dependent variables yB-100*exp(EBoth*Bt)*Bconc/(BKm+Bconc)+rnorm(30,0,1) yA-75*exp(EBoth*At)*Aconc/(AKm+Aconc)+rnorm(30,0,1) #The separate models BModel-nls(Response~lev*exp(Ev*t)*conc/(Km+conc),data=list(Response=yB,t=Bt,conc=Bconc),start=list(lev=90,Ev=-0.5,Km=0.8),trace=TRUE) AModel-nls(Response~lev*exp(Ev*t)*conc/(Km+conc),data=list(Response=yA,t=At,conc=Aconc),start=list(lev=90,Ev=-0.5,Km=0.8),trace=TRUE) ## I want to obtain a combined model of the form: ## Y=Intercept[1:2]*exp(Eboth*t)*conc/(Km[1:2]+conc) ## where I have a common E but two intercepts and two ## Kms (one of them should in fact be zero) yBoth-c(yB,yA) concBoth-c(Bconc,Aconc) tBoth-c(At,Bt) AorB-as.factor(c(rep(0,length(yA)),rep(1,length(yB ## Amongst other things I've tried FullModel-nls(Response~lev[AorB]*exp(Ev*t)*conc/(Km[AorB]+conc),data=list(Response=yBoth,t=tBoth,conc=concBoth),start=list(lev=c(90,70),Ev=-0.5,Km=c(0.8,0)),trace=TRUE) ## but i get to a singular gradient Any other pointers, thanks Manuel --- Prof Brian Ripley [EMAIL PROTECTED] escribió: On Thu, 20 Apr 2006, Manuel Gutierrez wrote: Is it possible to include a factor in an nls formula? Yes. What do you intend by it? If you mean what it would mean for a lm formula, you need A[a] and starting values for A. There's an example on p.219 of MASS4. I've searched the help pages without any luck so I guess it is not feasible. I've given it a few attempts without luck getting the message: + not meaningful for factors in: Ops.factor(independ^EE, a) This is a toy example, my realworld case is much more complicated (and can not be solved linearizing an using lm) a-as.factor(c(rep(1,50),rep(0,50))) independ-rnorm(100) respo-rep(NA,100) respo[a==1]-(independ[a==1]^2.3)+2 respo[a==0]-(independ[a==0]^2.1)+3 nls(respo~independ^EE+a,start=list(EE=1.8),trace=TRUE) Any pointers welcomed Many Thanks, Manu __ 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 -- 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
Re: [R] www.r-project.org
Barry Rowlingson wrote:
The frame-based nature of the CRAN pages is slightly problematic,
since you click on a menu item and the URL doesn't change. Hence there's
no way to send someone a URL that gives them the same view as you'd get
if you go to the home page and then click on 'screenshots', for example.
Sure you can send them to:
http://www.r-project.org/screenshots/screenshots.html
but then they dont see the menu.
Moreover, the menu gets lost also when an element of the menu is opened
in a new window (or tab).
It probably wouldn't take long to bash out a serviceable replacement
using something like HTML::Mason but then you'd have to find a hosting
provider that supported it (or PHP or IYFTLH[1]). I dont think it
warrants a full-on CMS given the size of www.r-project.org (not
including CRAN stuff). I'd just hack up some m4 scripts and 'include'
the menu into a flat file.
I have a perl script that reads a number of HTML files looking for
include directories which instruct it to put at that point the content
of another file.
For example, file index.html may contain the line:
!--#include file=menu.htm--
The script replaces that line with the content of menu.htm.
Thus, the files from the site are processed before being uploaded to the
web site in order to create the final html pages. There is no need for
Apache directives, no scripts, no PHP, no CMS. And more important, no
frames!
Here is the script. You may consider it GPL ;)
#!/usr/bin/perl -w
@files = `ls *.html`;
$outdir = ../;
`mkdir -p $outdir`; #or die Cannot create $outdir;
FILE: foreach $file (@files) {
$file =~ s/\n//g;
open(INPUT, $file) or (warn(*** Cannot read $file\n) and next
FILE);
@buffer = INPUT;
close INPUT;
LINE: foreach $line (@buffer) {
if($line =~ /!--#include file=([\w\.]+)--/) {
open(INPUT, $1) or
(warn(*** Cannot read $1 included in $file\n) and next LINE);
$temp = join(, INPUT);
close INPUT;
$line =~ s[!--#include file=([\w\.]+)--]
[!-- $1 --\n$temp];
print \t$1 included in $file\n;
}
}
open(OUTPUT, $outdir/$file) or
(warn(*** Cannot write $outdir/$file\n) and next FILE);
print OUTPUT @buffer;
close OUTPUT;
print $file processed\n;
}
__
LLama Gratis a cualquier PC del Mundo.
Llamadas a fijos y móviles desde 1 céntimo por minuto.
http://es.voice.yahoo.com
__
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] nls and factor
Is it possible to include a factor in an nls formula? I've searched the help pages without any luck so I guess it is not feasible. I've given it a few attempts without luck getting the message: + not meaningful for factors in: Ops.factor(independ^EE, a) This is a toy example, my realworld case is much more complicated (and can not be solved linearizing an using lm) a-as.factor(c(rep(1,50),rep(0,50))) independ-rnorm(100) respo-rep(NA,100) respo[a==1]-(independ[a==1]^2.3)+2 respo[a==0]-(independ[a==0]^2.1)+3 nls(respo~independ^EE+a,start=list(EE=1.8),trace=TRUE) Any pointers welcomed Many Thanks, Manu __ 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
Thanks Andrew. I am now trying but without much success. I don't now how to give start values for the factor?. Could you give me an example solution with my toy example? a-as.factor(c(rep(1,50),rep(0,50))) independ-1:100 respo-rep(NA,100) respo[a==1]-(independ[a==1]^2.3)+2 respo[a==0]-(independ[a==0]^2.1)+3 library(nlme) gnls(respo~independ^b+a,start=list(b=1.8)) Many thanks. Manu --- Andrew Robinson [EMAIL PROTECTED] escribió: Manuel, I don't think that it works very easily. Instead, try gnls() in the nlme package. Cheers Andrew On Thu, Apr 20, 2006 at 11:18:02AM +0200, Manuel Gutierrez wrote: Is it possible to include a factor in an nls formula? I've searched the help pages without any luck so I guess it is not feasible. I've given it a few attempts without luck getting the message: + not meaningful for factors in: Ops.factor(independ^EE, a) This is a toy example, my realworld case is much more complicated (and can not be solved linearizing an using lm) a-as.factor(c(rep(1,50),rep(0,50))) independ-rnorm(100) respo-rep(NA,100) respo[a==1]-(independ[a==1]^2.3)+2 respo[a==0]-(independ[a==0]^2.1)+3 nls(respo~independ^EE+a,start=list(EE=1.8),trace=TRUE) Any pointers welcomed Many Thanks, Manu __ 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 -- Andrew Robinson Department of Mathematics and Statistics Tel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 Email: [EMAIL PROTECTED] http://www.ms.unimelb.edu.au __ 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] what happen?
zhang jian wrote: * n=PIKO[status=snag,] Error in [.data.frame(PIKO, status = snag, ) : unused argument(s) (status ...) = is the assignment operator, you should use == as in: n=PIKO[status==snag,] __ LLama Gratis a cualquier PC del Mundo. Llamadas a fijos y móviles desde 1 céntimo por minuto. http://es.voice.yahoo.com __ 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] permutation of rows of a matrix
Dear John, I understand what you mean. However, when someone is learning R for the first time or have little experience, such examples help to understand the connection of different parts of the language. Moreover, things that make sense once you know them, can be difficult to relate in the first place. For example, it would be interesting to know how many new R users don't know that there is a manual page for [. I hope you can understand my point of view (you may disagree, though.) Regards, Manuel. John Fox wrote: Dear Manuel, Although ?sample doesn't specifically describe permuting the rows of a matrix, it does say that sample(x) generates a random permutation of the elements of x (or 1:x). Indexing the rows of the matrix by a permutation of 1:x (where x is the number of rows) doesn't seem to be much of a leap. Regards, John John Fox Department of Sociology McMaster University Hamilton, Ontario Canada L8S 4M4 905-525-9140x23604 http://socserv.mcmaster.ca/jfox -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Manuel López-Ibáñez Sent: Saturday, April 15, 2006 9:44 AM To: r-help@stat.math.ethz.ch Subject: Re: [R] permutation of rows of a matrix help(sample) does not say anything about randomly permuting the rows of a matrix M by using M[sample(m,m),]. Perhaps it could be added as an example of use. John Fox wrote: Dear Jose, M[sample(m, m),] will randomly permute the rows of M. [You probably could have figured this out via help.search(permutation), which would have led you to sample().] Regards, John John Fox Department of Sociology McMaster University Hamilton, Ontario Canada L8S 4M4 905-525-9140x23604 http://socserv.mcmaster.ca/jfox -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of javargas Sent: Saturday, April 15, 2006 7:53 AM To: r-help@stat.math.ethz.ch Subject: [R] permutation of rows of a matrix How can I generate a random permutation between rows of a matrix M (of m rows and n columns)? Thanks for your help, Jose __ 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 __ LLama Gratis a cualquier PC del Mundo. Llamadas a fijos y móviles desde 1 céntimo por minuto. http://es.voice.yahoo.com __ 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 __ LLama Gratis a cualquier PC del Mundo. Llamadas a fijos y móviles desde 1 céntimo por minuto. http://es.voice.yahoo.com __ 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] permutation of rows of a matrix
help(sample) does not say anything about randomly permuting the rows of a matrix M by using M[sample(m,m),]. Perhaps it could be added as an example of use. John Fox wrote: Dear Jose, M[sample(m, m),] will randomly permute the rows of M. [You probably could have figured this out via help.search(permutation), which would have led you to sample().] Regards, John John Fox Department of Sociology McMaster University Hamilton, Ontario Canada L8S 4M4 905-525-9140x23604 http://socserv.mcmaster.ca/jfox -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of javargas Sent: Saturday, April 15, 2006 7:53 AM To: r-help@stat.math.ethz.ch Subject: [R] permutation of rows of a matrix How can I generate a random permutation between rows of a matrix M (of m rows and n columns)? Thanks for your help, Jose __ 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 __ LLama Gratis a cualquier PC del Mundo. Llamadas a fijos y móviles desde 1 céntimo por minuto. http://es.voice.yahoo.com __ 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] Multiple error bar plots
Hi Young-Jin, On Fri, 2006-03-24 at 10:42 -0500, Young-Jin Lee wrote: Dear R-lister I have a question about how to create multiple error bar plots. I found that I can use an errbar function from Hmisc package to create one error bar plot in which there are multiple data points (x, y) with the error bars. Thus, I know that I can get one error bar plot which consists of many data points. However, I did not find a way to put multiple error bar plots into one. I wrote my own package to deal with exactly this scenario. You can download it from: http://mutualism.williams.edu/sciplot/ The relevant function is bargraph.CI. BTW, I did this mostly to learn about writing packages in R, but it may be helpful to you. In other words, my data looks like this: X1 (group index), Y11 (measurement), error of Y11 (error in measurement) X2, Y21, error of Y21 ... Xn, Yn1, error of Yn1 // End of the first set of measurements X1, Y12, error of Y12, X2, Y22, error of Y22, ... Xn, Yn2, error of Yn2 // End of the second set of measurements X1, Y13, error of Y13, X2, Y23, error of Y23, ... Xn, Yn3, error of Yn3 // end of the third set of measurements BTW, you'll need to format you're data like: X Y_i Measurement 1 1 1.1 1 1 1.2 ... 2 1 1.2 2 1 1.5 ... 3 1 1.3 3 1 1.1 ... 1 2 1.1 ... In other words, you need all the measurements for each group+measurement category rather than means and errors because the function calculates the mean and CI for you. HTH Manuel __ 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] issue with plot (type=h)
Hi Gašper On Wed, 2006-02-22 at 14:12 +0100, Gasper Cankar wrote: Hello everyone. For reasons too long to explain I wanted to do plots similar to histograms with plot(type=h). I ran into a problem - if I set line width too high, histogram isn't accurate anymore. For example: par(lend=2) plot(c(2,4,3,2),ylim=c(0,5), type=h) abline(h=3) Column 3 appears just as high as it should. But if I do par(lend=2) plot(c(2,4,3,2),ylim=c(0,5), type=h,lwd=100) abline(h=3) then columns become too high. Can I correct the problem or is there another way to display my data correctly? You need to use lend=1 or lend=butt in your par() statement. In my view, it would be nice to change the default to use lend=1 for plot type = h, or at least to include a warning when square is used, since the effect of increasing the lwd may not always be obvious. __ 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] R xyplot and background color
On Tue, 2006-02-21 at 12:49 -0800, Bryan Sykes wrote: Hi: I have tried (unsuccessfully) to change the default background color for my xyplot. I have used trellis.device(bg = white, new = F) and par(bg=white) before my xyplot command. Yet the color of the background has not changed. Is there something else I need to do? I am using a windows based machine to do this. Try xyplot(..., par.settings=list(background=white)) HTH, Manuel __ 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] Nesting Functions
On Thu, 2006-01-26 at 21:55 -0500, Duncan Murdoch wrote:
On 1/26/2006 9:45 PM, Manuel Morales wrote:
Dear list members,
I'm looking for a way to write nested functions similar to the
function Nest or NestList in Mathematica.
E.g.,
f-function(x) x+2*x
f(f(f(2)))
might instead be written as nest(f, 2, 3)
read as, nest function f 3 times with 2 as the initial value.
It's easy enough using a for loop:
nest - function(f, initial, reps) {
result - initial
for (i in seq(len=reps)) result - f(result)
result
}
Duncan Murdoch
That works, thanks! But what if I want to apply the function to a set of
vectors.
init.values-c(3,10,20)
rep.values-c(0,1,2)
nest(f,init.values,rep.values) fails because only the first value is
used in a for loop. The following works, but it's clunky and doesn't
scale with variation in the number of reps.
nest.vectorize-function(f, initial, reps)
ifelse(reps==0,initial,
ifelse(reps==1,f(initial),f(f(initial
nest.vectorize(f,init.values,rep.values)
Any suggestions?
Thanks again,
Manuel
__
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] Nesting Functions
Dear list members, I'm looking for a way to write nested functions similar to the function Nest or NestList in Mathematica. E.g., f-function(x) x+2*x f(f(f(2))) might instead be written as nest(f, 2, 3) read as, nest function f 3 times with 2 as the initial value. Thanks! Manuel __ 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.fit in predict.nls
Sorry to be so persistent but I really need to get
some measure of the error in the predictions of my nls
model.
I've tried to find out what predict.nls does and I've
got down to
MiModelo$m$predict
function (newdata = list(), qr = FALSE)
{
eval(form[[3]], as.list(newdata), env)
}
environment: 0x88a076c
But I can not find what is form.
Any help, please.
Manuel
Manuel Gutierrez wrote:
The option se.fit in predict.nls is currently
ignored.
Is there any other function available to calculate
the
error in the predictions?
Thanks,
Manuel
__
LLama Gratis a cualquier PC del Mundo.
Llamadas a fijos y móviles desde 1 céntimo por minuto.
__
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] se.fit in predict.nls
The option se.fit in predict.nls is currently ignored. Is there any other function available to calculate the error in the predictions? Thanks, Manuel __ LLama Gratis a cualquier PC del Mundo. Llamadas a fijos y móviles desde 1 céntimo por minuto. __ 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 propagation
Are there any functions to do error propagation in R? I have done a search with little success. Any pointers to read about this topic are greatly welcomed. My specific problem is: I use a linear model (lm) to predict the biomass of an individual in a population, then I add up the biomass of all individuals to calculate the total mass of the population. I want to calculate the error in the final estimate of the total biomass. Also, I have another case where I use a nls model instead of the lm, but in nls se.fit is currently ignored as stated in the help page. Is there an alternative way to calculate the errors in the predictions from nls? Thanks, Manuel __ 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] addition of error terms
Are there any functions to do error propagation in R? I have done a search with little success. Any pointers to read about this topic are greatly welcomed. My specific problem is: I use a linear model (lm) to predict the biomass of an individual in a population, then I add up the biomass of all individuals to calculate the total mass of the population. I want to calculate the error in the final estimate of the total biomass. Also, I have another case where I use a nls model instead of the lm, but in nls se.fit is currently ignored as stated in the help page. Is there an alternative way to calculate the errors in the predictions from nls? Thanks, Manuel __ 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] R moodle module
Hello, I have recently started to use moodle (a course management system) for my classes. It has turn out to be a very interesting tool. I woke up this morning wondering about the possibility of adding a module so that the students could use R without needing to install it, just using a server session into moodle. I know that there are already several on-going projects about R and the web, but, the question is : Is anybody already developing such a module for moodle? If the answer is no, I shall put it in my want-to-do list, although I can not promise any inmediate result due to may work overload. Thanks in advance for your patience and your good efforts! Best wishes, Manuel Castejón Limas E-mail: manuel.castejon at unileon.es Área de Proyectos de Ingeniería Departamento de Ingeniería Eléctrica y Electrónica Universidad de León Spain PS: You may read more about moodle at http://moodle.org __ 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] Attributing values to matrix according to names
Dear R-helpers Apologies for the basic question, but I just got stuck: I would like to write values from a vector into array cells with the same names count[1:10] 10010 10014 10015 10017 10030 10080 10100 10230 10250 10280 0 0 0 0 0 1 1 0 2 0 data[1:10,,1] [,1] [,2] [,3] [,4] [,5] 10010 NA NA NA NA NA 10014 NA NA NA NA NA 10015 NA NA NA NA NA 10016 NA NA NA NA NA 10017 NA NA NA NA NA 10100 NA NA NA NA NA 10140 NA NA NA NA NA 10150 NA NA NA NA NA 10160 NA NA NA NA NA 10170 NA NA NA NA NA length(count) [1] 2842 dim(data) [1] 4667 5 10 My operation should result in data[1:10,,1] [,1] [,2] [,3] [,4] [,5] 100100 NA NA NA NA 100140 NA NA NA NA 100150 NA NA NA NA 10016 NA NA NA NA NA 100170 NA NA NA NA 101001 NA NA NA NA 10140 NA NA NA NA NA 10150 NA NA NA NA NA 10160 NA NA NA NA NA 10170 NA NA NA NA NA data[10014,1,1]-count[10014] works but data[names(count),1,1]-count[names(count)] Fails with Error: indexing outside limits. Many thanks for any help Manuel __ 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] Attributing values to matrix according to names
Dear Peter and Uwe Thanks for your suggestions. However, data[names(count),1,1] - count[names(count)] Still gives the indexing problem, guess because not all element of count can be found in data. Found a way round this by temp-rep(NA, times=as.numeric(names(count[length(count)]))) temp[as.numeric(names(count))]-count rwname-rownames(data) for (i in 1:dim(data)[1]) data[i,1,1]-temp[as.numeric(rwname[i])] What works for me but I am convinced there is a far more elegant way. Kind regards Manuel -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Peter Dalgaard Sent: Thursday, October 20, 2005 4:29 PM To: Schneider, Manuel Cc: r-help@stat.math.ethz.ch Subject: Re: [R] Attributing values to matrix according to names Schneider, Manuel [EMAIL PROTECTED] writes: data[10014,1,1]-count[10014] works but data[names(count),1,1]-count[names(count)] Fails with Error: indexing outside limits. Well, you don't have a row named names(count) now do you? Try dropping the quotes. -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ 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-standard characters in Ascii-Files
Dear Henrik, dear Peter
many thanks for your reply. I've now done tests on other computers and yes,
this is a local problem on my computer. Unfortunately, my second HD produces
the same flaw so it must be the disk controller or whatever. When correcting
the files with vedit, the non-standard characters sometimes seem to jump around
to other places. Funny bug, but definitely nothing to do with R.
Best regards
Manuel
-Ursprüngliche Nachricht-
Von: Henrik Bengtsson [mailto:[EMAIL PROTECTED]
Gesendet: Montag, 29. August 2005 15:43
An: Peter Dalgaard
Cc: Schneider, Manuel; r-help@stat.math.ethz.ch
Betreff: Re: [R] Non-standard characters in Ascii-Files
Peter Dalgaard wrote:
Schneider, Manuel [EMAIL PROTECTED] writes:
Dear R-list
In R 2.1.1 under Win XP on a P4 with 2GB Ram when typing
temp-matrix(c(1:1600),4000,4000)
write(file=temp.txt, temp)
scan(temp.txt)
I receive:
Error in scan(temp.txt) : scan() expected 'a real', received
'414851'
The motivation for evoquing this meassage is that I am getting the
same meassage with exported Ascii-Files from the GIS. The files
contain very few, randomly scattered non-standard Ascii-characters.
This seems to be a local problem on my machine but I do not have a
clue on the reason (OS, Memory, HD?) nor who to ask. So, my apologies
for misusing this list and many thanks for any suggestion.
I tried this on a Linux box (with a somewhat outdated R version
though), and apart from eating up memory and disk space, nothing
untoward seems to happen:
temp-matrix(c(1:1600),4000,4000)
write(file=/tmp/temp.txt, temp)
dummy - scan(/tmp/temp.txt)
Read 1600 items
and on my R v2.1.1 patched (2005-08-25) on WinXP Pro SP2 (sic!), I get
temp-matrix(c(1:1600),4000,4000)
write(file=temp.txt, temp)
file.info(temp.txt)$size
[1] 136088897
rm(temp)
dummy - scan(temp.txt)
Read 1600 items
I'd suspect your harddisk or the disk controller...
I second this, check the file with an external application or try the
following ad-hoc code:
zcan - function(filename) {
fh - file(filename, open=r);
on.exit(close(fh));
count - 0;
while(TRUE) {
s - readChar(fh, n=1024);
if (nchar(s) == 0)
break;
count - count + nchar(s);
if (gsub([\n 0-9]*, , s) != )
stop(Error after reading , count, characters: , s);
}
}
Cheers
Henrik
__
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] Non-standard characters in Ascii-Files
Dear R-list In R 2.1.1 under Win XP on a P4 with 2GB Ram when typing temp-matrix(c(1:1600),4000,4000) write(file=temp.txt, temp) scan(temp.txt) I receive: Error in scan(temp.txt) : scan() expected 'a real', received '414851' The motivation for evoquing this meassage is that I am getting the same meassage with exported Ascii-Files from the GIS. The files contain very few, randomly scattered non-standard Ascii-characters. This seems to be a local problem on my machine but I do not have a clue on the reason (OS, Memory, HD?) nor who to ask. So, my apologies for misusing this list and many thanks for any suggestion. Kind regards Manuel __ 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] path analysis
Someone knows if it is possible to perform a path analysis with sem package (or any other) to explain a dependent *dichotomus* variable? Thanks, Manel Salamero __ 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] Console not found
I played around with memory limits in R 2.1.0 under XP in order to be able to work with large matrixes (3600x4100). Among several things I tried was to alter console settings and saving them. Since then, I can't restart Rgui. It says several times 'Console not found' with pieces of the text that usually appears in the console and then crashes. Rterm.exe works fine. I've now unistalled R 2.1.0 and installed R 2.1.1 with no effect, still console is not found. Any clues on this? Best regards Manuel °°° Manuel Schneider Eawag Environmental chemistry Ueberlandstr. 133 8600 Dübendorf Phone +41 44 823 51 18 __ 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(): Levenberg-Marquardt, Gauss-Newton, plinear - PI curve fitting
On Tue, 2005-06-21 at 06:57 -0400, Gabor Grothendieck wrote:
On 6/21/05, Christfried Kunath [EMAIL PROTECTED] wrote:
Hello,
i have a problem with the function nls().
This are my data in k:
V1V2
[1,]0 0.367
[2,] 85 0.296
[3,] 122 0.260
[4,] 192 0.244
[5,] 275 0.175
[6,] 421 0.140
[7,] 603 0.093
[8,] 831 0.068
[9,] 1140 0.043
With the nls()-function i want to fit following formula whereas a,b, and c
are variables: y~1/(a*x^2+b*x+c)
With the standardalgorithm Newton-Gauss the fitted curve contain an peak
near the second x,y-point.
This peak is not correct for my purpose. The fitted curve should descend
from the maximum y to the minimum y given in my data.
The algorithm plinear give me following error:
phi function(x,y) {
k.nls-nls(y~1/(a*(x^2)+b*x+c),start=c(a=0.0005,b=0.02,c=1.5),alg=plinear)
coef(k.nls)
}
phi(k[,1],k[,2])
Error in qr.solve(QR.B, cc) : singular matrix `a' in solve
I have found in the mailinglist
https://stat.ethz.ch/pipermail/r-help/2001-July/012196.html; that is if t
he data are artificial. But the data are from my measurment.
The commercial software Origin V.6.1 solved this problem with the
Levenberg-Marquardt algorithm how i want.
The reference results are: a = 9.6899E-6, b = 0.00689, c = 2.72982
What are the right way or algorithm for me to solve this problem and what
means this error with alg=plinear?
Thanks in advance.
This is not a direct answer to your question but log(y) looks nearly linear
in x when plotting them together and log(y) ~ a + b*x or
y ~ a*exp(b*x) will always be monotonic. Also, this model uses only 2
rather than 3 parameters.
If you want to use the original model, you could also use optim() and
specify your own minimization function. The default algorithm is slow
but has worked well for me where the faster algorithms in nls() have
choked. I adapted one of my functions for your data below:
minimize.fn-function(parms) {
fit.model-function(parms) {
y.pred={}; for (i in 1:9) {
a=parms[1];b=parms[2];c=parms[3];
y=1/(a*V1[i]^2+b*V1[i]+c)
y.pred=append(y.pred,y)}
list(fitted.values=y.pred)}
minimize.model-function(parms) {
Resids=V2-fit.model(parms)$fitted.values
ss=sum(Resids^2)
list(c(ss))}
fit=optim(c(parms[1],parms[2],parms[3]),minimize.model,
control=list(trace=1,maxit=1))
results=fit.model(c(fit$par[1],fit$par[2],fit$par[3]))
list(parms=fit$par,SS=fit$value,residuals=results$fitted.values-V2,
fitted.values=results$fitted.values)}
result-minimize.fn(c(0.0005,.02,1.5))
result-minimize.fn(result$parms)
This gives the results:
result$parms
[1] 1.184172e-05 4.878992e-03 3.045663e+00
result$SS
[1] 0.005436355
Which is very similar to what you want.
HTH,
Manuel
__
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] Increasing Console Paste Buffer
Hello list. I'm using R from the gnome-terminal in Fedora. My preference is to write programs in VIM, and then source the file from R, or copy and paste the lines into the console. I'm wondering if there is a way to increase the paste buffer as an alternative to sourcing large analyses. As was mentioned in a recent thread on Linux GUI's, I find that if I paste in a large amount of text, the lines end up getting cut off at some point. I wonder if this is an R restriction, because it seems like I am able to paste substantially more text in other console-based programs. Is there any way to increase the amount of text that I can paste into an R session? Thanks! Manuel __ 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] Iterative process for reading in text files
Hi Meredith,
When I've wanted to do this, I put all my files (group1.txt,
group2.txt ...) in a separate directory. Then, running R from that
directory:
files-list(files)
for (i in 1:length(files)) {
group-read.table(files[i],header=T)
do stats here
}
On Fri, 2005-04-29 at 15:17 +1000, Briggs, Meredith M wrote:
Hello
Instead of reading in group1.txt I want to read in groups1 for the first
iteration of i, then groups2 for the second and so on. Obviously I can't use
groups(i) but assume there is a way to do this.
group-read.table(C:/Data/April 2005/group1.txt,header=T)
thanks in advance
Meredith
__
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] Problem with R-2.1.0: install.packages() doesn't work
On Fri, 2005-04-22 at 15:42 -0500, Marc Schwartz wrote: On Fri, 2005-04-22 at 15:44 -0400, Manuel Morales wrote: On Fri, 2005-04-22 at 17:48 +0200, Uwe Ligges wrote: Waichler, Scott R wrote: I installed R-2.1.0 from source on a Linux box running Red Hat Enterprise Linux WS release 4 but install.packages() wouldn't work (see below). install.packages(rgenoud) --- Please select a CRAN mirror for use in this session --- Error in inherits(x, factor) : Object res not found Quite probably you have no X11 connection to this machine. R tries to ask you which CRAN mirror you are going to choose, and it fails to present the tcltk window. You might want to call chooseCRANmirror(graphics=FALSE) and setRepositories(graphics=FALSE) prior to install.packages(). Uwe Ligges I have the same problem after building R-2.1.0 from source on Fedora Core 3. The suggestion above fixes this, but what do you mean by Quite probably you have no X11 connection on this machine? I'm guessing you don't mean that X11 is not running (I use Gnome for my desktop). Manuel For both Scott and Manuel, Can you post back with the output of: capabilities() This is what I got: capabilities() jpeg pngtcltk X11 http/ftp sockets libxml fifo TRUE TRUEFALSE TRUE TRUE TRUE TRUE TRUE cledit IEEE754iconv TRUE TRUE TRUE I recompiled after downloading the tc and tk development packages, which gave tcltk as TRUE. update.packages() works if tcltk is enabled, but it seems not to revert to the non-graphical interface otherwise... __ 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] Problem with R-2.1.0: install.packages() doesn't work
On Fri, 2005-04-22 at 17:48 +0200, Uwe Ligges wrote: Waichler, Scott R wrote: I installed R-2.1.0 from source on a Linux box running Red Hat Enterprise Linux WS release 4 but install.packages() wouldn't work (see below). install.packages(rgenoud) --- Please select a CRAN mirror for use in this session --- Error in inherits(x, factor) : Object res not found Quite probably you have no X11 connection to this machine. R tries to ask you which CRAN mirror you are going to choose, and it fails to present the tcltk window. You might want to call chooseCRANmirror(graphics=FALSE) and setRepositories(graphics=FALSE) prior to install.packages(). Uwe Ligges I have the same problem after building R-2.1.0 from source on Fedora Core 3. The suggestion above fixes this, but what do you mean by Quite probably you have no X11 connection on this machine? I'm guessing you don't mean that X11 is not running (I use Gnome for my desktop). Manuel __ 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] large matrix
Dear R-users
I need to convert a matrix with three columns in a new array with
multiple columns.
For example,
oldmatrix
1 4 5
1 54 52
1 9 43
2 32 5
2 54 6
2 76 6
3 54 54
3 543 7
3 54 6
newmatrix
5 5 54
52 6 7
43 6 6
if the first column have a new value then add a column to the new
matrix and the new[i,j] - old[,3][i]
I write this code, but my initial matrix is very large and the
convertion is very slow. How I can optimise that code?
converter - function(oldmatrix) {
nr - sapply(oldmatrix[,1], nrow)
temporario - oldmatrix[,1][1]
j - 1
nn-1
for (i in seq(along=nr)) {
if (temporario == oldmatrix[,1][i]) {
nn - nn +1
}
else
{
j - j + 1
}
temporario - oldmatrix[,1][i]
}
nn-i/j
newmatrix-array(dim=c(nn, j))
temporario - oldmatrix[,1][1]
j - 1
i-1
for (i in seq(along=nr)) {
if (temporario == oldmatrix[,1][i]) {
newmatrix[,j][i] - oldmatrix[,3][i]
}
else
{
j - j + 1
}
temporario - oldmatrix[,1][i]
}
return(newmatrix)
}
Thanks a lot
Jorge
[[alternative text/enriched 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] dealing with multicollinearity
I have a linear model y~x1+x2 of some data where the coefficient for x1 is higher than I would have expected from theory (0.7 vs 0.88) I wondered whether this would be an artifact due to x1 and x2 being correlated despite that the variance inflation factor is not too high (1.065): I used perturbation analysis to evaluate collinearity library(perturb) P-perturb(A,pvars=c(x1,x2),prange=c(1,1)) summary(P) Perturb variables: x1 normal(0,1) x2 normal(0,1) Impact of perturbations on coefficients: mean s.d. min max (Intercept) -26.0670.270 -27.235 -25.481 x1 0.7260.0250.6720.882 x2 0.0600.0110.0370.082 I get a mean for x1 of 0.726 which is closer to what is expected. I am not an statistical expert so I'd like to know if my evaluation of the effects of collinearity is correct and in that case any solutions to obtain a reliable linear model. Thanks, Manuel Some more detailed information: A-lm(y~x1+x2) summary(A) Call: lm(formula = y ~ x1 + x2) Residuals: Min1QMedian3Q Max -4.221946 -0.484055 -0.004762 0.397508 2.542769 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) -27.234720.27996 -97.282 2e-16 *** x10.882020.02475 35.639 2e-16 *** x20.081800.01239 6.604 2.53e-10 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Residual standard error: 0.823 on 241 degrees of freedom Multiple R-Squared: 0.8411, Adjusted R-squared: 0.8398 F-statistic: 637.8 on 2 and 241 DF, p-value: 2.2e-16 cor.test(x1,x2) Pearson's product-moment correlation data: x1 and x2 t = -3.9924, df = 242, p-value = 8.678e-05 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: -0.3628424 -0.1269618 sample estimates: cor -0.248584 __ 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] nlsList vs nlme with lsoda function
Hello list members,
I'm trying to get nlme to work with lsoda from odesolve. Currently, I
can use nlsList to fit a simple logistic growth model to some simulated
data, but I get the following error message with nlme: Error in
model.frame(formula, rownames, variables, varnames, extras, extranames,
: invalid variable type
The commands I'm trying are below. Any suggestions for how to get this
to work?
Thanks!
This works:
fit.nlsList-nlsList(xvals~lsoda(x0,times,Logist,
+c(r=r,K=K,x0=x0))[,2]|group,
+start=list(r=1,K=500,x0=2),data=data.group)
This returns the error message above:
fit.nlme-nlme(xvals~(lsoda(x0,times,Logist,
+c(r=r,K=K,x0=x0))[,2]),
+start=list(r=1,K=500,x0=2),data=data.group,
+fixed=r+K+x0~1,random=x0~1)
The model function is:
Logist=function(t,x,parms) {
N1-x[1]
with(as.list(parms),{
dN1=r*N1*(1-(N1/K))
list(c(dN1))
})}
My data set looks like:
Grouped Data: xvals ~ times | group
xvals times group
1-0.44543969 0 1
238.86411972 3 1
3 310.77106961 6 1
4 486.78653327 9 1
5 0.40613173 0 2
634.94635643 3 2
7 307.81884870 6 2
8 486.46098417 9 2
etc...
__
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] Re: testing slopes different than a given value
Thanks to all, Yes, I meant a single test for both coefficients. Peter's reply is what I wanted. I've tried with linear.hypothesis but I must confess that with my limited statistical experience and without the car book at hand, the nomenclature for the function was a bit difficult to understand for me. A toy example of linear.hypothesis for my case would be great. Thanks, Manuel --- Peter Dalgaard [EMAIL PROTECTED] escribió: John Fox [EMAIL PROTECTED] writes: Dear Vito, Since Manuel says that he wants to obtain a test and not obtain two tests, I assume that he's interested in the F-test for the hypothesis that both coefficients are simultaneously equal to the specified values rather than in the t-tests for the individual hypotheses. Regards, John ...in which case one answer is this: y-3+0.6*x1+0.3*x2 + rnorm(100,sd=.1) # as meant, no? fm-lm(y~x1+x2) anova(fm, lm(y~offset(0.6*x1+0.3*x2))) Analysis of Variance Table Model 1: y ~ x1 + x2 Model 2: y ~ offset(0.6 * x1 + 0.3 * x2) Res.Df RSS Df Sum of Sq F Pr(F) 1 97 1.06118 2 99 1.06184 -2 -0.00066 0.0302 0.9703 -- O__ Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ 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 slopes different than a given value
In a multiple linear regression with two independent variables is there any function in R to test for the coefficients being different than some given values? Example: x1-rnorm(100) x2-rnorm(100) y-3+0.6*x1+0.3*x2 lm(y~x1+x2) Obtain a test for the coefficients for x1 being different than 0.6 and for x2 different than 0.3 Thanks Manuel __ 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] Rare Cases and SOM
I am trying to understand how the SOM algorithm works using library(class) SOM function. I have a 1000*10 matrix and I want to be able to summarize the different types of 10-element vectors. In my real world case it is likely that most of the 1000 values are of one kind the rest of other (this is an oversimplification). Say for example: InputA-matrix(cos(1:10),nrow=900,ncol=10,byrow=TRUE) InputB-matrix(sin(5:14),nrow=100,ncol=10,byrow=TRUE) Input-rbind(InputA,InputB) I though that a small grid of 3*3 would be enough to extract the patterns in such simple matrix : GridWidth-3 GridLength-3 gr - somgrid(xdim=GridWidth,ydim=GridLength,topo = hexagonal) test.som - SOM(Input, gr) par(mfrow=c(GridLength,GridWidth)) for(i in 1:(GridWidth*GridLength)) plot(test.som$codes[i,],type=l) Only when I use a larger grid (say for example 7*3 ) I get some of the representatives for the sin pattern. This must have something to do with the initialization of the grid, as the sin is so rare it is unlikely that I get it as a reference vector. Afterwards, because the selection for the training is also random it is also unlikely they are picked. I've been trying to modify some of the other parameters for the SOM also, but I would appreciatte some input to keep me going until I receive the reference books from my bookstore. Are my suspictions right? Should I be using the SOM for my study or should I look somewhere else? NOTE: I have no prior knowledge of whether the datasets I want to analyse will have rare cases or not or where they will be located. Thanks, Manuel __ 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] memory and swap space in ncdf
I've a linux system with 2Gb of memory which is not enough for reading a 446Mb netcdf file using ncdf: library(ncdf) ncold - open.ncdf(gridone.grd) Error: cannot allocate vector of size 1822753 Kb When I look at the free memory in my system I can see that none of the Swap space is being used by R. I am a newbie in linux and R, I've read the Memory help pages but still have some questions: can I use the swap space in R to solve my problem of lack of memory? if not, are there any ways to read the data apart from buying more RAM? Thanks, M __ 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] kde2d and borders
Hallo, I want to use kde2d to visualize data on a sphere given in spherical coordinates. Now the problem is, that phi == 2*pi = 0, so in principal I have to connect (in a graphical view) the left and right border of my plot (and the bottom and top). Has anyone any idea how to do that ? Thanks, Manuel -- - Manuel Metz Sternwarte der Universitaet Bonn Auf dem Huegel 71 (room 3.06) D - 53121 Bonn E-Mail: [EMAIL PROTECTED] Phone: (+49) 228 / 73-3660 Fax:(+49) 228 / 73-3672 __ 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] unlist kills R
When I try to unlist a very large list, R is killed without any other warning: A-as.list(as.data.frame(matrix(1:21639744,nrow=3578,ncol=6048))) with AA-unlist(A) or AA-c(A,recursive=TRUE) I get a R terminado (killed) and the end of the session. I think I'll need to get more RAM (now 1Gb, any other solutions welcomed) to be able to do this but, shouldn't I get a more gentle warning than the kill message? Thanks, Manuel platform i686-pc-linux-gnu arch i686 os linux-gnu system i686, linux-gnu status major2 minor0.1 year 2004 month11 day 15 language R __ 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] test multiple objects for being equal length
I could not find any help pages on How to test many objects for being of equal length Something like identical for more than two objects? x-1:6 y-1:10 z-3:5 ## For two objects I can do: identical(length(x),length(y)) ## For more than two I currently can do: length(unique(c(length(x),length(y),length(z==1 but there must be a better way. Thanks, M __ [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] augPred with lme(...,subset=...)
Hello, Is there a way to get augPred to work with lme if a subset statement is used? For example, if I modify the example from ?augPred.lme to include a subset statement, I get the following error: fm1 - lme(Orthodont, random = ~1, subset=distance19) augPred(fm1, length.out = 2, level = c(0,1)) Error in tapply(object[[nm]], groups, FUN[[numeric]], ...) : arguments must have same length Thanks! Manuel __ [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] substitute accents
I have an openoffice spreadsheet with a column of
character strings.
Some of them contain accents.
I want to read it in R so I have saved it as a csv
file using Western Europe (ISO-8859-1) character set
(the default, I've tried other sets but it doesn't
help).
R reads it fine with
CharMatrix-read.csv(test.csv,header=F,sep=,,as.is=TRUE);
Say I wan't to replace the 'o' with accent in the
first cell.
I've tried:
gsub('ó','o', CharMatrix[1,1])
But, It doesn't make any substitution
Trying to find a solution I input the character string
in R and do the substitution:
CharMatrix[1,1]-hóla
gsub('ó','o', CharMatrix[1,1])
And it works. I think the difference is that when I
now print the content of CharMatrix I get a \201
before the ó while I didn't get it with the openoffice
imported csv file.
I'm sure it is a problem with my understanding of how
accents can be specified. Can someone give me any
solutions / references?
Thanks,
M
_
platform i686-pc-linux-gnu
arch i686
os linux-gnu
system i686, linux-gnu
status
major2
minor0.0
year 2004
month10
day 04
language R
__
Renovamos el Correo Yahoo!: ¡100 MB GRATIS!
Nuevos servicios, más seguridad
__
[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] substitute accents
$ locale
LANG=en_GB
LC_CTYPE=en_GB
LC_NUMERIC=en_GB
LC_TIME=en_GB
LC_COLLATE=en_GB
LC_MONETARY=en_GB
LC_MESSAGES=en_GB
LC_PAPER=en_GB
LC_NAME=en_GB
LC_ADDRESS=en_GB
LC_TELEPHONE=en_GB
LC_MEASUREMENT=en_GB
LC_IDENTIFICATION=en_GB
LC_ALL=
$ locale charmap
ISO-8859-1
I have tried changing the locales with no difference.
Is this fine?
And, no, I didn't get any warning message.
My sistem is a debian sid under kde 3.3.
Thanks,
M
--- Prof Brian Ripley [EMAIL PROTECTED]
escribió:
Can you please tell us what locale you are working
in?
This looks as if the problem might be the use of a
UTF-8 locale, which R
does not currently support and which some Linux
distros have made their
default. However, R does issue a warning -- so did
you get one?
On Thu, 25 Nov 2004, Manuel Gutierrez wrote:
I have an openoffice spreadsheet with a column of
character strings.
Some of them contain accents.
I want to read it in R so I have saved it as a csv
file using Western Europe (ISO-8859-1) character
set
(the default, I've tried other sets but it doesn't
help).
R reads it fine with
CharMatrix-read.csv(test.csv,header=F,sep=,,as.is=TRUE);
Say I wan't to replace the 'o' with accent in the
first cell.
I've tried:
gsub('ó','o', CharMatrix[1,1])
But, It doesn't make any substitution
Trying to find a solution I input the character
string
in R and do the substitution:
CharMatrix[1,1]-hóla
gsub('ó','o', CharMatrix[1,1])
And it works. I think the difference is that when
I
now print the content of CharMatrix I get a \201
before the ó while I didn't get it with the
openoffice
imported csv file.
I'm sure it is a problem with my understanding of
how
accents can be specified. Can someone give me any
solutions / references?
Thanks,
M
_
platform i686-pc-linux-gnu
arch i686
os linux-gnu
system i686, linux-gnu
status
major2
minor0.0
year 2004
month10
day 04
language R
__
Renovamos el Correo Yahoo!: ¡100 MB GRATIS!
Nuevos servicios, más seguridad
__
[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
--
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
__
Renovamos el Correo Yahoo!: ¡100 MB GRATIS!
Nuevos servicios, más seguridad
__
[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] ML vs. REML with gls()
Hello listmembers,
I've been thinking of using gls in the nlme package to test for serial
correlation in my data set. I've simulated a sample data set and have
found a large discrepancy in the results I get when using the default
method REML vs. ML.
The data set involves a response that is measured twice a day (once for
each level of a treatment factor). In my simulated data set, I have a
mean day effect but no serial correlation between days, otherwise.
However, if I include an AR1 structure in my model, I find significant
evidence of serial correlation when I fit the model using maximum
likelihood! I don't see the same result when I use restricted maximum
likelihood. I'm not a statistician, but I'm wondering if this is
expected. I also realize that this particular example would work well
with day as a random effect, but then I'm not sure how I would test for
serial correlation across days. I have Pinheiro and Bate's extremely
useful book, but I can't seem to find anything that addresses this. I'm
using R 1.9.1 with version 3.1-50 of nlme.
Thanks very much for any help. I've included my sample code below.
Manuel Morales
#Below I generate a day variable (note two trials per day).
day-rep(1:25,each=2)
#Below I generate a random day effect
day.effect-c(); for(i in 1:25) {
tmp-rnorm(1);
day.effect-append(rep(tmp,2),day.effect)}
#Below I randomize the application of the treatment effect w/i days.
trt-c(); for(i in 1:25) {
if (rnorm(1)0) trt-append(c(0,1),trt)
else trt-append(c(1,0),trt)}
#Below I generate the response variable. Trt increases the response.
resp-trt+day.effect+rnorm(50)
#Below I analyze the data using REML or ML w/ and w/o AR1 errors.
library(nlme)
base.gls-gls(resp~factor(trt)+factor(day))
ar.gls-gls(resp~factor(trt)+factor(day),corr=corAR1())
base.gls.ml-gls(resp~factor(trt)+factor(day),method=ML)
ar.gls.ml-gls(resp~factor(trt)+factor(day),corr=corAR1(),method=ML)
#Below I compare models using REML or ML
anova(base.gls,ar.gls); anova(base.gls.ml,ar.gls.ml)
__
[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] interaction plot with error bars based on TukeyHSD intervals
Hello,
maybe the previous message was not very clear. Thus, I will try to
explain the problem in a better way.
I would like to have an interaction plot with error bars based on
TukeyHSD intervals.
In an experiment with two factors A and B, each of them with two
levels A={a1,a2} and B={b1,b2}, I would like to do a plot like the
following:
attached image prueba.png
I calculate the corresponding interval with:
tk - TukeyHSD(aov(X$response ~ (factor(X$A) + factor(X$B))2, data=X) ,
conf.level=0.95)
HSDfactor -
max(abs(tk$factor(X$A):factor(X$B)[,2]-tk$factor(X$A):factor(X$B)[,3]))
Finally, I modified interaction.plot() to add error bars with:
..
++ ylim - c(min(cells)-(HSDfactor*0.5), max(cells)+(HSDfactor*0.5))
matplot(xvals, cells, ..., type = type, xlim = xlim, ylim = ylim,
xlab = xlab, ylab = ylab, axes = axes, xaxt = n,
col = col, lty = lty, pch = pch)
++ ly - cells[,1]+(HSDfactor*0.5)
++ uy - cells[,1]-(HSDfactor*0.5)
++ errbar(xvals,cells[,1],ly,uy,add=TRUE, lty=3, cap=0, lwd=2)
++ ly - cells[,2]+(HSDfactor*0.5)
++ uy - cells[,2]-(HSDfactor*0.5)
++ errbar(xvals,cells[,2],ly,uy,add=TRUE, lty=3, cap=0, lwd=2)
if(legend) {
yrng - diff(ylim)
yleg - ylim[2] - 0.1 * yrng
.
However, the resulting intervals are much bigger than they should be,
thus I think I did something wrong.
I have searched on Google, CRAN and the R mail archive but I only found
related problems [1, 2], but not any answer...
Andrew Robinson [1] shows that there is a problem when using TukeyHSD
for interaction terms. However, nobody suggested any work-around.
In another different thread, Martin Henry H. Stevens [2] suggests a
work-around for the problem, but he thinks that the results obtained are
slightly inaccurate. As well, there is not feedback to solve the problem.
Any idea?
Thank you very much.
Manuel.
[1] http://finzi.psych.upenn.edu/R/Rhelp02/archive/12926.html
[2] http://finzi.psych.upenn.edu/R/Rhelp02/archive/32849.html
inline: prueba.png__
[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] interaction plot with intervals based on TukeyHSD
Hi,
The problem is that I would like to do an interaction plot with
intervals based on Tukey's honestly significant difference (HSD)
procedure, but I do not know how to do it in R.
I have 3 factors A, B and C and a response variable response.
I would like to study a model where there are main effects and second
order interaction effects.
For instance, for factors B and C, I would like to plot a line for each
level of C connecting the least square means for every level of B.
Additionally, I would like to include the 95,0% HSD intervals for
the means in such a way that any two intervals which do not overlap
correspond to a
pair of means which have a statistically significant difference.
For comparison, in STATGRAPHICS Plus, I choose Analysis of Variance -
Multifactor ANOVA...
Then I plot the second order interactions and I have the option to plot
intervals. I choose to plot Tukey HSD intervals at 95% and I obtain an
interaction plot like using interaction.plot on R, but addittionally an
interval is plotted on every point.
How can I do this on R?
I have already looked in Google, in the mail archive and in some books,
but I didn´t find the answer.
I tried to calculate the intervals using:
|tk - TukeyHSD(aov(X$response ~ (factor(X$A) + factor(X$B) +
factor(X$C))^2, data=X) , conf.level=0.95)
HSDfactor1 -
max(abs(tk$factor(X$A):factor(X$B)[,2]-tk$factor(X$A):factor(X$B)[,3]))
HSDfactor2 -
max(abs(tk$factor(X$A):factor(X$C)[,2]-tk$factor(X$A):factor(X$C)[,3]))
HSDfactor3 -
max(abs(tk$factor(X$B):factor(X$C)[,2]-tk$factor(X$B):factor(X$C)[,3]))
And I modified the function interaction.plot() adding the following
lines of code:
...
*++ ylim - c(min(cells)-(HSDfactor*0.5), max(cells)+(HSDfactor*0.5))*
matplot(xvals, cells, ..., type = type, xlim = xlim, ylim = ylim,
xlab = xlab, ylab = ylab, axes = axes, xaxt = n,
col = col, lty = lty, pch = pch)
*++ ly - cells[,1]+(HSDfactor*0.5)
++ uy - cells[,1]-(HSDfactor*0.5)*
*++ errbar(xvals,cells[,1],ly,uy,add=TRUE, lty=3, cap=0, lwd=2)
++ ly - cells[,2]+(HSDfactor*0.5)
++ uy - cells[,2]-(HSDfactor*0.5)
++ errbar(xvals,cells[,2],ly,uy,add=TRUE, lty=3, cap=0, lwd=2)*
if(legend) {
yrng - diff(ylim)
yleg - ylim[2] - 0.1 * yrng
.
|Finally, I call this modified function as:
|interaction.plot2( factor(X$A), factor(X$B), response,las=3, type=b,
HSDfactor = HSDfactor1, lwd=3)
|However, the resulting intervals are much bigger than the ones
calculated by STATGRAPHICS, thus I think I did something wrong.
I am not an expert in Statistics nor in R, thus if anyone has any
suggestion...
Thank you very much,
Manuel.
[[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] persp(), axis font size
Is there a way to adjust the font size for axis labels when using persp()? The parameter cex works for adjusting the global font size, but I can't seem to make cex.lab or cex.axis work for adjusting these values independently. Or, is there a preferred method for making surface plots in R? I'm using R version 1.8. Thanks, Manuel __ [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] intervals in lme() and ill-defined models
There has been some recent discussion on this list about the value of using intervals with lme() to check for whether a model is ill-defined. My question is, what else can drive very large confidence intervals for the variance components (or cause the error message Error in intervals.lme(Object) : Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance). To illustrate my question, I use examples from the book Mixed-Effects-Models in S and S-PLUS by Pinheiro and Bates, and from an analysis of my own data. In chapter 1, Pinheiro and Bates show that if you use a model with an interaction in the random effects term without appropriate replication in the data, the model will appear to fit but the confidence intervals for the variance components will be very large. They suggest using intervals() as a check that the model is appropriately defined: test.1 - lme(effort~Type, data=ergoStool, random=~1|Subject) test.2 - lme(effort~Type, data=ergoStool, random=~1|Subject/Type) intervals(test.2) Random Effects: Within-group standard error: lower est.upper 1.054760e-07 4.599834e-01 2.005999e+06 In fact, using anova() to compare these two models shows that nothing is gained by adding the interaction: anova(test.1,test.2) Model df AIC BIClogLik Test L.Ratio p-value test.1 1 6 133.1308 141.9252 -60.56539 test.2 2 7 135.1308 145.3909 -60.56539 1 vs 2 0 1 HOWEVER, for the example in chapter 5.3 of the book in which an autoregressive structure is used for the within group errors, I get the following error: test - lme(follicles~sin(2*pi*Time)+cos(2*pi*Time),data=Ovary,random=pdDiag(~sin(2* pi*Time))) test.ar1 - lme(follicles~sin(2*pi*Time)+cos(2*pi*Time),data=Ovary,random=pdDiag(~sin(2* pi*Time)),correlation=corAR1()) intervals(test.ar1) Error in intervals.lme(test.ar1) : Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance BUT, anova appropriately selects the autoregressive model as best: anova(test,test.ar1) Model df AIC BIClogLik Test L.Ratio p-value test1 6 1638.082 1660.404 -813.0409 test.ar127 1564.445 1590.487 -775.2224 1 vs 2 75.637 .0001 In the book, intervals() DOES appear to work, but the authors are using S-PLUS. My concern is that when I try to fit the following two models to my own data, I get very large confidence intervals for the within-subject error even thought AIC selects the autoregressive model as best: result - lme(log(T1+1)~factor(trt1)*factor(trt2)*factor(Census),data=data,random=~1|B lock/Subject) result.ar1 - lme(log(T1+1)~factor(trt1)*factor(trt2)*factor(Census),data=data,random=~1|B lock/Subject,correlation=corAR1()) intervals(result.ar1) Random Effects: lower est. upper sd((Intercept)) 3.491934e-13 0.01461032 611299013 Correlation structure: lower est. upper Phi 0.6543028 0.7574806 0.8329729 anova(result,result.ar1) Model df AIC BIClogLik Test L.Ratio p-value result 1 19 518.1501 581.5633 -240.0750 result.ar1 2 20 475.9776 542.7283 -217.9888 1 vs 2 44.17249 .0001 Why are my within-subject errors so large, and why doesn't intervals() work for the autoregressive errors example in the book. I am using R version 1.8 on Windows 2000. Any insight would be greatly appreciated! Manuel Manuel A. Morales Assistant Professor, Biology Williams College Williamstown, MA 01267 ph: 413-597-2983 | fax: 413-597-3495 http://mutualism.williams.edu __ [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] PROC MIXED vs. lme()
I'm trying to learn how to do a repeated measures ANOVA in R using lme(). A data set that comes from the book Design and Analysis has the following structure: Measurements (DV) were taken on 8 subjects (SUB) with two experimental levels (GROUP) at four times (TRIAL). In SAS, I use the code: PROC MIXED DATA=[data set below]; CLASS sub group trial; MODEL dv = group trial group*trial; REPEATED trial / SUBJECT=sub TYPE=CS; run; which gives the results: Tests of Fixed Effects SourceNDF DDF Type III F Pr F GROUP 1 62.51 0.1645 TRIAL 318 22.34 0.0001 GROUP*TRIAL 3180.58 0.6380 In R, I'm trying the code: results.cs - lme(DV ~ factor(GROUP)*factor(TRIAL), data=[data set below], random= ~factor(TRIAL)|SUB, correlation=corCompSymm() ) anova(results.cs) which gives the results: numDF denDF F-value p-value (Intercept) 118 3383.953 .0001 factor(GROUP) 1 64.887 0.0691 factor(TRIAL) 318 239.102 .0001 factor(GROUP):factor(TRIAL) 3181.283 0.3103 Why are these results different? I'm a newbie to R, have the book Mixed Effects Models in S and S-Plus, but can't seem to get this analysis to work. Any suggestions? Thanks! Manuel Data: SUB GROUP DV TRIAL 1 1 3 1 1 1 4 2 1 1 7 3 1 1 3 4 2 1 6 1 2 1 8 2 2 1 12 3 2 1 9 4 3 1 7 1 3 1 13 2 3 1 11 3 3 1 11 4 4 1 0 1 4 1 3 2 4 1 6 3 4 1 6 4 5 2 5 1 5 2 6 2 5 2 11 3 5 2 7 4 6 2 10 1 6 2 12 2 6 2 18 3 6 2 15 4 7 2 10 1 7 2 15 2 7 2 15 3 7 2 14 4 8 2 5 1 8 2 7 2 8 2 11 3 8 2 9 4 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] VIP and pls
Hi everyone, does anyone know a way to calculate VIP (Variable Importance in the Projection) from R pls.pcr package procedures? Thanks, -- Manuel Martin 23 rue des gâtines 75020 Paris tel: 0143664965 _ Envie de discuter en live avec vos amis ? Télécharger MSN Messenger http://www.ifrance.com/_reloc/m la 1ère messagerie instantanée de France __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Re: Mixed effects with nlme
Hi, R-users: Last week I send a request for help to this list. I have receive until now two kindly responses from Spencer Graves and Renauld Lancelot. They both point interesting things to fit an adequate model to my data but unfortunately it persists without a satisfactory solution. I propose again the same problem but using with a different dataset (Assay), taken from Pinheiro and Bates' book on page 163, that is relevant with crossed random effects. I have fitted the same model (p. 165) fmAssay - lme(logDens ~ sample * dilut, Assay, random=pdBlocked(list(, pdIdent(~1), pdIdent(~sample-1),pdIdent(~dilut-1))) ) and the results with anova function (p. 166) are numDF denDF F-value p-value (Intercept) 129 537.6294 .0001 sample 529 11.2103 .0001 dilut429 420.5458 .0001 sample:dilut2029 1.6072 0.1192 The problem is that with this approach one obtains correct F-values, but using a common residual term for DenDF that is a combination of (Block + Block:sample + Block:dilut). Then probability values are different to that obtained when we used the classical AOV funtion to fit the same model, because in this case each term is mapped with a error term (so sample uses Block:sample, dilut uses Block:dilut, and sample:dilut uses Block:sample:dilut): mod-aov(logDens ~ sample*dilut + Error(Block+Block/sample+Block/dilut), data=Assay) summary(mod) Error: Block DfSum Sq Mean Sq F value Pr(F) Residuals 1 0.0083115 0.0083115 Error: Block:sample Df Sum Sq Mean Sq F value Pr(F) sample 5 0.276153 0.055231 11.213 0.009522 Residuals 5 0.024627 0.004925 Error: Block:dilut Df Sum Sq Mean Sq F valuePr(F) dilut 4 3.7491 0.9373 420.79 1.684e-05 Residuals 4 0.0089 0.0022 Error: Within Df Sum Sq Mean Sq F value Pr(F) sample:dilut 20 0.055525 0.002776 1.6069 0.1486 Residuals20 0.034555 0.001728 Obviously, the interest on linear mixed effects is open with the possibility of modeling with correlated and/or heterocedastic errors, and this end cannot be pursued with AOV function. Summarizing, the problem is that I have not found a way to obtain with NLME the same results (DF, F-ratios and probabilities) that I get with AOV and multistratum errors. Is this an inconvenience of program?, probably due to the impossibility to use multiple nested arguments as random(~1|Block/sample+dilut) or random(~1|Block/sample*dilut) SAS MIXED can also fit these data and obtain correct results by means of a combination of random and repeated arguments: model = sample dilut sample*dilut; random = Block sample*Block dilut*Block; repeated /type=cs Sub=Block; Type 3 Tests of Fixed Effects Num Den Effect DF DFF ValuePr F sample 5 5 11.210.0095 dilut 4 4 420.79.0001 sample*dilut 20 20 1.610.1486 May be possible obtain the same results with NLME function? I would appreciate any kind of help. Best regards, Manuel Ato University of Murcia Spain e-mail: [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] mixed effects with nlme
Dear R users: I have some difficulties analizing data with mixed effects NLME and the last version of R. More concretely, I have a repeated measures design with a single group and 2 experimental factors (say A and B) and my interest is to compare additive and nonadditive models. suj rvA B 1 s1 4 a1 b1 2 s1 5 a1 b2 3 s1 7 a1 b3 4 s1 1 a2 b1 5 s1 4 a2 b2 6 s1 2 a2 b3 7 s2 6 a1 b1 8 s2 8 a1 b2 9 s2 10 a1 b3 10 s2 3 a2 b1 11 s2 6 a2 b2 12 s2 6 a2 b3 13 s3 1 a1 b1 14 s3 6 a1 b2 15 s3 5 a1 b3 16 s3 3 a2 b1 17 s3 5 a2 b2 18 s3 4 a2 b3 19 s4 2 a1 b1 20 s4 10 a1 b2 21 s4 12 a1 b3 22 s4 1 a2 b1 23 s4 4 a2 b2 24 s4 7 a2 b3 25 s5 5 a1 b1 26 s5 10 a1 b2 27 s5 10 a1 b3 28 s5 5 a2 b1 29 s5 6 a2 b2 30 s5 5 a2 b3 31 s6 1 a1 b1 32 s6 7 a1 b2 33 s6 8 a1 b3 34 s6 2 a2 b1 35 s6 8 a2 b2 36 s6 7 a2 b3 It is very easy to fit these data with base R function AOV: NonAdditive model: aov(rv ~ A*B + Error(suj+suj/A+suj/B) Additive model: aov(rv ~ A*B + Error(suj) and also easy with SAS MIXED (I missed some obvious lines): NonAdditive model model vr = A B A*B; random suj A*suj B*suj; repeated / type=cs subj=suj; Additive model; model vr = A B A*B /ddfm=satterth; repeated / type=cs subj=suj; Using LME I do not find any problems to fit the additive model with lme(vr~A*B, random=~1|suj, cor=corCompSymm()) but I have found some difficulties fitting the nonadditive model. Can anyone help me? Thanks in advance. Manuel Ato Dpto. Psic.Básica y Metodología Apartado 4021 30080 MURCIA (Spain) e-mail: [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] rterm
I would like to use Rterm but i don`t know its parameters. I have searched about this issue but i haven´t found anything. thank you. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] re: box counting method and other landscape ecology measures
Dear Rohan, Have a look at the fdim library, it may be of interest to you as far as fractal dimension (or box counting if you prefer) is concerned. Best wishes, Manuel Castejon __ [EMAIL PROTECTED] mailing list http://www.stat.math.ethz.ch/mailman/listinfo/r-help
