[R] How Does One Use the Value of Density Function?
How do people usually use the result of density function (e.g. dnorm)? Especially when its value can be greater than 1. What do they do with such density 1? dnorm(2.02,2,.24) [1] 1.656498 - G.V. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] List of tags in roxygen and use for S4 classes?
On Mon, Sep 7, 2009 at 9:41 AM, Peter Danenbergp...@roxygen.org wrote: I would also wish for a better (online) documentation, as I think the general idea of roxygen is great. I agree completely. Good call; the vignette is terse and outdated. Manuel and I are in the process of preparing a paper based on our DSC talks; that should fill in some of the details w.r.t. S4. That sounds exciting. I am really looking forward to it - because roxygen is the way to go for documenting. Rainer -- Rainer M. Krug, Centre of Excellence for Invasion Biology, Stellenbosch University, South Africa __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Scan and read.table
Dear R experts.. I am trying to read data-sections in a large consolidated dataset, containing section headers and the data . There are many options available to implement, I was wondering what optimal function, to extract section headers and data (w/ columns), could be used on the dataset that looks like as provided at the end of this email? In each section of a dataset, 1st line of the section is a table title, followed by names of the columns and rows of data. TABLE NO. 4: MCMC Bayesian Analysis ITERATIONTHETA1 THETA2 THETA3THETA4 SIGMA(1,1) OMEGA(1,1) OMEGA(2,2) OBJ -1 1.63523E+00 1.56116E+00 7.51601E-01 2.35158E+00 5.71097E-02 1.66941E-01 1.39843E-01 -2573 - 1.60770E+00 1.48763E+00 7.25607E-01 2.41005E+00 4.15829E-02 1.75023E-01 1.14078E-01 -2588 -9998 1.67015E+00 1.50197E+00 8.04380E-01 2.32958E+00 4.60430E-02 1.68910E-01 1.70382E-01 -2548 -9997 1.60714E+00 1.56161E+00 7.36944E-01 2.37716E+00 4.96144E-02 1.35797E-01 1.62153E-01 -2539 Thanks Santosh [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] List of tags in roxygen and use for S4 classes?
On Fri, Sep 4, 2009 at 4:47 PM, Bernd Bischlbernd_bis...@gmx.net wrote: Rainer M Krug wrote: Hi is there a list of all roxygen tags which are available? I couldn't find them. I am asking specifically towards the use of roxygen in documenting S4 classes - is that implemented yet (i am using roxygen 0.1 from CRAN at the moment)? Thanks Rainer I am using it do document a now rather large S4 package - with the new development version on rforge (its 0.1-1). It basically works, _but_ : - It took me quite some time to write everything the way roxygen wants it - Error messages are often not very informative. - Many (legal) ways to specify S4 classes / signatures will produce errors, you need to be very verbose. Thanks a lot for the info. - Even disregarding there are still some clear bugs left IMHO. I would also wish for a better (online) documentation, as I think the general idea of roxygen is great. I agree completely. . So it's worth a try but expect some obstacles. There's also a roxygen devel list, where you could go with questions. Good idea. I will post there with questions. I got some S4 examples from one of the recent roxygen talks, If you want them I could forward those. Yes - that would be great. It would give me an idea abut the tags. Cheers, Rainer Bernd -- Rainer M. Krug, Centre of Excellence for Invasion Biology, Stellenbosch University, South Africa __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] linear mixed model question
Dear Wen, Since each worker only works on one machine, your model fm2 does not make sense. Your random effects tries to model how the effect of each worker differs between machines. But you don't have that kind of information if a work only works one machine. HTH, Thierry ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest Cel biometrie, methodologie en kwaliteitszorg / Section biometrics, methodology and quality assurance Gaverstraat 4 9500 Geraardsbergen Belgium tel. + 32 54/436 185 thierry.onkel...@inbo.be www.inbo.be To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher The plural of anecdote is not data. ~ Roger Brinner The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey -Oorspronkelijk bericht- Van: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] Namens Wen Huang Verzonden: zondag 6 september 2009 17:49 Aan: r-help@r-project.org Onderwerp: [R] linear mixed model question Hello, I wanted to fit a linear mixed model to a data that is similar in terms of design to the 'Machines' data in 'nlme' package except that each worker (with triplicates) only operates one machine. I created a subset of observations from 'Machines' data such that it looks the same as the data I wanted to fit the model with (see code below). I fitted a model in which 'Machine' was a fixed effect and 'Worker' was random (intercept), which ran perfectly. Then I decided to complicate the model a little bit by fitting 'Worker' within 'Machine', which was saying variation among workers was nested within each machine. The model could be fitted by 'lme', but when I tried to get confidence intervals by 'intervals(fm2)' it gave me an error: Error in intervals.lme(fm2) : Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance I am wondering if this is because it is impossible to fit a model like 'fm2' or there is some other reasons? Thanks a lot! Wen # library(nlme) data(Machines) new.data = Machines[c(1:6, 25:30, 49:54), ] fm1 = lme(score ~ Machine, random = ~1|Worker, data = new.data) fm1 fm2 = lme(score ~ Machine, random = ~Machine-1|Worker, data = new.data) fm2 intervals(fm2) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Druk dit bericht a.u.b. niet onnodig af. Please do not print this message unnecessarily. Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document. The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Using anova(f1, f2) to compare lmer models yields seemingly erroneous Chisq = 0, p = 1
Are the two models (f1 and f2) actually nested? Aside from that, it is strange that the output is exactly the same after you used REML=FALSE. The log likelihoods should have changed. Best, -- Wolfgang Viechtbauer Department of Methodology and Statistics School for Public Health and Primary Care University of Maastricht, The Netherlands http://www.wvbauer.com/ Original Message From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Matt Killingsworth Sent: Friday, September 04, 2009 22:29 To: Bert Gunter Cc: r-help@r-project.org Subject: Re: [R] Using anova(f1, f2) to compare lmer models yields seemingly erroneous Chisq = 0, p = 1 Hi Bert, Thank you for your note! I tried changing the REML default, and it still produces the same result (see below). Is that what you meant for me to try? Incidentally, I am using lmer() not lme() ### ORIGINAL ### f1 - (lmer(outcome ~ predictor.1 + (1 | person), data=i)) f2 - (lmer(outcome ~ predictor.2 + (1 | person), data=i)) anova(f1, f2) Data: i Models: f1: outcome ~ predictor.1 + (1 | person) f2: outcome ~ predictor.2 + (1 | person) DfAIC BIClogLik Chisq Chi Df Pr(Chisq) f1 6 45443 45489 -22715 f2 25 47317 47511 -23633 0 19 1 ### DO NOT USE REML ### f1 - (lmer(outcome ~ predictor.1 + (1 | person), data=i, REML = FALSE)) f2 - (lmer(outcome ~ predictor.2 + (1 | person), data=i, REML = FALSE)) anova(f1, f2) Data: i Models: f1: outcome ~ predictor.1 + (1 | person) f2: outcome ~ predictor.2 + (1 | person) DfAIC BIClogLik Chisq Chi Df Pr(Chisq) f1 6 45443 45489 -22715 f2 25 47317 47511 -23633 0 19 1 On Fri, Sep 4, 2009 at 4:18 PM, Bert Gunter gunter.ber...@gene.com wrote: My guess would be: Likelihood comparisons are not meaningful for objects fit using restricted maximum likelihood and with different fixed effects. (from ?anova.lme in the nlme package). Are you using the REML = TRUE default? Bert Gunter Genentech Nonclinical Statistics -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of rapton Sent: Friday, September 04, 2009 9:10 AM To: r-help@r-project.org Subject: [R] Using anova(f1, f2) to compare lmer models yields seemingly erroneous Chisq = 0, p = 1 Hello, I am using R to analyze a large multilevel data set, using lmer() to model my data, and using anova() to compare the fit of various models. When I run two models, the output of each model is generated correctly as far as I can tell (e.g. summary(f1) and summary(f2) for the multilevel model output look perfectly reasonable), and in this case (see below) predictor.1 explains vastly more variance in outcome than predictor.2 (R2 = 15% vs. 5% in OLS regression, with very large N). What I am utterly puzzled by is that when I run an anova comparing the two multilevel model fits, the Chisq comes back as 0, with p = 1. I am pretty sure that fit #1 (f1) is a much better predictor of the outcome than f2, which is reflected in the AIC, BIC , and logLik values. Why might anova be giving me this curious output? How can I fix it? I am sure I am making a dumb error somewhere, but I cannot figure out what it is. Any help or suggestions would be greatly appreciated! -Matt f1 - (lmer(outcome ~ predictor.1 + (1 | person), data=i)) f2 - (lmer(outcome ~ predictor.2 + (1 | person), data=i)) anova(f1, f2) Data: i Models: f1: outcome ~ predictor.1 + (1 | person) f2: outcome ~ predictor.2 + (1 | person) DfAIC BIClogLik Chisq Chi Df Pr(Chisq) f1 6 45443 45489 -22715 f2 25 47317 47511 -23633 0 19 1 -- __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How Does One Use the Value of Density Function?
On 07-Sep-09 06:42:06, Gundala Viswanath wrote: How do people usually use the result of density function (e.g. dnorm)? Especially when its value can be greater than 1. What do they do with such density 1? dnorm(2.02,2,.24) [1] 1.656498 - G.V. The point is that it is a *density* of probability. The greater the density, the more compactly a given amount of probability is distributed over a given range of X (or, the narrower the range of X ove which a given amount of probability is distributed); the lower the density, the more widely it is dispersed. The concept of density is, in effect, the amount of probability per unit length in an interval, so division by length is part of it. To convert probability density to probability, multiply by a length. To obtain (a close approximation to) the amount of probability within a short range, such as 2.01 to 2.03, when X has your Normal distribution with mean 2.02 and SD 0.24, multiply the value of the density function at the midpoint 2.02 by the length 0.02 of the interval: dnorm(2.02,2,.24)*0.02 # = 0.03312996 and compare it with the amount of probability calculated from pnorm: pnorm(2.03,2,.24) - pnorm(2.01,2,.24) # = 0.03312044 The approximation is even closer for shorter intervals. Consider the ratio between the approximate and the true probabilities for an interval of length 0.0002: (dnorm(2.02,2,.24)*0.0002)/(pnorm(2.0201,2,.24) - pnorm(2.0199,2,.24)) # = 1 In fact, 1 - (dnorm(2.02,2,.24)*0.0002)/ (pnorm(2.0201,2,.24) - pnorm(2.0199,2,.24)) # = -2.873419e-08 Ted. E-Mail: (Ted Harding) ted.hard...@manchester.ac.uk Fax-to-email: +44 (0)870 094 0861 Date: 07-Sep-09 Time: 10:14:00 -- XFMail -- __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] getting name from function?
Hi
I have the following function which should return the name of FUN:
myFUN - function(FUN) {
return( THE_NAME_OF_FUN(FUN))
}
Is it possible? What do I have tio use here instead of THE_NAME_OF_FUN?
Thanks,
Rainer
--
Rainer M. Krug, Centre of Excellence for Invasion Biology, Stellenbosch
University, South Africa
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] getting name from function?
Hi,
Try this,
myFUN - function(FUN) {
return(deparse(substitute(FUN)))
}
HTH,
baptiste
2009/9/7 Rainer M Krug r.m.k...@gmail.com
Hi
I have the following function which should return the name of FUN:
myFUN - function(FUN) {
return( THE_NAME_OF_FUN(FUN))
}
Is it possible? What do I have tio use here instead of THE_NAME_OF_FUN?
Thanks,
Rainer
--
Rainer M. Krug, Centre of Excellence for Invasion Biology, Stellenbosch
University, South Africa
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
_
Baptiste Auguié
School of Physics
University of Exeter
Stocker Road,
Exeter, Devon,
EX4 4QL, UK
http://newton.ex.ac.uk/research/emag
__
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] getting name from function?
On Mon, Sep 7, 2009 at 11:26 AM, baptiste auguie
baptiste.aug...@googlemail.com wrote:
Hi,
Try this,
myFUN - function(FUN) {
return(deparse(substitute(FUN)))
}
Thanks - that's it
Rainer
HTH,
baptiste
2009/9/7 Rainer M Krug r.m.k...@gmail.com
Hi
I have the following function which should return the name of FUN:
myFUN - function(FUN) {
return( THE_NAME_OF_FUN(FUN))
}
Is it possible? What do I have tio use here instead of THE_NAME_OF_FUN?
Thanks,
Rainer
--
Rainer M. Krug, Centre of Excellence for Invasion Biology, Stellenbosch
University, South Africa
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
_
Baptiste Auguié
School of Physics
University of Exeter
Stocker Road,
Exeter, Devon,
EX4 4QL, UK
http://newton.ex.ac.uk/research/emag
__
--
Rainer M. Krug, Centre of Excellence for Invasion Biology, Stellenbosch
University, South Africa
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] OT - Banker's Algorithum
Hi folks, Can R-Project be used to perform Banker's Algorithum? http://en.wikipedia.org/wiki/Banker%27s_algorithm Pointer would be appreciated. TIA B.R. Stephen L [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Rmetrics: Problem with align
Hi there! I'm stuck with a problem aligning financial timeseries and haven't found a cue how to fix it... When I run that simple script, everything goes well until the align-command: -- rm(list=ls()) x - yahooSeries(^GDAXI) head(x) xAligned - align(x = x, by = 1d, method = before, include.weekends = FALSE) -- Here's the command-line dialog: -- rm(list=ls()) x - yahooSeries(^GDAXI) versuche URL 'http://chart.yahoo.com/table.csv?s=^GDAXIa=b=c=d=8e=07f=2009g=dx=.csv' Content type 'text/csv' length unknown URL geöffnet .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. downloaded 256 Kb head(x) GMT ^GDAXI.Open ^GDAXI.High ^GDAXI.Low ^GDAXI.Close ^GDAXI.Volume 2009-09-04 5328.59 5404.305320.28 5384.43 28535100 2009-09-03 5332.50 5364.035278.50 5301.42 25969000 2009-09-02 5325.94 5337.845263.11 5319.84 30401300 2009-09-01 5479.35 5507.485327.29 5327.29 31245400 2009-08-31 5474.62 5495.565437.92 5458.04 10442600 2009-08-28 5512.17 5573.895502.46 5517.35 23480600 ^GDAXI.Adj.Close 2009-09-04 5384.43 2009-09-03 5301.42 2009-09-02 5319.84 2009-09-01 5327.29 2009-08-31 5458.04 2009-08-28 5517.35 xAligned - align(x = x, by = 1d, method = before, include.weekends = FALSE) Fehler in seq.default(u[1] + offset, u[length(u)], by = by) : wrong sign in 'by' argument - Besides, Example-Data (SWX) from the fPortfolio-package is aligned nicely by exactly the same command. Does anybody have an idea what the problem could be? I'm grateful for any help! Kind regards, Gero __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Using anova(f1, f2) to compare lmer models yields seemingly erroneous Chisq = 0, p = 1
On Mon, 2009-09-07 at 10:23 +0200, Viechtbauer Wolfgang (STAT) wrote: Are the two models (f1 and f2) actually nested? Aside from that, it is strange that the output is exactly the same after you used REML=FALSE. The log likelihoods should have changed. I might be completely misremembering, but I recall a thread in R-SIG-Mixed where this was discussed and it was pointed out that anova(...) on mer objects extracts the ML information even if fitted using REML = TRUE. The log likelihoods of the models supplied to 'anova' are being extracted using REML = FALSE. So, if the above is correct, it does not surprise me that there is no difference. 'anova' was doing the right thing in both cases. See ?mer-class for more details, then try: logLik(f1, REML = FALSE) logLik(f1, REML = TRUE) logLik(f2, REML = FALSE) logLik(f2, REML = TRUE) 'anova' is calling logLik with REML = FALSE regardless of what you define in your model fitting call. HTH G Best, -- Wolfgang Viechtbauer Department of Methodology and Statistics School for Public Health and Primary Care University of Maastricht, The Netherlands http://www.wvbauer.com/ Original Message From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Matt Killingsworth Sent: Friday, September 04, 2009 22:29 To: Bert Gunter Cc: r-help@r-project.org Subject: Re: [R] Using anova(f1, f2) to compare lmer models yields seemingly erroneous Chisq = 0, p = 1 Hi Bert, Thank you for your note! I tried changing the REML default, and it still produces the same result (see below). Is that what you meant for me to try? Incidentally, I am using lmer() not lme() ### ORIGINAL ### f1 - (lmer(outcome ~ predictor.1 + (1 | person), data=i)) f2 - (lmer(outcome ~ predictor.2 + (1 | person), data=i)) anova(f1, f2) Data: i Models: f1: outcome ~ predictor.1 + (1 | person) f2: outcome ~ predictor.2 + (1 | person) DfAIC BIClogLik Chisq Chi Df Pr(Chisq) f1 6 45443 45489 -22715 f2 25 47317 47511 -23633 0 19 1 ### DO NOT USE REML ### f1 - (lmer(outcome ~ predictor.1 + (1 | person), data=i, REML = FALSE)) f2 - (lmer(outcome ~ predictor.2 + (1 | person), data=i, REML = FALSE)) anova(f1, f2) Data: i Models: f1: outcome ~ predictor.1 + (1 | person) f2: outcome ~ predictor.2 + (1 | person) DfAIC BIClogLik Chisq Chi Df Pr(Chisq) f1 6 45443 45489 -22715 f2 25 47317 47511 -23633 0 19 1 On Fri, Sep 4, 2009 at 4:18 PM, Bert Gunter gunter.ber...@gene.com wrote: My guess would be: Likelihood comparisons are not meaningful for objects fit using restricted maximum likelihood and with different fixed effects. (from ?anova.lme in the nlme package). Are you using the REML = TRUE default? Bert Gunter Genentech Nonclinical Statistics -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of rapton Sent: Friday, September 04, 2009 9:10 AM To: r-help@r-project.org Subject: [R] Using anova(f1, f2) to compare lmer models yields seemingly erroneous Chisq = 0, p = 1 Hello, I am using R to analyze a large multilevel data set, using lmer() to model my data, and using anova() to compare the fit of various models. When I run two models, the output of each model is generated correctly as far as I can tell (e.g. summary(f1) and summary(f2) for the multilevel model output look perfectly reasonable), and in this case (see below) predictor.1 explains vastly more variance in outcome than predictor.2 (R2 = 15% vs. 5% in OLS regression, with very large N). What I am utterly puzzled by is that when I run an anova comparing the two multilevel model fits, the Chisq comes back as 0, with p = 1. I am pretty sure that fit #1 (f1) is a much better predictor of the outcome than f2, which is reflected in the AIC, BIC , and logLik values. Why might anova be giving me this curious output? How can I fix it? I am sure I am making a dumb error somewhere, but I cannot figure out what it is. Any help or suggestions would be greatly appreciated! -Matt f1 - (lmer(outcome ~ predictor.1 + (1 | person), data=i)) f2 - (lmer(outcome ~ predictor.2 + (1 | person), data=i)) anova(f1, f2) Data: i Models: f1: outcome ~ predictor.1 + (1 | person) f2: outcome ~ predictor.2 + (1 | person) DfAIC BIClogLik Chisq Chi Df Pr(Chisq) f1 6 45443 45489 -22715 f2 25 47317 47511 -23633 0 19 1 -- __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% Dr.
Re: [R] Color index in image function
On Sat, 2009-09-05 at 04:14 -0700, FMH wrote: Dear All, I was looking for the color index in image function, such as from topo.colors(n) and etc. but still never found it. For instance, from the help menu. ### # Volcano data visualized as matrix. Need to transpose and flip # matrix horizontally. image(t(volcano)[ncol(volcano):1,]) # A prettier display of the volcano x - 10*(1:nrow(volcano)) y - 10*(1:ncol(volcano)) image(x, y, volcano, col = terrain.colors(100), axes = FALSE) contour(x, y, volcano, levels = seq(90, 200, by = 5), add = TRUE, col = peru) axis(1, at = seq(100, 800, by = 100)) axis(2, at = seq(100, 600, by = 100)) box() title(main = Maunga Whau Volcano, font.main = 4) # From the script above, it yields a beautiful image of volcano with variety of colors but i have to list down the color index that could show the meaning of each color in my thesis. Could someone please help me to extract this color index? Thank you Fir If I understand your question you need change the Palette of image plot. So you need use colorRampPalette look my example Brazilan.Pallete - colorRampPalette(c(green,yellow, blue)) image(x, y, volcano, col = Brazilan.Pallete(50), axes = FALSE) -- Bernardo Rangel Tura, M.D,MPH,Ph.D National Institute of Cardiology Brazil __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Matrix regression
Dear friends, I would like to solve the following regression problem: y=c1 x1 + c2 x2 + + cn xn where the y, xi are all matrices and the ci are constants that need to be determined. The y, xi are distance matrices (symmetric). ci should be forced to positive or null (i.e. non negative). Any suggestion? I will be more than happy to share the results of my quest with the list or developers. Regards -- Corrado Topi Global Climate Change Biodiversity Indicators Area 18,Department of Biology University of York, York, YO10 5YW, UK Phone: + 44 (0) 1904 328645, E-mail: ct...@york.ac.uk __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Fwd: trouble installing gtools package in local directory
Hello, I posted my question to r-sig-fedora a few days ago, but noticed that list is not very active. I hope sending my question here is okay, so here goes. I've recently transitioned to using R in Linux. My OS/installation versions are: $ Red Hat Enterprise Linux AS release 4 (Nahant Update 6) $ Linux lx-chgmqsd05 2.6.9-67.0.1.ELsmp #1 SMP Fri Nov 30 11:57:43 EST 2007 x86_64 x86_64 x86_64 GNU/Linux I am trying to create a repository of packages under ~/R/library which I've created since I do not have admin access on my system. From my research on the web, this is the standard way of installing and using packages from CRAN in linux. I am ultimately trying to install the package gplots which depends on gtools. I have tried to install gtools using both of the 2 ways I know: i) Download gtools_2.6.1.tar.gz from CRAN and then running R CMD INSTALL gtools_2.6.1.tar.gz ii) install.packages(gtools,lib=~/R/library) directly at the R prompt I get the same result with both and not sure what I'm missing here. I have a slight suspicion that I don't have R_LIBS or R_HOME set correctly but I've messed around with those, to no avail. I'm including the result of attempting the install below. Any help? $ R CMD INSTALL gtools_2.6.1.tar.gz WARNING: ignoring environment value of R_HOME * Installing to library '/home/ssaeed/R/library/' * Installing *source* package 'gtools' ... ** libs gcc -std=gnu99 -I/tp64/r-project/2.9.0/lib64/R/include -I/usr/local/include -fpic -O2 -m64 -march=nocona -mmmx -msse -msse2 -msse3 -c setTCPNoDelay.c -o setTCPNoDelay.o gcc -std=gnu99 -shared -L/usr/local/lib64 -o gtools.so setTCPNoDelay.o -L/tp64/r-project/2.9.0/lib64/R/lib -lR ** R ** data ** inst ** preparing package for lazy loading ** help *** installing help indices Building/Updating help pages for package 'gtools' Formats: text html latex example Can't use an undefined value as filehandle reference at /tp64/r-project/2.9.0/lib64/R/share/perl/R/Rdconv.pm line 81. ERROR: building help failed for package 'gtools' * Removing '/home/ssaeed/R/library//gtools' I'm returning to this group after a long time so please let me know if I'm posting at the wrong place. Thanks, Saad [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How many attributes are there of a variable?
Peng, based on a suggestion, Frank made years ago (18.7.2006), I use one attribute that contains all further attributes, I want to assign to variables. It's necessary to create your own class and subsetting method, so that this attribute does not get lost. Together with some functions I use labels for variables, value.labels, missing.value definitions etc. It seems, without protection by your own class and the corresponding subsetting method, you can never be sure, if an attribute survives subsetting. Heinz At 23:21 06.09.2009, Frank E Harrell Jr wrote: Peng, You can create all the attributes you want, with one headache: R does not keep attributes across subsetting operations so you need to write classes and [.something methods when attributions need to be kept or adjusted upon subsetting rows. The Hmisc package uses attributes such as label, units, imputed. You might look at the code to see how it did that. For example, label(x) will use attr(x, 'label') to fetch the 'label' attribute. There are attribute-setting functions there too. Frank Peng Yu wrote: Hi, According to the example below this email, attr(x,names) is the same as names(x). I am wondering how many attributes there are of a given variable. How to find out what they are? Can I always use some_attribute(x) instead of attr(x, some_attribute)? Regards, Peng x=c(1,2,3) attr(x,names)=c(a,b,c) x a b c 1 2 3 y=c(1,2,3) names(y)=c(a,b,c) y a b c 1 2 3 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Frank E Harrell Jr Professor and Chair School of Medicine Department of Biostatistics Vanderbilt University __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problems with Boxplot
Hi r-help-boun...@r-project.org napsal dne 05.09.2009 04:59:41: Hi Petr, Thanks for these comments. I'm sorry that my post was not clear. I was referring to the questions in my original post/code/file uploads, but I had forgotten to include an updated file (now attached http://www.nabble.com/file/p25304663/Post%2Btrial%2Bdata.csv Post+trial+data.csv ) to work with the new code: testdata- c(C:\\Files\\R\\Sample R code\\Post trial data.csv) new_data- read.table(testdata, skip = 0, sep = ,, na.strings = na,header = TRUE) x11(width=16, height=7, pointsize=14) boxplot(new_data,outline = FALSE, col = c(lightblue, salmon), las =1, boxwex = 0.5) legend(top, c(Label for blue boxes,Label for red boxes), cex=1.5, lty=1:2, fill=c(lightblue, salmon), bty=n); title(main=Chart title text, cex.main = 1.8) grid() I'm still not clear how I can get the number format showing #,###. E.g. with this code and attached file, the scale shows as 2000, 1 etc. I don't know how to show 2,000. 10,000 etc. I have looked through sprintf (thanks for suggesting that - I'd spent hours looking without finding it) and it seems incredibly flexible, but the formats shown are more scientific in focus. I still haven't been able to find a way of getting a comma style. AFAIK you can not format these in boxplot directly. You need to plot without y axis and in axis you can use formating with prettyNum. I found quite easily from sprintf and formatC help pages (I did not do it before so I learned it now:-) x-rnorm(100)+1 bbb-boxplot(x, axes=F) axis(2, at= pretty(x), labels=prettyNum(pretty(x), big.mark=,)) Regards Petr Thanks again Guy Petr Pikal wrote: Hi it is rather difficult to understand what you mean by your questions/answers without real reproducible code. r-help-boun...@r-project.org napsal dne 03.09.2009 13:41:11: I'd be interested if anyone has a quick way to get percentages and additionally, how do I get numbers in the 0,000 format along the x or y-axis? In the meantime, I can live with this. plot(1:10,1:10, axes=F) axis(2, at=c(2,3,7,9), labels=c(1.2, 2.38, 13.54, 16.8)) the same applies with boxplot. by bbb- boxplot() you obtain an object which is used by bxp. See help page for boxplot, section See also ... See also par for graphic options, format and or sprintf for formating numbers Regards Petr -- View this message in context: http://www.nabble.com/Problems-with-Boxplot- tp25256461p25304663.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] [Hmisc] Latex to pdf
Hi Jeroen,
From: Jeroen Ooms j.c.l.o...@uu.nl
To: r-help@r-project.org
Date: Sun, 6 Sep 2009 04:49:08 -0700 (PDT)
Subject: [R] [Hmisc] Latex to pdf
I would like to print some tables and figures to a PDF device on a CentOS 5
vps. However, I cannot seem to get the latex function from Hmisc working. I
followed the example, and got an error: sh: xdvi: command not found. I tried
installing the 'tetex-xdvi' linux package, and now it returns: Error: Can't
open display. I guess the reason for this is that the machine is a VPS
terminal, so it has no display, however I dont understand why it does not
write to the PDF device. Some code:
Yep, the VSP terminal is likely causing the problem.
pdf(test.pdf)
No, don't call latex() inside a pdf device like this. The latex()
command produces a .tex file, which can then be processed by the
system latex or pdflatex etc.
x - matrix(1:6, nrow=2, dimnames=list(c('a','b'),c('c','d','this that')))
latex(x) # creates x.tex in working directory
And what happens when you call latex x.tex from the system console
or system(latex x.tex) from R? Does it work? You can avoid having
the latex() function try to create and display the dvi by setting the
file name: latex(x, file=abcThisThat.tex).
Hope it helps,
Ista
snip
--
Ista Zahn
Graduate student
University of Rochester
Department of Clinical and Social Psychology
http://yourpsyche.org
Jeroen Ooms * Dept. of Methodology and Statistics * Utrecht University
Visit http://www.jeroenooms.com www.jeroenooms.com to explore some of my
current projects.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Andrews plot
Dear all Colleague of mine ask me if R is capable of Andrews plot like andrewsplot(x) in Matlab. Quick search did not reveal anything but before I start to write any routine I would like to ask this ingenious audience if there is any implementation of Andrews plots somewhere. I know about parallel coordinate plots in lattice (although I do not use them as I am not sure what the plot tells me :-). Regards Petr __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Andrews plot
On 09/07/2009 03:36 PM, Petr PIKAL wrote: Dear all Colleague of mine ask me if R is capable of Andrews plot like andrewsplot(x) in Matlab. Quick search did not reveal anything but before I start to write any routine I would like to ask this ingenious audience if there is any implementation of Andrews plots somewhere. I know about parallel coordinate plots in lattice (although I do not use them as I am not sure what the plot tells me :-). Regards Petr Hi, The example I find on the matlab help page (http://tr.im/y5mf) looks like this : http://addictedtor.free.fr/graphiques/search.php?q=andrewengine=RGG Romain -- Romain Francois Professional R Enthusiast +33(0) 6 28 91 30 30 http://romainfrancois.blog.free.fr |- http://tr.im/xMdt : update on the ant package |- http://tr.im/xHLs : R capable version of ant `- http://tr.im/xHiZ : Tip: get java home from R with rJava __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Andrews plot
On Mon, Sep 7, 2009 at 8:36 AM, Petr PIKALpetr.pi...@precheza.cz wrote:
Dear all
Colleague of mine ask me if R is capable of Andrews plot like
andrewsplot(x) in Matlab.
Quick search did not reveal anything but before I start to write any
routine I would like to ask this ingenious audience if there is any
implementation of Andrews plots somewhere.
Here's the code I use in the tourr package:
#' Compute Andrews' curves
#'
#' This function takes a numeric vector of input, and returns a function which
#' allows you to compute the value of the Andrew's curve at every point along
#' its path from -pi to pi.
#'
#' @param x input a new parameter
#' @return a function with single argument, theta
#'
#' @examples
#' a - andrews(1:2)
#' a(0)
#' a(-pi)
#' grid - seq(-pi, pi, length = 50)
#' a(grid)
#'
#' plot(grid, andrews(1:2)(grid), type = l)
#' plot(grid, andrews(runif(5))(grid), type = l)
andrews - function(x) {
n - length(x)
y - rep(x[1] / sqrt(2), length(t))
function(t) {
for(i in seq(2, n, by = 1)) {
val - i %/% 2 * t
y - y + x[i] * (if(i %% 2 == 0) sin(val) else cos(val))
}
y / n
}
}
I know about parallel coordinate plots in lattice (although I do not use
them as I am not sure what the plot tells me :-).
If you don't understand parallel coordinates plots, I think you're
going to find Andrew's curves even harder to understand.
Hadley
--
http://had.co.nz/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Andrews plot
Thank you.
hadley wickham h.wick...@gmail.com napsal dne 07.09.2009 15:50:03:
On Mon, Sep 7, 2009 at 8:36 AM, Petr PIKALpetr.pi...@precheza.cz
wrote:
Dear all
Colleague of mine ask me if R is capable of Andrews plot like
andrewsplot(x) in Matlab.
Quick search did not reveal anything but before I start to write any
routine I would like to ask this ingenious audience if there is any
implementation of Andrews plots somewhere.
Here's the code I use in the tourr package:
#' Compute Andrews' curves
#'
#' This function takes a numeric vector of input, and returns a function
which
#' allows you to compute the value of the Andrew's curve at every point
along
#' its path from -pi to pi.
#'
#' @param x input a new parameter
#' @return a function with single argument, theta
#'
#' @examples
#' a - andrews(1:2)
#' a(0)
#' a(-pi)
#' grid - seq(-pi, pi, length = 50)
#' a(grid)
#'
#' plot(grid, andrews(1:2)(grid), type = l)
#' plot(grid, andrews(runif(5))(grid), type = l)
andrews - function(x) {
n - length(x)
y - rep(x[1] / sqrt(2), length(t))
function(t) {
for(i in seq(2, n, by = 1)) {
val - i %/% 2 * t
y - y + x[i] * (if(i %% 2 == 0) sin(val) else cos(val))
}
y / n
}
}
I know about parallel coordinate plots in lattice (although I do not
use
them as I am not sure what the plot tells me :-).
If you don't understand parallel coordinates plots, I think you're
going to find Andrew's curves even harder to understand.
As I said, the question was from my colleague. Maybe he knows what to
expect from such plots. I just told him that I did not find any direct
implementation but that it seems to me rather straightforward to make such
a plot if I knew what shall be on x axis and on y axis.
Regards
Petr
Hadley
--
http://had.co.nz/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] R-crash when loading workspace - Windows
Dear all, One day when I tried to load an existing workspace (when opening R or by load()), R crashed without any error notification. The day before I had worked and saved my workspace without any trouble. At first I though it was a memory problem (workspace reaching 180Mo) or related to a particular script or command, so I start a new workspace. Everything was ok, that script and others working. Then I saved the workspace (55Mo) and tried to open it, without any result : R crashes without any notification again. This occurs only with Windows. Does someone know how to solve that problem? Regards, Edwige. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Color index in image function
Thank you for the tips. I have manage to run your script, but was still never get the way to include the color index beside the image which could explain the intensity of the color from the lower index(green) to the higher index(blue). This color index might be represented by an increasing of color index in another table beside the image, started from green followed by green-yellow, yellow, yellow-blue and blue? Could someone please advice on this matter? Cheers Fir - Original Message From: Bernardo Rangel Tura t...@centroin.com.br To: FMH kagba2...@yahoo.com Sent: Monday, September 7, 2009 11:13:12 AM Subject: Re: [R] Color index in image function On Sat, 2009-09-05 at 04:14 -0700, FMH wrote: Dear All, I was looking for the color index in image function, such as from topo.colors(n) and etc. but still never found it. For instance, from the help menu. ### # Volcano data visualized as matrix. Need to transpose and flip # matrix horizontally. image(t(volcano)[ncol(volcano):1,]) # A prettier display of the volcano x - 10*(1:nrow(volcano)) y - 10*(1:ncol(volcano)) image(x, y, volcano, col = terrain.colors(100), axes = FALSE) contour(x, y, volcano, levels = seq(90, 200, by = 5), add = TRUE, col = peru) axis(1, at = seq(100, 800, by = 100)) axis(2, at = seq(100, 600, by = 100)) box() title(main = Maunga Whau Volcano, font.main = 4) # From the script above, it yields a beautiful image of volcano with variety of colors but i have to list down the color index that could show the meaning of each color in my thesis. Could someone please help me to extract this color index? Thank you Fir If I understand your question you need change the Palette of image plot. So you need use colorRampPalette look my example Brazilan.Pallete - colorRampPalette(c(green,yellow, blue)) image(x, y, volcano, col = Brazilan.Pallete(50), axes = FALSE) -- Bernardo Rangel Tura, M.D,MPH,Ph.D National Institute of Cardiology Brazil __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Rmetrics: Problem with align
Hi all! I've found the solution - the ordering of the series is important. It's expected that the most recent entries are located at the bottom of the data matrix. Here's the sorting command: x - sort(x, decreasing = FALSE) Kind regards: Gero Gero Schwenk schrieb: Hi there! I'm stuck with a problem aligning financial timeseries and haven't found a cue how to fix it... When I run that simple script, everything goes well until the align-command: -- rm(list=ls()) x - yahooSeries(^GDAXI) head(x) xAligned - align(x = x, by = 1d, method = before, include.weekends = FALSE) -- Here's the command-line dialog: -- rm(list=ls()) x - yahooSeries(^GDAXI) versuche URL 'http://chart.yahoo.com/table.csv?s=^GDAXIa=b=c=d=8e=07f=2009g=dx=.csv' Content type 'text/csv' length unknown URL geöffnet .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. downloaded 256 Kb head(x) GMT ^GDAXI.Open ^GDAXI.High ^GDAXI.Low ^GDAXI.Close ^GDAXI.Volume 2009-09-04 5328.59 5404.305320.28 5384.43 28535100 2009-09-03 5332.50 5364.035278.50 5301.42 25969000 2009-09-02 5325.94 5337.845263.11 5319.84 30401300 2009-09-01 5479.35 5507.485327.29 5327.29 31245400 2009-08-31 5474.62 5495.565437.92 5458.04 10442600 2009-08-28 5512.17 5573.895502.46 5517.35 23480600 ^GDAXI.Adj.Close 2009-09-04 5384.43 2009-09-03 5301.42 2009-09-02 5319.84 2009-09-01 5327.29 2009-08-31 5458.04 2009-08-28 5517.35 xAligned - align(x = x, by = 1d, method = before, include.weekends = FALSE) Fehler in seq.default(u[1] + offset, u[length(u)], by = by) : wrong sign in 'by' argument - Besides, Example-Data (SWX) from the fPortfolio-package is aligned nicely by exactly the same command. Does anybody have an idea what the problem could be? I'm grateful for any help! Kind regards, Gero __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] linear mixed model question
Hello Wen: On 2009.09.06 10:49:03, Wen Huang wrote: Hello, I wanted to fit a linear mixed model to a data that is similar in terms of design to the 'Machines' data in 'nlme' package except that each worker (with triplicates) only operates one machine. I created a subset of observations from 'Machines' data such that it looks the same as the data I wanted to fit the model with (see code below). I fitted a model in which 'Machine' was a fixed effect and 'Worker' was random (intercept), which ran perfectly. Then I decided to complicate the model a little bit by fitting 'Worker' within 'Machine', which was saying variation among workers was nested within each machine. The model could be fitted by 'lme', but when I tried to get confidence intervals by 'intervals(fm2)' it gave me an error: Error in intervals.lme(fm2) : Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance I am wondering if this is because it is impossible to fit a model like 'fm2' or there is some other reasons? The problem doesn't seem to be the model specification but is most likely the result of estimating a more complicated model with very little data. Using the complete Machines dataset with the same model specification seems to work fine: # - # fm3 - lme(score ~ Machine, random = ~ Machine - 1 | Worker, data = Machines) intervals(fm3) Approximate 95% confidence intervals Fixed effects: lower est.upper (Intercept) 48.972459 52.36 55.73865 MachineB 3.093747 7.97 12.83959 MachineC10.816607 13.916667 17.01673 attr(,label) [1] Fixed effects: Random Effects: Level: Worker lower est. upper sd(MachineA)2.1702468 4.0792807 7.6675752 sd(MachineB)4.6301082 8.6252908 16.0677975 sd(MachineC)2.3387870 4.3894795 8.2382579 cor(MachineA,MachineB) 0.1992744 0.8027499 0.9647702 cor(MachineA,MachineC) -0.1702480 0.6225047 0.9260744 cor(MachineB,MachineC) 0.1235115 0.7708309 0.9579666 Within-group standard error: lower est. upper 0.7629124 0.9615766 1.2119736 # - # With the restricted dataset, there are only 18 observations in 6 groups. This is probably too little data for the (restricted) maximum likelihood technique used by lme(). Hope that helps, ~Jason -- Jason W. Morgan Graduate Student Department of Political Science *The Ohio State University* 154 North Oval Mall Columbus, Ohio 43210 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Tinn-R setup
I recently installed R 2.9.2 on a new Windows platform. Everything seemed to installed OK. I then downloaded the latest Tinn-R (2.3.2.3 I think) and as I have always done I selected R - Configure - Permanent. I was greeted with a dialog box asking me for a mirror site. I don't remember this prompt before but I decided to play along and I select a mrror site. Then the process takes off installing what appears to be every package that has ever been conceived for R (translate - alot of packages). Is this normal? I repeated this about three times because at the end it gave me an error indicating sus-and-such library was not found. Each time it was a different library that couldn't be found. Finally on the third try it seemed to complete without error. Again I have never had to go through so many steps to get Tinn-R installed and configured. Did I do something wrong? Is there a bug in this package or a problem with the integration between R 2.9.2 and Tin-R 2.3.2.3? Thank you. Kevin __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Averaging rows if a condition is true.
Dear All,
I have matrix (5 X 60) of subjects and their responses to a set of
questions. All responses are classified into categories (500). I would like to
average all subject's responses for each category. I wrote a code using a for
loop but is not working. Could please tell me what's wrong with the code? I
guess, there is a elegant R way of doing the same thing.
Thanks in advance.
Kind regards,
Ezhil
j - 1; n - dim(dat)[1]; cat - as.character(dat[,1]);
row - matrix(nrow=nrow(dat), ncol=ncol(dat));
for(i in 1:n-1) {
if(cat[i] != cat[i+1]) {row[j, ] - dat[j, ]}
else {
start - j;
end - i;
}
row[j, ] - colMeans(dat[j:i, ]);
j+1;
}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Tinn-R setup
I recently installed R 2.9.2 on a new Windows platform. Everything seemed to installed OK. I then downloaded the latest Tinn-R (2.3.2.3 I think) and as I have always done I selected R - Configure - Permanent. I was greeted with a dialog box asking me for a mirror site. I don't remember this prompt before but I decided to play along and I select a mrror site. Then the process takes off installing what appears to be every package that has ever been conceived for R (translate - alot of packages). Is this normal? I repeated this about three times because at the end it gave me an error indicating sus-and-such library was not found. Each time it was a different library that couldn't be found. Finally on the third try it seemed to complete without error. Again I have never had to go through so many steps to get Tinn-R installed and configured. Did I do something wrong? Is there a bug in this package or a problem with the integration between R 2.9.2 and Tin-R 2.3.2.3? Thank you. Kevin __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Averaging rows if a condition is true.
Hi,
Try to use aggregate function
RSiteSearch (aggregate) #for help
Regards
ML
A Ezhil a écrit :
Dear All,
I have matrix (5 X 60) of subjects and their responses to a set of
questions. All responses are classified into categories (500). I would like to
average all subject's responses for each category. I wrote a code using a for
loop but is not working. Could please tell me what's wrong with the code? I
guess, there is a elegant R way of doing the same thing.
Thanks in advance.
Kind regards,
Ezhil
j - 1; n - dim(dat)[1]; cat - as.character(dat[,1]);
row - matrix(nrow=nrow(dat), ncol=ncol(dat));
for(i in 1:n-1) {
if(cat[i] != cat[i+1]) {row[j, ] - dat[j, ]}
else {
start - j;
end - i;
}
row[j, ] - colMeans(dat[j:i, ]);
j+1;
}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Mohamed Lajnef
INSERM Unité 955.
40 rue de Mesly. 94000 Créteil.
Courriel : mohamed.laj...@inserm.fr
tel. : 01 49 81 31 31 (poste 18470)
Sec : 01 49 81 32 90
fax : 01 49 81 30 99
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] finding the minimum value
Hi all, I'm using a certain procedure to calculate the value of some variable(Bayes risk),B. So I got the values B1, B2, , B1000, each under certain input values and using a long procedure. Now, I want to put the values I got in a nummerical vector and find their minimum value. I think c( ) should work.For example if I have only 10 values I could have used c(B1,B2,B3,B4,B5,B6,B7,B8,B9,B10) But how can I do this for 1000 values? I think the question is really trivial, but I tried to work on it and couldn't reach anythg. Thanks. Maram [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] finding the minimum value
On Mon, 7 Sep 2009, maram salem wrote: Hi all, I'm using a certain? procedure to calculate the value of some variable(Bayes risk),B. So I got the values B1, B2, , B1000, each under certain input values and using a long procedure. Now, I want to put the values I got in a nummerical vector?and find their minimum value. I?think c( ) should work.For example if I have only 10 values I could have used c(B1,B2,B3,B4,B5,B6,B7,B8,B9,B10) But how can I do this for 1000 values? sapply( paste( B, 1:1000, sep='' ), get ) but I wonder why you didn't store the values in a list to start with. HTH, Chuck I think the question is really trivial, but I tried to work on it and couldn't reach anythg. Thanks. Maram [[alternative HTML version deleted]] Charles C. Berry(858) 534-2098 Dept of Family/Preventive Medicine E mailto:cbe...@tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Using anova(f1, f2) to compare lmer models yields seemingly erroneous Chisq = 0, p = 1
rapton wrote: Hello, I am using R to analyze a large multilevel data set, using lmer() to model my data, and using anova() to compare the fit of various models. When I run two models, the output of each model is generated correctly as far as I can tell (e.g. summary(f1) and summary(f2) for the multilevel model output look perfectly reasonable), and in this case (see below) predictor.1 explains vastly more variance in outcome than predictor.2 (R2 = 15% vs. 5% in OLS regression, with very large N). What I am utterly puzzled by is that when I run an anova comparing the two multilevel model fits, the Chisq comes back as 0, with p = 1. I am pretty sure that fit #1 (f1) is a much better predictor of the outcome than f2, which is reflected in the AIC, BIC , and logLik values. Why might anova be giving me this curious output? How can I fix it? I am sure I am making a dumb error somewhere, but I cannot figure out what it is. Any help or suggestions would be greatly appreciated! -Matt f1 - (lmer(outcome ~ predictor.1 + (1 | person), data=i)) f2 - (lmer(outcome ~ predictor.2 + (1 | person), data=i)) anova(f1, f2) Data: i Models: f1: outcome ~ predictor.1 + (1 | person) f2: outcome ~ predictor.2 + (1 | person) DfAIC BIClogLik Chisq Chi Df Pr(Chisq) f1 6 45443 45489 -22715 f2 25 47317 47511 -23633 0 19 1 Your models are nest nestedit doesn't make sense to do. Alain - Dr. Alain F. Zuur First author of: 1. Analysing Ecological Data (2007). Zuur, AF, Ieno, EN and Smith, GM. Springer. 680 p. 2. Mixed effects models and extensions in ecology with R. (2009). Zuur, AF, Ieno, EN, Walker, N, Saveliev, AA, and Smith, GM. Springer. 3. A Beginner's Guide to R (2009). Zuur, AF, Ieno, EN, Meesters, EHWG. Springer Statistical consultancy, courses, data analysis and software Highland Statistics Ltd. 6 Laverock road UK - AB41 6FN Newburgh Email: highs...@highstat.com URL: www.highstat.com -- View this message in context: http://www.nabble.com/Using-anova%28f1%2C-f2%29-to-compare-lmer-models-yields-seemingly-erroneous-Chisq-%3D-0%2C-p-%3D-1-tp25297254p25333120.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Using anova(f1, f2) to compare lmer models yields seemingly erroneous Chisq = 0, p = 1
rapton wrote: Hello, I am using R to analyze a large multilevel data set, using lmer() to model my data, and using anova() to compare the fit of various models. When I run two models, the output of each model is generated correctly as far as I can tell (e.g. summary(f1) and summary(f2) for the multilevel model output look perfectly reasonable), and in this case (see below) predictor.1 explains vastly more variance in outcome than predictor.2 (R2 = 15% vs. 5% in OLS regression, with very large N). What I am utterly puzzled by is that when I run an anova comparing the two multilevel model fits, the Chisq comes back as 0, with p = 1. I am pretty sure that fit #1 (f1) is a much better predictor of the outcome than f2, which is reflected in the AIC, BIC , and logLik values. Why might anova be giving me this curious output? How can I fix it? I am sure I am making a dumb error somewhere, but I cannot figure out what it is. Any help or suggestions would be greatly appreciated! -Matt f1 - (lmer(outcome ~ predictor.1 + (1 | person), data=i)) f2 - (lmer(outcome ~ predictor.2 + (1 | person), data=i)) anova(f1, f2) Data: i Models: f1: outcome ~ predictor.1 + (1 | person) f2: outcome ~ predictor.2 + (1 | person) DfAIC BIClogLik Chisq Chi Df Pr(Chisq) f1 6 45443 45489 -22715 f2 25 47317 47511 -23633 0 19 1 ** NOT ** nested sorrythe brain is going faster than the fingers. - Dr. Alain F. Zuur First author of: 1. Analysing Ecological Data (2007). Zuur, AF, Ieno, EN and Smith, GM. Springer. 680 p. 2. Mixed effects models and extensions in ecology with R. (2009). Zuur, AF, Ieno, EN, Walker, N, Saveliev, AA, and Smith, GM. Springer. 3. A Beginner's Guide to R (2009). Zuur, AF, Ieno, EN, Meesters, EHWG. Springer Statistical consultancy, courses, data analysis and software Highland Statistics Ltd. 6 Laverock road UK - AB41 6FN Newburgh Email: highs...@highstat.com URL: www.highstat.com -- View this message in context: http://www.nabble.com/Using-anova%28f1%2C-f2%29-to-compare-lmer-models-yields-seemingly-erroneous-Chisq-%3D-0%2C-p-%3D-1-tp25297254p25333148.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Creating mixed line and point graphs with xyplot
OOPS Skitts' law ( assuming I'm spelling it correctly. --- On Sun, 9/6/09, David Winsemius dwinsem...@comcast.net wrote: From: David Winsemius dwinsem...@comcast.net Subject: Re: [R] Creating mixed line and point graphs with xyplot To: John Kane jrkrid...@yahoo.ca Cc: Paul Sweeting m...@paulsweeting.co.uk, jim holtman jholt...@gmail.com, r-help@r-project.org Received: Sunday, September 6, 2009, 3:24 PM On Sep 6, 2009, at 3:03 PM, John Kane wrote: I think you have a couple of typos. Should it not be par(new=True) points(x,b) Probably not True --- On Sat, 9/5/09, jim holtman jholt...@gmail.com wrote: From: jim holtman jholt...@gmail.com Subject: Re: [R] Creating mixed line and point graphs with xyplot To: Paul Sweeting m...@paulsweeting.co.uk Cc: r-help@r-project.org Received: Saturday, September 5, 2009, 6:37 PM I would use the base plot routine plot(x,c, type='l', ylim=range(a,c)) points(x,a) park(new=TRUE) plot(x,d,type='l', ylim=range(b,d), axes=FALSE,ylab='', xlab='') pints(x,b) axis(4) On Fri, Sep 4, 2009 at 11:28 AM, Paul Sweetingm...@paulsweeting.co.uk wrote: Hi Well, I think the title says it all! I've looked through the documentation but I can't find a way of doing this. The situation is that I have 4 series, say a, b, c and d. Series a and c are plotted on the lh y axis, series b and d are plotted on the rh (secondary) y axis. I've worked out how to do this. However, I need to plot series a and b a points (symbols only, no line), whislt c and d need plotting as lines (with no symbols). What is the easiest way to do this in xyplot? Or should I be using something else? Thanks! Paul __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Jim Holtman Cincinnati, OH +1 513 646 9390 What is the problem that you are trying to solve? __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ Looking for the perfect gift? Give the gift of Flickr! http://www.flickr.com/gift/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide..html and provide commented, minimal, self-contained, reproducible code. David Winsemius, MD Heritage Laboratories West Hartford, CT __ Looking for the perfect gift? Give the gift of Flickr! http://www.flickr.com/gift/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Averaging rows if a condition is true.
Hi, On Mon, Sep 7, 2009 at 11:46 AM, A Ezhilezhi...@yahoo.com wrote: Dear All, I have matrix (5 X 60) of subjects and their responses to a set of questions. All responses are classified into categories (500). I would like to average all subject's responses for each category. I wrote a code using a for loop but is not working. Could please tell me what's wrong with the code? I guess, there is a elegant R way of doing the same thing. The dlply function in the plyr library provides a very easy and intuitive way to do this ... if you're having problems using it, please post a small representative matrix of your data that we can paste into our own R session to slice and dice for you as an example. -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] calling combinations of variable names
R-2.9.1, Windows7 Dear list, I have a question to you that seems very simple to me, but I just can't figure it out. I have a dataframe called ratings which contains the following variables: evalR1, evalR2, evalR3, evalR4, scoreR1, scoreR2, scoreR3, scoreR4, opinionR1, opinionR2, opinionR3, opinionR4. (there are more variables, but this gives an idea of the data structure). What I want is run several analyses on all 3 or 4-combinations of a given variable. So, for example, I want to compute the following ICC's (function from the psych package): ICC(cbind(evalR1,evalR2, evalR3)) ICC(cbind(evalR1,evalR2, evalR4)) ICC(cbind(evalR1, evalR3, evalR4)) ICC(cbind(evalR2, evalR3, evalR4)) ICC(cbind(evalR1, evalR2, evalR3, eval4)). I create a matrix containing the 3-combinations by combn(4,3). Now I need to call the variables into the function. First, I tried paste as follows: combis - combn(4,3) # this gives the 3-combinations attach(ratings) eval - paste(evalR,combis[1,1],,evalR,combis[2,1],,evalR,combis[3,1],sep =) (this is of course just for 1 combination, as an example) the output of this is evalR1,evalR2,evalR3, but when I run ICC(cbind(eval)), an error message is given which is not given when I enter ICC(cbind(evalR1,evalR2, evalR3)) manually. The function appears not to recognize the variable names. It also does not work to type ICC(cbind(unquote(eval))). Alternatively, I have tried the cat function, but also here ICC does not recognize the input as variable names. What am I doing wrong? How can I automatically construct the set of variable names such that a function recognizes them as variable names? ICC is one example, but there are also other computations to be run and the set of variables is pretty large, so typing the combinations of variable names manually is really unattractive. What am I missing? It seems to me that there probably is a very simple solution in R, but which? Thank you, Peter Verbeet [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] linear mixed model question
Wen, to follow up on Thierry, your workers are nested in machines (since each worker only works one machine). Consider fitting a nested model. Though, with few observations, you might run into the same problem. Further, if you have observation triplets, and you expect systematic differences between each trial, then you would have to include a trial effect in some way. Have you plotted the data? This can be very informative. Finally, you may want to consider a fixed-effects approach. Throw in worker dummies only (without intercept) and test the linear hypothesis that the sum of the worker dummy coefficients is equal between two machines. The following example does that for two machines (if you have n machines you would have binomial coefficient (n,2) linear hypotheses): #simulate data machines=rep(c(0,1),each=9) workers=rep(c(1:6),each=3) workers.re=rep(rnorm(6),each=3) error=rnorm(18,0.3) output=2*machines+workers.re+error #code machines and workers as factors/dummies machines=as.factor(machines) #run OLS (fixed effects) of output on worker dummies without intercept fm4=lm(output~-1+workers) summary(fm4) #test the linear hypothesis that coefficients on the worker dummies are equal for both machines #the equivalent formulation used is: is the sum of the coefficients across machines equal to zero? library(car) linear.hypothesis(fm4,hypothesis.matrix=c(1,1,1,-1,-1,-1),rhs=c(0)) hope that helps, daniel workers=as.factor(workers) Wen Huang-3 wrote: Hello, I wanted to fit a linear mixed model to a data that is similar in terms of design to the 'Machines' data in 'nlme' package except that each worker (with triplicates) only operates one machine. I created a subset of observations from 'Machines' data such that it looks the same as the data I wanted to fit the model with (see code below). I fitted a model in which 'Machine' was a fixed effect and 'Worker' was random (intercept), which ran perfectly. Then I decided to complicate the model a little bit by fitting 'Worker' within 'Machine', which was saying variation among workers was nested within each machine. The model could be fitted by 'lme', but when I tried to get confidence intervals by 'intervals(fm2)' it gave me an error: Error in intervals.lme(fm2) : Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance I am wondering if this is because it is impossible to fit a model like 'fm2' or there is some other reasons? Thanks a lot! Wen # library(nlme) data(Machines) new.data = Machines[c(1:6, 25:30, 49:54), ] fm1 = lme(score ~ Machine, random = ~1|Worker, data = new.data) fm1 fm2 = lme(score ~ Machine, random = ~Machine-1|Worker, data = new.data) fm2 intervals(fm2) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- View this message in context: http://www.nabble.com/linear-mixed-model-question-tp25318961p25333956.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Writing R Scripts and passing command line arguments
Hi Guys I am Abhishek, primarily a bioinformatician. I have recently started using a lot of R thanks to some excellent packages available. Lately I have felt the need to batch process few of the R scripts I have been working with and strangely enough I am not able to find a good resource on how to best do this. I did find few old threads on the archives but none convinced me much. So here I am asking the same thing again hoping to get a good solution. 1. What's the best way to pass command line arguments to R scripts ? 2. How to execute R scripts from command line ? When I use R CMD BATCH I see no output on the screen 3. What does R --slave --vanilla do ? Thanks, -Abhi [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Writing R Scripts and passing command line arguments
See ?commandArgs, the getopt package and ?Rscript On Mon, Sep 7, 2009 at 1:47 PM, Abhishek Pratapabhishek@gmail.com wrote: Hi Guys I am Abhishek, primarily a bioinformatician. I have recently started using a lot of R thanks to some excellent packages available. Lately I have felt the need to batch process few of the R scripts I have been working with and strangely enough I am not able to find a good resource on how to best do this. I did find few old threads on the archives but none convinced me much. So here I am asking the same thing again hoping to get a good solution. 1. What's the best way to pass command line arguments to R scripts ? 2. How to execute R scripts from command line ? When I use R CMD BATCH I see no output on the screen 3. What does R --slave --vanilla do ? Thanks, -Abhi [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Adding to the middle of a vector
Colleagues, Assume that I have a vector with N elements Vector - 101:110). I would like to interpose an element into the vector at a position other than the start or end: c(101:104, 200, 105:110). I could do the following: NewVector - c(Vector[1:4], NewElement, Vector[5:10]) However, I suspect that there is a more elegant solution. Any proposals? Dennis Dennis Fisher MD P (The P Less Than Company) Phone: 1-866-PLessThan (1-866-753-7784) Fax: 1-415-564-2220 www.PLessThan.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] model simulation
Dear useRs, Is there a package or a function able to simulate models with sets of differential equations? Where we could input our model and give R some value to start with and it would generate the graphs? Regards, Rafael. [[elided Yahoo spam]] [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Adding to the middle of a vector
have a look at ?append(), e.g., x - 101:110 append(x, 200, after = 4) I hope it helps. Best, Dimitris Dennis Fisher wrote: Colleagues, Assume that I have a vector with N elements Vector- 101:110). I would like to interpose an element into the vector at a position other than the start or end: c(101:104, 200, 105:110). I could do the following: NewVector- c(Vector[1:4], NewElement, Vector[5:10]) However, I suspect that there is a more elegant solution. Any proposals? Dennis Dennis Fisher MD P (The P Less Than Company) Phone: 1-866-PLessThan (1-866-753-7784) Fax: 1-415-564-2220 www.PLessThan.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Dimitris Rizopoulos Assistant Professor Department of Biostatistics Erasmus University Medical Center Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands Tel: +31/(0)10/7043478 Fax: +31/(0)10/7043014 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] calling combinations of variable names
Hi, It's not enough to create a string with your instructions, it also needs to be evaluated as such. If you really wanted to evaluate your string, you'd need something like, a - b - cc - 1 # dummy example eval(parse(text = cbind(a, b, cc))) #library(fortunes) #fortune(parse) but fortune-ately, you probably oughtn't do that. Instead, I think you could use this example, d = data.frame(a1=1:10, b1=sin(1:10), c1 = cos(1:10)) # some data.frame variables = paste(letters[1:3],1,sep=) # its variable names (alternatively, names(d) would do) ICC = function(x) x # your function ICC(d[ variables[ c(1, 3) ] ]) # apply the function to a subset of the data, by selecting the required columns (If I understood correctly your question) HTH, baptiste 2009/9/7 Helter Two helter...@care2.com R-2.9.1, Windows7 Dear list, I have a question to you that seems very simple to me, but I just can't figure it out. I have a dataframe called ratings which contains the following variables: evalR1, evalR2, evalR3, evalR4, scoreR1, scoreR2, scoreR3, scoreR4, opinionR1, opinionR2, opinionR3, opinionR4. (there are more variables, but this gives an idea of the data structure). What I want is run several analyses on all 3 or 4-combinations of a given variable. So, for example, I want to compute the following ICC's (function from the psych package): ICC(cbind(evalR1,evalR2, evalR3)) ICC(cbind(evalR1,evalR2, evalR4)) ICC(cbind(evalR1, evalR3, evalR4)) ICC(cbind(evalR2, evalR3, evalR4)) ICC(cbind(evalR1, evalR2, evalR3, eval4)). I create a matrix containing the 3-combinations by combn(4,3). Now I need to call the variables into the function. First, I tried paste as follows: combis - combn(4,3) # this gives the 3-combinations attach(ratings) eval - paste(evalR,combis[1,1],,evalR,combis[2,1],,evalR,combis[3,1],sep =) (this is of course just for 1 combination, as an example) the output of this is evalR1,evalR2,evalR3, but when I run ICC(cbind(eval)), an error message is given which is not given when I enter ICC(cbind(evalR1,evalR2, evalR3)) manually. The function appears not to recognize the variable names. It also does not work to type ICC(cbind(unquote(eval))). Alternatively, I have tried the cat function, but also here ICC does not recognize the input as variable names. What am I doing wrong? How can I automatically construct the set of variable names such that a function recognizes them as variable names? ICC is one example, but there are also other computations to be run and the set of variables is pretty large, so typing the combinations of variable names manually is really unattractive. What am I missing? It seems to me that there probably is a very simple solution in R, but which? Thank you, Peter Verbeet [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- _ Baptiste Auguié School of Physics University of Exeter Stocker Road, Exeter, Devon, EX4 4QL, UK http://newton.ex.ac.uk/research/emag __ [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Writing R Scripts and passing command line arguments
Abhishek Pratap wrote: 1. What's the best way to pass command line arguments to R scripts ? As Gabor mentioned, the commandArgs function and the getopt package provide some excellent starting points for this. Abhishek Pratap wrote: 2. How to execute R scripts from command line ? When I use R CMD BATCH I see no output on the screen I believe R CMD BATCH dumps all of it's output to a file ending in .Rout. If you want more control over input and output to your script then Rscript is the utility to use. Abhishek Pratap wrote: 3. What does R --slave --vanilla do ? If I recall correctly, --vanilla makes it so that R does not waste time trying to load a previously saved session and also makes it so that R does not try to save history and environment variables when it exits. Vanilla also disables the loading of options from profile files such as ~/.Rprofile. I think --slave makes R shut up about it's self and run more quietly than it normally does. Hope that helps! -Charlie - Charlie Sharpsteen Undergraduate Environmental Resources Engineering Humboldt State University -- View this message in context: http://www.nabble.com/Writing-R-Scripts-and-passing-command-line-arguments-tp25334067p25334648.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] model simulation
Rafael Moral rafa_moral2004 at yahoo.com.br writes: Dear useRs, Is there a package or a function able to simulate models with sets of differential equations? Where we could input our model and give R some value to start with and it would generate the graphs? Regards, Rafael. install.packages(deSolve) it won't generate the graphs automatically, you have to use standard R tools (plot, matplot, etc.) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] model simulation
On Sep 7, 2009, at 1:55 PM, Rafael Moral wrote: Dear useRs, Is there a package or a function able to simulate models with sets of differential equations? Where we could input our model and give R some value to start with and it would generate the graphs? Your request seems a bit on the under-determined side, but perhaps this Petzold article will get you started: R as a Simulation Platform in Ecological Modelling http://www.cas.muohio.edu/%7Estevenmh/Rnews_2003-3.pdf ## pages 8-16 And then for more worked examples, the JSS volume that he edited: http://www.jstatsoft.org/v22 -- David Winsemius, MD Heritage Laboratories West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] linear mixed model question
The workers=as.factor(workers) codeline in my post dropped below my name. It should be in the code before the command line for the linear model. Daniel Malter wrote: Wen, to follow up on Thierry, your workers are nested in machines (since each worker only works one machine). Consider fitting a nested model. Though, with few observations, you might run into the same problem. Further, if you have observation triplets, and you expect systematic differences between each trial, then you would have to include a trial effect in some way. Have you plotted the data? This can be very informative. Finally, you may want to consider a fixed-effects approach. Throw in worker dummies only (without intercept) and test the linear hypothesis that the sum of the worker dummy coefficients is equal between two machines. The following example does that for two machines (if you have n machines you would have binomial coefficient (n,2) linear hypotheses): #simulate data machines=rep(c(0,1),each=9) workers=rep(c(1:6),each=3) workers.re=rep(rnorm(6),each=3) error=rnorm(18,0.3) output=2*machines+workers.re+error #code machines and workers as factors/dummies machines=as.factor(machines) #run OLS (fixed effects) of output on worker dummies without intercept fm4=lm(output~-1+workers) summary(fm4) #test the linear hypothesis that coefficients on the worker dummies are equal for both machines #the equivalent formulation used is: is the sum of the coefficients across machines equal to zero? library(car) linear.hypothesis(fm4,hypothesis.matrix=c(1,1,1,-1,-1,-1),rhs=c(0)) hope that helps, daniel workers=as.factor(workers) Wen Huang-3 wrote: Hello, I wanted to fit a linear mixed model to a data that is similar in terms of design to the 'Machines' data in 'nlme' package except that each worker (with triplicates) only operates one machine. I created a subset of observations from 'Machines' data such that it looks the same as the data I wanted to fit the model with (see code below). I fitted a model in which 'Machine' was a fixed effect and 'Worker' was random (intercept), which ran perfectly. Then I decided to complicate the model a little bit by fitting 'Worker' within 'Machine', which was saying variation among workers was nested within each machine. The model could be fitted by 'lme', but when I tried to get confidence intervals by 'intervals(fm2)' it gave me an error: Error in intervals.lme(fm2) : Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance I am wondering if this is because it is impossible to fit a model like 'fm2' or there is some other reasons? Thanks a lot! Wen # library(nlme) data(Machines) new.data = Machines[c(1:6, 25:30, 49:54), ] fm1 = lme(score ~ Machine, random = ~1|Worker, data = new.data) fm1 fm2 = lme(score ~ Machine, random = ~Machine-1|Worker, data = new.data) fm2 intervals(fm2) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- View this message in context: http://www.nabble.com/linear-mixed-model-question-tp25318961p25333981.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Confused - better empirical results with error in data
Hi, I have a strange one for the group. We have a system that predicts probabilities using a fairly standard svm (e1017). We are looking at probabilities of a binary outcome. The input data is generated by a perl script that calculates a bunch of things, fetches data from a database, etc. We train the system on 30,000 examples and then test the system on an unseen set of 5,000 records. The real world results on the test set looked VERY good. We were really happy with our model. The, we noticed that there was a big error in our data generation script and one of the values (an average of sorts.) was being calculated incorrectly. (The perl script failed to clear two iterators, so they both grew with every record.) As an quick experiment, we removed that item from our data set and re-ran the process. The results were not very good. Perhaps 75% as good as training with the wrong factor included. So, this is really a philosophical question. Do we: 1) Shrug and say, who cares, the SVM figured it out and likes that bad data item for some inexplicable reason 2) Tear into the math and try to figure out WHY the SVM is predicting more accurately Any opinions?? Thanks! __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Confused - better empirical results with error in data
Predicting whilst confused is unlikely to produce sound predictions...
my vote is for finding out why before believing anything.
Noah Silverman n...@smartmediacorp.com 09/07/09 8:33 PM
Hi,
I have a strange one for the group.
We have a system that predicts probabilities using a fairly standard svm
(e1017). We are looking at probabilities of a binary outcome.
The input data is generated by a perl script that calculates a bunch of
things, fetches data from a database, etc.
We train the system on 30,000 examples and then test the system on an
unseen set of 5,000 records.
The real world results on the test set looked VERY good. We were
really happy with our model.
The, we noticed that there was a big error in our data generation script
and one of the values (an average of sorts.) was being calculated
incorrectly. (The perl script failed to clear two iterators, so they
both grew with every record.)
As an quick experiment, we removed that item from our data set and
re-ran the process. The results were not very good. Perhaps 75% as
good as training with the wrong factor included.
So, this is really a philosophical question. Do we:
1) Shrug and say, who cares, the SVM figured it out and likes
that bad data item for some inexplicable reason
2) Tear into the math and try to figure out WHY the SVM is
predicting more accurately
Any opinions??
Thanks!
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
***
This email and any attachments are confidential. Any use...{{dropped:8}}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Omnibus test for main effects in the face of an interaction containing the main effects.
R 2.9.1
Windows XP
I am fitting a random effects ANOVA with two factors Group which has two levels
and Time which has three levels:
fita-lme(Post~Time+factor(Group)+factor(Group)*Time,
random=~1|SS,data=blah$alldata)
I want to get the omnibus significance tests for each factor and the
interaction. I believe I can get the omnibus test for the interaction by
running the model:
fitb-lme(Post~Time+factor(Group), random=~1|SS,data=blah$alldata)
followed by
anova(fita,fitb).
How do I get the omnibus test for the main effects i.e. for Time and
factor(Group)? I could drop each from the model, i.e.
fitc-lme(Post~ factor(Group)+factor(Group)*Time,
random=~1|SS,data=blah$alldata)
fitd-lme(Post~Time+factor(Group)*Time,
random=~1|SS,data=blah$alldata)
and then run
anova(fita,fitc)
anova(fita,fitd)
but I don't like this option as it will have in interaction that contains a
factor that is not included in the model as a main effect. How then do I get
the omnibus test for Time and factor(Group)?
Thanks
John
John David Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
Confidentiality Statement:
This email message, including any attachments, is for th...{{dropped:6}}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] xyplot {lattice} are different types possible for each panel?
Hello R Folks...
Using the example below, I¹d like two of the panels to be plotted with type
= ³p² but the third to be done with type = ³h². I can¹t use type = c(³p²,
³p², ³h²) because this syntax applies all given types to every panel. I
don¹t think I can use groups and distribute.type because these are intended
for different styles of plotting within a single panel. As you can see, I
tried to do a panel function following something I saw in the Lattice book,
but this has no effect at all. Looks like it may have to be more elaborate,
but I¹m stuck. Any suggestions appreciated!
Thanks, Bryan
*
Bryan Hanson
Professor of Chemistry Biochemistry
DePauw University, Greencastle IN USA
y - rnorm(100)
x - rnorm(100)
names - rep(c(Set 1, Set 2, Set 3), 4)
df - data.frame(y = y, x = y, names = as.factor(names))
p - xyplot(y ~ x | names,
layout = c(1, 3),
panel = function(...) {
panel.xyplot(...)
if (panel.number() == 1) type = h
})
plot(p)
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Omnibus test for main effects in the face of an interaction containing the main effects.
R 2.9.1
Windows XP
UPDATE,
Even my first suggestion
anova(fita,fitb) is probably not appropriate as the fixed effects are different
in the two model, so I don't even know how to perform the ombnibus test for the
interaction!
I am fitting a random effects ANOVA with two factors Group which has two levels
and Time which has three levels:
fita-lme(Post~Time+factor(Group)+factor(Group)*Time,
random=~1|SS,data=blah$alldata)
I want to get the omnibus significance tests for each factor and the
interaction. I believe I can get the omnibus test for the interaction by
running the model:
fitb-lme(Post~Time+factor(Group), random=~1|SS,data=blah$alldata)
followed by
anova(fita,fitb).
How do I get the omnibus test for the main effects i.e. for Time and
factor(Group)? I could drop each from the model, i.e.
fitc-lme(Post~ factor(Group)+factor(Group)*Time,
random=~1|SS,data=blah$alldata)
fitd-lme(Post~Time+factor(Group)*Time,
random=~1|SS,data=blah$alldata)
and then run
anova(fita,fitc)
anova(fita,fitd)
but I don't like this option as it will have in interaction that contains a
factor that is not included in the model as a main effect. How then do I get
the omnibus test for Time and factor(Group)?
Thanks
John
John David Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
Confidentiality Statement:
This email message, including any attachments, is for th...{{dropped:6}}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Confused - better empirical results with error in data
On Mon, Sep 7, 2009 at 12:33 PM, Noah Silvermann...@smartmediacorp.com wrote: SNIP So, this is really a philosophical question. Do we: 1) Shrug and say, who cares, the SVM figured it out and likes that bad data item for some inexplicable reason 2) Tear into the math and try to figure out WHY the SVM is predicting more accurately Any opinions?? Thanks! Boy, I'd sure think you'd want to know why it worked with the 'wrong' calculations. It's not that the math is wrong, really, but rather that it wasn't what you thought it was. I cannot see why you wouldn't want to know why this mistake helped. Won't future project benefit? Just my 2 cents, Mark __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R-crash when loading workspace - Windows
On 07/09/2009 10:34 AM, sebed1110-div...@yahoo.fr wrote: Dear all, One day when I tried to load an existing workspace (when opening R or by load()), R crashed without any error notification. The day before I had worked and saved my workspace without any trouble. At first I though it was a memory problem (workspace reaching 180Mo) or related to a particular script or command, so I start a new workspace. Everything was ok, that script and others working. Then I saved the workspace (55Mo) and tried to open it, without any result : R crashes without any notification again. This occurs only with Windows. Does someone know how to solve that problem? You need to come up with a simple scheme to reproduce it, and then someone else will debug it, or you need to debug it yourself. If you're willing to let me download the saved workspace, I could give it a try. (I can't accept an email that big, you'll need to put it up for http or ftp download.) Duncan Murdoch __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R-crash when loading workspace - Windows
On Mon, Sep 7, 2009 at 1:08 PM, Duncan Murdochmurd...@stats.uwo.ca wrote: On 07/09/2009 10:34 AM, sebed1110-div...@yahoo.fr wrote: Dear all, One day when I tried to load an existing workspace (when opening R or by load()), R crashed without any error notification. The day before I had worked and saved my workspace without any trouble. At first I though it was a memory problem (workspace reaching 180Mo) or related to a particular script or command, so I start a new workspace. Everything was ok, that script and others working. Then I saved the workspace (55Mo) and tried to open it, without any result : R crashes without any notification again. This occurs only with Windows. Does someone know how to solve that problem? You need to come up with a simple scheme to reproduce it, and then someone else will debug it, or you need to debug it yourself. If you're willing to let me download the saved workspace, I could give it a try. (I can't accept an email that big, you'll need to put it up for http or ftp download.) Duncan Murdoch Wouldn't urt to uninstall and then reinstall Rgui. (I'm assuming this is Rgui?) Maybe a file got corrupted and a reinstall will clear things up? Cheers, Mark __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] xyplot {lattice} are different types possible for each panel?
Hi,
Something like this perhaps,
p - xyplot(y ~ x | names,
layout = c(1, 3),
panel = function(...,type=p) {
if (panel.number() == 1) {
panel.xyplot(...,type = h)
} else {
panel.xyplot(...,type = type)
}
})
plot(p)
HTH,
baptiste
2009/9/7 Bryan Hanson han...@depauw.edu
Hello R Folks...
Using the example below, IÄ
d like two of the panels to be plotted with type
= ÅpË but the third to be done with type = ÅhË. I canÄ
t use type =
c(ÅpË,
ÅpË, ÅhË) because this syntax applies all given types to every panel. I
donÄ
t think I can use groups and distribute.type because these are intended
for different styles of plotting within a single panel. As you can see, I
tried to do a panel function following something I saw in the Lattice book,
but this has no effect at all. Looks like it may have to be more
elaborate,
but IÄ
m stuck. Any suggestions appreciated!
Thanks, Bryan
*
Bryan Hanson
Professor of Chemistry Biochemistry
DePauw University, Greencastle IN USA
y - rnorm(100)
x - rnorm(100)
names - rep(c(Set 1, Set 2, Set 3), 4)
df - data.frame(y = y, x = y, names = as.factor(names))
p - xyplot(y ~ x | names,
layout = c(1, 3),
panel = function(...) {
panel.xyplot(...)
if (panel.number() == 1) type = h
})
plot(p)
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
_
Baptiste Auguié
School of Physics
University of Exeter
Stocker Road,
Exeter, Devon,
EX4 4QL, UK
http://newton.ex.ac.uk/research/emag
__
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Confused - better empirical results with error in data
You both make good points. Ideally, it would be nice to know WHY it works. Without digging into too much verbiage, the system is designed to predict the outcome of certain events. The broken model predicts outcomes correctly much more frequently than one with the broken data withheld. So, to answer Mark's question, we say it's better because we see much better results with our broken model when applied to real-world data used for testing. I have one theory. The data is listed in our CSV file from newest to oldest. We are supposed to calculated a valued that is an average of some items. We loop through some queries to our database and increment two variables - $total_found and $total_score. The final value is simply $total_score / $total_found. Our programmer forgot to reset both $total_score and $total_found back to zero for each record we process. So both grow. I think that this may, in a way, be some warped form of a recency weighted score. The newer records will have a score more affected by their contribution to the wrongly growing totals. A record that is closer to the end of the data set will be starting with HUGE values for $total_score and $total_found, so addition of its values will have very little effect. We've done the following so far today (Note, scores are just relative to indicate performance. Higher is better) 1) Run with bad data = 6.9 2) Run with bad data missing = 5.5 3) Run with correct data = ?? (We're running now, will take a few hours to compute.) I might also try to plot the bad data. It would be interesting to see what shape it has... On 9/7/09 1:05 PM, Mark Knecht wrote: On Mon, Sep 7, 2009 at 12:33 PM, Noah Silvermann...@smartmediacorp.com wrote: SNIP So, this is really a philosophical question. Do we: 1) Shrug and say, who cares, the SVM figured it out and likes that bad data item for some inexplicable reason 2) Tear into the math and try to figure out WHY the SVM is predicting more accurately Any opinions?? Thanks! Boy, I'd sure think you'd want to know why it worked with the 'wrong' calculations. It's not that the math is wrong, really, but rather that it wasn't what you thought it was. I cannot see why you wouldn't want to know why this mistake helped. Won't future project benefit? Just my 2 cents, Mark __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Using anova(f1, f2) to compare lmer models yields seemingly erroneous Chisq = 0, p = 1
On Mon, Sep 7, 2009 at 10:34 AM, Alain Zuurhighs...@highstat.com wrote: rapton wrote: Hello, I am using R to analyze a large multilevel data set, using lmer() to model my data, and using anova() to compare the fit of various models. When I run two models, the output of each model is generated correctly as far as I can tell (e.g. summary(f1) and summary(f2) for the multilevel model output look perfectly reasonable), and in this case (see below) predictor.1 explains vastly more variance in outcome than predictor.2 (R2 = 15% vs. 5% in OLS regression, with very large N). What I am utterly puzzled by is that when I run an anova comparing the two multilevel model fits, the Chisq comes back as 0, with p = 1. I am pretty sure that fit #1 (f1) is a much better predictor of the outcome than f2, which is reflected in the AIC, BIC , and logLik values. And, unless I'm missing something, also by the (misspecified) test. A large p-value indicates you have no evidence that the additional 19 parameters in f2 improve fit, which matches what the other methods suggested. However, as has been pointed out, the lack of nesting makes this a faulty LRT. This is made apparent by the fact that you get a test statistic outside the support of the chi-squared distribution (positive reals) (lambda - (-2)*(-22715 - (-23633))) [1] -1836 and since the test is uses right-tail probability, anova is not changing anything by moving the statistic to 0. pchisq(lambda, 19, lower = FALSE) [1] 1 pchisq(0, 19, lower = FALSE) [1] 1 To do the test properly the restricted (null) model must be a special case of the general (alternative) model (e.g., with the additional 19 parameters set to zero) which will result in the null model having a smaller likelihood, leading to a positive tests statistic. When that statistic is small you get a large p-value indicating a lack of evidence that the additional parameters improve fit... hth, Kingsford Why might anova be giving me this curious output? How can I fix it? I am sure I am making a dumb error somewhere, but I cannot figure out what it is. Any help or suggestions would be greatly appreciated! -Matt f1 - (lmer(outcome ~ predictor.1 + (1 | person), data=i)) f2 - (lmer(outcome ~ predictor.2 + (1 | person), data=i)) anova(f1, f2) Data: i Models: f1: outcome ~ predictor.1 + (1 | person) f2: outcome ~ predictor.2 + (1 | person) Df AIC BIC logLik Chisq Chi Df Pr(Chisq) f1 6 45443 45489 -22715 f2 25 47317 47511 -23633 0 19 1 ** NOT ** nested sorrythe brain is going faster than the fingers. - Dr. Alain F. Zuur First author of: 1. Analysing Ecological Data (2007). Zuur, AF, Ieno, EN and Smith, GM. Springer. 680 p. 2. Mixed effects models and extensions in ecology with R. (2009). Zuur, AF, Ieno, EN, Walker, N, Saveliev, AA, and Smith, GM. Springer. 3. A Beginner's Guide to R (2009). Zuur, AF, Ieno, EN, Meesters, EHWG. Springer Statistical consultancy, courses, data analysis and software Highland Statistics Ltd. 6 Laverock road UK - AB41 6FN Newburgh Email: highs...@highstat.com URL: www.highstat.com -- View this message in context: http://www.nabble.com/Using-anova%28f1%2C-f2%29-to-compare-lmer-models-yields-seemingly-erroneous-Chisq-%3D-0%2C-p-%3D-1-tp25297254p25333148.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Size of plots in pdf files#can it be smaller?
Hi all, I have to produce arrangements of 25 simple plots of the type plot(x,y,pch=.) where there are typically on the order of 2 points. So, overall, I have about 50 points. When I use the pdf device, I get file sizes (on a Windows machine) of about 10 MB. When I then zip the files, I'm down to about 0.5MB, so the original pdf files were created with a lot of 'air'. I'm wondering why the pdf files are so large and whether there is an alternative to produce them smaller. Any ideas? Have a nice afternoon, Chris. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] avoid NA in list
Henrique and Peter, thank you very much for your advices :) -- View this message in context: http://www.nabble.com/avoid-NA-in-list-tp25321020p25330422.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Usage of OCaml/R binding.
Hello. I've been pulling together a Debian package out of Maxence Guesdon's OCaml bindings for R. Will be available from my website as soon as I get my router to obey me. Here's Maxence's bindings: http://pauillac.inria.fr/~guesdon/ocaml-r.en.html The purpose of this software is to access R from OCaml programs. However, my issue is that after having pulled things to a Debian package, I am completely unfamiliar with the API, and I'm having trouble initialising this module. You'll find below a list of the functions accessible from Maxence Guesdon's R module (there's also a Rmath module for the math library which I'm not yet concerned about). I'd specifically appreciate information on how to get started using these R API functions. If you're scared about OCaml, simply look at the two declarations below in the following way: external init_r : string array - int = init_r 'external' means that this value is to be fetched by linking C code identified by the '= init_r' ending statement. Typically, 'external' values wrap the API functions. So this value is called init_r, and takes an array of strings as first argument, and yields back an integer. val init : ?argv:string array - unit - int 'val' means that it is a standard OCaml value. It's called init. And its type is 'I take an array of strings as first argument, then a unit ( think of it a C void value containing nothing), and I get back an integer. The '?argv:' stuff is concerned with optional arguments, so no real big deal. Documentation and examples of how the API works to embed R in an application would be highly appreciated. Please forward this email to R-dev if you feel it belongs there. All the best, Guillaume Yziquel. yziq...@seldon:~/sandbox/repo/debian/debian-ocaml/ocaml-r$ ocaml-batteries Objective Caml version 3.11.1 _ | | | | [| + | | Batteries Included - | |___|_|___| _ | | | | | -Type '#help;;' | | + |] |___|_|___| # #rectypes;; # #require R;; # module X = R;; module X : sig external init_r : string array - int = init_r val init : ?argv:string array - unit - int external terminate : unit - unit = end_r type sexp = R.sexp type symbol = string type arg = [ `Anon of sexp | `Named of sexp * symbol ] external sexp : string - sexp = r_sexp_of_string external sexp_of_symbol : symbol - sexp = r_sexp_of_symbol external set_var : symbol - sexp - unit = r_set_var external r_print_value : sexp - unit = r_print_value external exec : string - arg array - unit = r_exec external current_test : unit - unit = r_current_test external to_bool : sexp - bool = bool_of_sexp external to_int : sexp - int = int_of_sexp external to_float : sexp - float = float_of_sexp external to_string : sexp - string = string_of_sexp external of_bool : bool - sexp = sexp_of_bool external of_int : int - sexp = sexp_of_int external of_float : float - sexp = sexp_of_float external of_string : string - sexp = sexp_of_string external to_bool_array : sexp - bool array = bool_array_of_sexp external to_int_array : sexp - int array = int_array_of_sexp external to_float_array : sexp - float array = float_array_of_sexp external to_string_array : sexp - string array = string_array_of_sexp external of_bool_array : bool array - sexp = sexp_of_bool_array external of_int_array : int array - sexp = sexp_of_int_array external of_float_array : float array - sexp = sexp_of_float_array external of_string_array : string array - sexp = sexp_of_string_array external get_attrib : sexp - string - sexp = r_get_attrib val dim : sexp - int array val dimnames : sexp - string array end # -- Guillaume Yziquel http://yziquel.homelinux.org/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [R-pkgs] amer: generalized additive mixed models with lme4
Dear R-users, I'd like to announce the release of the amer-package that adds the capability to fit generalized additive mixed models to lme4. It includes a vignette with real data examples and a brief summary of the theory behind the implementation. Best, Fabian [[alternative HTML version deleted]] ___ R-packages mailing list r-packa...@r-project.org https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] How can I change characteristics of a cca biplot in R
Hi, I´m doing cca for a community data set in R and I have made a biplot for my data. Otherwise everything seems to be allright but the biplot is so messy I can´t read it well enough or publish it. I would like to get the row numbers out of the plot: I want the species position and the environmental variables in it as vectors but not the station numbers. How do I get them out? I have a raw species data set and a raw environmental data set and I don´t have station names or numbers in the data set but R puts the row names in anyway. I understand they are necessary but I just don´t want to show them in the biplot. Thank you! Yours sincerely Heta Rousi ( Finnish environment institute, Marine center) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Sort an array
Hi, Would anybody know how to sort an array in order? I basically store the results from an analysis in an array and would like to organise it at the end of my loop with the lowest result at the index 1 and the highest result at the last index. Thanks. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] List of tags in roxygen and use for S4 classes?
I would also wish for a better (online) documentation, as I think the general idea of roxygen is great. I agree completely. Good call; the vignette is terse and outdated. Manuel and I are in the process of preparing a paper based on our DSC talks; that should fill in some of the details w.r.t. S4. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] How to reduce memory demands in a function?
I've written a function that regularly throws the cannot allocate vector of size X Kb error, since it contains a loop that creates large numbers of big distance matrices. I'd be very grateful for any simple advice on how to reduce the memory demands of my function. Besides increasing memory.size to the maximum available, I've tried reducing my dist objects to 3 sig.. fig.s (not sure if that made any difference), I've tried the distance function daisy() from package cluster instead of dist(), and I've avoided storing unnecessary intermediary objects as far as possible by nesting functions in the same command. I've even tried writing each of my dist() objects to a text file, one line for each, and reading them in again one at a time as and when required, using scan() - and although this seemed to avoid the memory problem, it ran so slowly that it wasn't much use for someone with deadlines to meet... I don't have formal training in programming, so if there's something handy I should read, do let me know. Thanks, Richard Gunton. Postdoctoral researcher in arable weed ecology, INRA Dijon. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] using an array of strings with strsplit, issue when including a space in split criteria
Dear all, I'm having a problem understanding why a split does not occur with in the 2nd use of the function strsplit below: # text strings txt - c(sales to 23 August 2008 published 29 August, + sales to 6 September 2008 published 11 September) # first use strsplit(txt, 'published', fixed=TRUE) [[1]] [1] sales to 23 August 2008 29 August [[2]] [1] sales to 6 September 2008 11 September # second use, but with a space ' ' in the split strsplit(txt, 'published ', fixed=TRUE) [[1]] [1] sales to 23 August 2008 29 August [[2]] [1] sales to 6 September 2008 published 11 September Thank you kindly for any help in advance. Tony O/S: Win Vista Ultimate sessionInfo() R version 2.9.2 (2009-08-24) i386-pc-mingw32 locale: LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United Kingdom. 1252;LC_MONETARY=English_United Kingdom. 1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] RODBC_1.3-0 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Confused - better empirical results with error in data
You both make good points. Ideally, it would be nice to know WHY it works. Without digging into too much verbiage, the system is designed to predict the outcome of certain events. The broken model predicts outcomes correctly much more frequently than one with the broken data withheld. So, to answer Mark's question, we say it's better because we see much better results with our broken model when applied to real-world data used for testing. I have one theory. The data is listed in our CSV file from newest to oldest. We are supposed to calculated a valued that is an average of some items. We loop through some queries to our database and increment two variables - $total_found and $total_score. The final value is simply $total_score / $total_found. Our programmer forgot to reset both $total_score and $total_found back to zero for each record we process. So both grow. I think that this may, in a way, be some warped form of a recency weighted score. The newer records will have a score more affected by their contribution to the wrongly growing totals. A record that is closer to the end of the data set will be starting with HUGE values for $total_score and $total_found, so addition of its values will have very little effect. We've done the following so far today (Note, scores are just relative to indicate performance. Higher is better) 1) Run with bad data = 6.9 2) Run with bad data missing = 5.5 3) Run with correct data = ?? (We're running now, will take a few hours to compute.) I might also try to plot the bad data. It would be interesting to see what shape it has... On 9/7/09 1:05 PM, Mark Knecht wrote: On Mon, Sep 7, 2009 at 12:33 PM, Noah Silvermann...@smartmediacorp.com wrote: SNIP So, this is really a philosophical question. Do we: 1) Shrug and say, who cares, the SVM figured it out and likes that bad data item for some inexplicable reason 2) Tear into the math and try to figure out WHY the SVM is predicting more accurately Any opinions?? Thanks! Boy, I'd sure think you'd want to know why it worked with the 'wrong' calculations. It's not that the math is wrong, really, but rather that it wasn't what you thought it was. I cannot see why you wouldn't want to know why this mistake helped. Won't future project benefit? Just my 2 cents, Mark __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how do I draw this surface -- hand drawn in the attachemtn
for example... x - y - seq(.1, 2, .1) ftn - function(x, y) x - .25*x^2 z - outer(x, y, ftn) persp(x, y, z, theta = 330, phi = 30) hth, Kingsford On Sun, Sep 6, 2009 at 8:28 AM, Steve Lianogloumailinglist.honey...@gmail.com wrote: Hi, On Sat, Sep 5, 2009 at 5:16 AM, gallon ligallon...@gmail.com wrote: Basially I have the observation for x, y, and z for each fixed z, the value of y moves from 0 to y0 as x moves from 0 to x0. I then move the curve for z from 0 to z0 and form a continuous surface. want to know if there is a way that R can replace my handdrawing. Do the examples in ?persp help? How about this example from the r graph gallery: http://addictedtor.free.fr/graphiques/graphcode.php?graph=42 I think lattice graphics also has functions to help plotting in 3d, but I've never used it. I've only seen it in actio in the lattics vs. ggplot series of posts on the learnr.wordpress.com blog (very helpful). -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Size of plots in pdf files#can it be smaller?
I have to produce arrangements of 25 simple plots of the type plot(x,y,pch=.) where there are typically on the order of 2 points. So, overall, I have about 50 points. When I use the pdf device, I get file sizes (on a Windows machine) of about 10 MB. When I then zip the files, I'm down to about 0.5MB, so the original pdf files were created with a lot of 'air'. I'm wondering why the pdf files are so large and whether there is an alternative to produce them smaller. Any ideas? As far as I know, R cannot directly produce compressed pdf files. You could postprocess your plots. E.g. pdftk is quite conveniant: pdftk foo.pdf cat output foo2.pdf compress Maybe in a loop if many files are involved. cu Philipp -- Dr. Philipp Pagel Lehrstuhl für Genomorientierte Bioinformatik Technische Universität München Wissenschaftszentrum Weihenstephan 85350 Freising, Germany http://webclu.bio.wzw.tum.de/~pagel/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] lean text label below barplot table
Hi, everyone: I am plotting an graph with bar plot, but the label after every bar is too long, I wanna if I can draw the label lean to an angle thanks -- Xiaogang Yang Sensorweb Research Laboratory http://sensorweb.vancouver.wsu.edu/ Washington State University Vancouver [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How Does One Use the Value of Density Function?
On Mon, Sep 7, 2009 at 12:42 AM, Gundala Viswanathgunda...@gmail.com wrote: How do people usually use the result of density function (e.g. dnorm)? Especially when its value can be greater than 1. What do they do with such density 1? dnorm(2.02,2,.24) [1] 1.656498 There are countless uses. E.g., when viewed as a function of the parameters it is a likelihood, allowing one to estimate parameters given data... hth, Kingsford - G.V. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] finding the minimum value
Hi,
If I correctly understand what you want to do, then this might work:
tmp - rep(0, times=1000)
for (i in 1:1000) tmp[i] - get(paste('B', i, sep=''))
mean(tmp)
HTH,
Vik
maram salem wrote:
Hi all,
I'm using a certain procedure to calculate the value of some
variable(Bayes risk),B.
So I got the values B1, B2, , B1000, each under certain input
values and using a long procedure.
Now, I want to put the values I got in a nummerical vector and find their
minimum value. I think c( ) should work.For example if I have only 10
values I could have used
c(B1,B2,B3,B4,B5,B6,B7,B8,B9,B10)
But how can I do this for 1000 values?
I think the question is really trivial, but I tried to work on it and
couldn't reach anythg.
Thanks.
Maram
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
View this message in context:
http://www.nabble.com/finding-the-minimum-value-tp25332756p25336385.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] finding the minimum value
Sorry there was a typo in the code. This should work:
tmp - rep(0, times=1000)
for (i in 1:1000) tmp[i] - get(paste('B', i, sep=''))
mean(tmp)
Viknesh wrote:
Hi,
If I correctly understand what you want to do, then this might work:
tmp - rep(0, times=1000)
for (i in 1:1000) tmp[i] - get(paste('B', i, sep=''))
mean(tmp)
HTH,
Vik
maram salem wrote:
Hi all,
I'm using a certain procedure to calculate the value of some
variable(Bayes risk),B.
So I got the values B1, B2, , B1000, each under certain input
values and using a long procedure.
Now, I want to put the values I got in a nummerical vector and find their
minimum value. I think c( ) should work.For example if I have only 10
values I could have used
c(B1,B2,B3,B4,B5,B6,B7,B8,B9,B10)
But how can I do this for 1000 values?
I think the question is really trivial, but I tried to work on it and
couldn't reach anythg.
Thanks.
Maram
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
View this message in context:
http://www.nabble.com/finding-the-minimum-value-tp25332756p25336425.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] lean text label below barplot table
On Sep 7, 2009, at 5:03 PM, Xiaogang Yang wrote: Hi, everyone: I am plotting an graph with bar plot, but the label after every bar is too long, I wanna if I can draw the label lean to an angle thanks It depends on the particular function and bar plot is insufficiently specified to be specific. If it's a Lattice function, you should be looking at the rot (sub) parameter within the scales parameters. If it's a graphics function, you could try looking at the las parameter to par(). -- David Winsemius, MD Heritage Laboratories West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] spplot modifications
http://www.nabble.com/file/p25336596/Conductivity1.jpeg I need a little help making modifications to the image included with this post. First, rather than using a linear color legend to display the output I would like to use a log-scale legend. Thus, the legend on the right would go from 1 to 1000, for example, with a classic log-scale gradation. What I hope to avoid is taking the log of the values and displaying the logged values using a linear scale. Second, at the top of the legend there is a random green bar that shouldn't be there. If you know what causes this or more importantly, how to get rid of it, please help. I've included the code and text files necessary to create this figure in the zip file. Thank you, Eric Morway http://www.nabble.com/file/p25336596/Morway_R_Qs.zip Morway_R_Qs.zip -- View this message in context: http://www.nabble.com/spplot-modifications-tp25336596p25336596.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Confused - better empirical results with error in data
On Mon, Sep 7, 2009 at 1:22 PM, Noah Silvermann...@smartmediacorp.com wrote: SNIP The data is listed in our CSV file from newest to oldest. We are supposed to calculated a valued that is an average of some items. We loop through some queries to our database and increment two variables - $total_found and $total_score. The final value is simply $total_score / $total_found. SNIP This does seem like it's rife with possibilities for non-causal action. (Assuming you process from newest toward oldest which is what I think you say you are doing...) I'm pretty sure that if I knew that the Dow was going to be higher 3 months from now then my day trading results would tend toward long vs short and I'd do better. Unfortunately I don't know where it will be and cannot really do that. Have you considered processing the data in the other direction. Not in R, but rather reversing the data frame or better yet writing the csv file in date order? Cheers, Mark __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Sort an array
?order Possibly something like A = A[order(A$Field1, A$Field2),] On Mon, Sep 7, 2009 at 3:22 AM, OLIVIER REGNIER-COUDERT (0509785)o.regnier-coud...@rgu.ac.uk wrote: Hi, Would anybody know how to sort an array in order? I basically store the results from an analysis in an array and would like to organise it at the end of my loop with the lowest result at the index 1 and the highest result at the last index. Thanks. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] spplot modifications
On Sep 7, 2009, at 5:29 PM, emorway wrote: http://www.nabble.com/file/p25336596/Conductivity1.jpeg I need a little help making modifications to the image included with this post. First, rather than using a linear color legend to display the output I would like to use a log-scale legend. Thus, the legend on the right would go from 1 to 1000, for example, with a classic log-scale gradation. What I hope to avoid is taking the log of the values and displaying the logged values using a linear scale. I don't know what you mean here, (classic not being a particularly precise term) but suggest you look at the chapter 8 examples Figs 8.3 - 8.5 in Sarkar's website: http://lmdvr.r-forge.r-project.org/figures/figures.html Second, at the top of the legend there is a random green bar that shouldn't be there. The random green bar is there because colors are recycled and you have more levels than colors. If you know what causes this or more importantly, how to get rid of it, please help. I've included the code and text files necessary to create this figure in the zip file. Thank you, Eric Morway http://www.nabble.com/file/p25336596/Morway_R_Qs.zip Morway_R_Qs.zip -- View this message in context: http://www.nabble.com/spplot-modifications-tp25336596p25336596.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius, MD Heritage Laboratories West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Confused - better empirical results with error in data
Interesting point. Our data is NOT continuous. Sure, some of the test examples are older than others, but there is no relationship between them. (More Markov like in behavior.) When creating a specific record, we actually account for this in our SQL queries which tend to be along the lines of: select x from table where id=1234 and date '2008-05-01' This way, whatever data we're looking at, we set things so the current and future data doesn't exist yet. My understanding was that an SVM wouldn't care about the order of the data input as long as the examples are independent. Regardless of all this, we look at real-world test for our evaluation. 1) We trained the system on examples prior to a certain date. 2) We test the system with unseen examples after that date. We take the approach of: If we had used this model, what would our portfolio be at the end of the test period. Sure, we also look at things like AUC and R2 (from applying the model to the TEST data.) Generally, we see a correlation between AUC, R2, and our final result, but not a perfect one. A model with a SLIGHTLY lower R2 actually produced better results in a few cases. This process should produce solid results as we are eliminating any chance of over-fitting when measuring performance. So, one could argue, that whatever gives the best results on the test data is the best model, regardless of the correctness of the theory. Just for fun, I'll see if I can schedule a few hours to run the same experiment with the training data order reversed. If I'm correct, the results should be the same. Thanks! -- N On 9/7/09 2:34 PM, Mark Knecht wrote: On Mon, Sep 7, 2009 at 1:22 PM, Noah Silvermann...@smartmediacorp.com wrote: SNIP The data is listed in our CSV file from newest to oldest. We are supposed to calculated a valued that is an average of some items. We loop through some queries to our database and increment two variables - $total_found and $total_score. The final value is simply $total_score / $total_found. SNIP This does seem like it's rife with possibilities for non-causal action. (Assuming you process from newest toward oldest which is what I think you say you are doing...) I'm pretty sure that if I knew that the Dow was going to be higher 3 months from now then my day trading results would tend toward long vs short and I'd do better. Unfortunately I don't know where it will be and cannot really do that. Have you considered processing the data in the other direction. Not in R, but rather reversing the data frame or better yet writing the csv file in date order? Cheers, Mark [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] xyplot {lattice} are different types possible for each panel?
Thanks Baptiste, your suggestion works wonderfully. Bryan
For anyone following along, the following line needs to replace the similar
one in my original example:
names - rep(c(Set 1, Set 2, Set 3, Set 4), 25)
Or the data lengths will be wrong.
On 9/7/09 4:19 PM, baptiste auguie baptiste.aug...@googlemail.com wrote:
Hi,
Something like this perhaps,
p - xyplot(y ~ x | names,
  layout = c(1, 3),
  panel = function(...,type=p) {
      if (panel.number() == 1) {
        panel.xyplot(...,type = h)
         } else {
        panel.xyplot(...,type = type)
        }
      })
plot(p)
HTH,
baptiste
2009/9/7 Bryan Hanson han...@depauw.edu
Hello R Folks...
Using the example below, IÄ
d like two of the panels to be plotted with type
= ÅpË but the third to be done with type = ÅhË. Â I canÄ
t use type =
c(ÅpË,
ÅpË, ÅhË) because this syntax applies all given types to every panel  I
donÄ
t think I can use groups and distribute.type because these are intended
for different styles of plotting within a single panel. Â As you can see, I
tried to do a panel function following something I saw in the Lattice book,
but this has no effect at all. Â Looks like it may have to be more elaborate,
but IÄ
m stuck. Â Any suggestions appreciated!
Thanks, Bryan
*
Bryan Hanson
Professor of Chemistry Biochemistry
DePauw University, Greencastle IN USA
y - rnorm(100)
x - rnorm(100)
names - rep(c(Set 1, Set 2, Set 3), 4)
df - data.frame(y = y, x = y, names = as.factor(names))
p - xyplot(y ~ x | names,
  layout = c(1, 3),
  panel = function(...) {
    panel.xyplot(...)
    if (panel.number() == 1) type = h
    })
plot(p)
    [[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Size of plots in pdf files#can it be smaller?
The reason for the size of the file is that when generating a PDF, commands are generated to plot each point. I generated a PDF file with 10,000 and there were 10,000 of lines similar to the following: 159.93 349.01 1.00 1.00 re f 343.19 283.07 1.00 1.00 re f 427.86 323.58 1.00 1.00 re f 431.68 230.08 1.00 1.00 re f 93.79 278.89 1.00 1.00 re f 425.78 332.10 1.00 1.00 re f 332.04 366.16 1.00 1.00 re f 78.55 305.61 1.00 1.00 re f 277.22 135.42 1.00 1.00 re f 423.47 101.07 1.00 1.00 re f 435.86 289.67 1.00 1.00 re f My question is why do you need so many points? Have you looked at 'hexbin' as way of plotting the data? Have you considered generating 'jpg' output which will be much smaller; e.g. a jpg file with 10,000 points was 70K and one with 100,000 was 90K -- it only created the pixels to be plotted. On Mon, Sep 7, 2009 at 7:43 AM, Christian Ritterchristian.rit...@uclouvain.be wrote: Hi all, I have to produce arrangements of 25 simple plots of the type plot(x,y,pch=.) where there are typically on the order of 2 points. So, overall, I have about 50 points. When I use the pdf device, I get file sizes (on a Windows machine) of about 10 MB. When I then zip the files, I'm down to about 0.5MB, so the original pdf files were created with a lot of 'air'. I'm wondering why the pdf files are so large and whether there is an alternative to produce them smaller. Any ideas? Have a nice afternoon, Chris. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Jim Holtman Cincinnati, OH +1 513 646 9390 What is the problem that you are trying to solve? __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Size of plots in pdf files#can it be smaller?
No (destructive) JPGs - they are evil (draw a line as see for yourself) and should be banned from publications (only useful for pictures/photos). Use PNGs for you plots if you don't like vector graphics. /H On Mon, Sep 7, 2009 at 5:00 PM, jim holtmanjholt...@gmail.com wrote: The reason for the size of the file is that when generating a PDF, commands are generated to plot each point. I generated a PDF file with 10,000 and there were 10,000 of lines similar to the following: 159.93 349.01 1.00 1.00 re f 343.19 283.07 1.00 1.00 re f 427.86 323.58 1.00 1.00 re f 431.68 230.08 1.00 1.00 re f 93.79 278.89 1.00 1.00 re f 425.78 332.10 1.00 1.00 re f 332.04 366.16 1.00 1.00 re f 78.55 305.61 1.00 1.00 re f 277.22 135.42 1.00 1.00 re f 423.47 101.07 1.00 1.00 re f 435.86 289.67 1.00 1.00 re f My question is why do you need so many points? Have you looked at 'hexbin' as way of plotting the data? Have you considered generating 'jpg' output which will be much smaller; e.g. a jpg file with 10,000 points was 70K and one with 100,000 was 90K -- it only created the pixels to be plotted. On Mon, Sep 7, 2009 at 7:43 AM, Christian Ritterchristian.rit...@uclouvain.be wrote: Hi all, I have to produce arrangements of 25 simple plots of the type plot(x,y,pch=.) where there are typically on the order of 2 points. So, overall, I have about 50 points. When I use the pdf device, I get file sizes (on a Windows machine) of about 10 MB. When I then zip the files, I'm down to about 0.5MB, so the original pdf files were created with a lot of 'air'. I'm wondering why the pdf files are so large and whether there is an alternative to produce them smaller. Any ideas? Have a nice afternoon, Chris. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Jim Holtman Cincinnati, OH +1 513 646 9390 What is the problem that you are trying to solve? __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Equivalence of Mann-Whitney test and Kruskal-Wallis test with k=2
Hi all, The Kruskal-Wallis test is a generalization of the two-sample Mann-Whitney test to *k* samples. That being the case, the Kruskal-Wallis test with *k*=2 should give an identical p-value to the Mann-Whitney test, should it not? x1-c(1:5) x2-c(6,8,9,11) a-wilcox.test(x1,x2,paired=FALSE) b-kruskal.test(list(x1,x2),paired=FALSE) a$p.value [1] 0.01587302 b$p.value [1] 0.01430588 The p-values are slightly different (note that there are no ties in the data, so computed p-values should be exact). Can anyone explain the discrepancy? It's been awhile since I studied nonparametric stats and this one has me scratching my head. Many thanks! Tom [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Size of plots in pdf files#can it be smaller?
Hi,
I've got a trial version of a thinScatter() function that
(down-)samples 2d-scatter plots while preserving the empirical density
distribution. You can grab it by:
source(http://www.braju.com/R/hbLite.R;);
hbLite(scatterPlots);
Example from example(thinScatter):
library(scatterPlots);
# Sample scatter data
n - 10e3
x - rnorm(n=n)
y - rnorm(n=n)
xy - cbind(x=x, y=sin(x)+y/5)
layout(matrix(1:4, nrow=2, byrow=TRUE))
par(mar=c(5,4,2,1))
# Plot data
plot(xy, pch=1)
# Thin scatter data by subsampling
rhos - c(1/3, 1/4, 1/6)
for (kk in seq(along=rhos)) {
xy2 - thinScatter(xy, size=rhos[kk])
points(xy2, pch=1, col=kk+1)
title - paste(sprintf(%.1f%%, 100*c(1, rhos)), collapse=, )
mtext(side=3, line=0, paste(Densities:, title))
}
for (kk in seq(along=rhos)) {
xy2 - thinScatter(xy, size=rhos[kk])
plot(xy2, pch=1, col=kk+1)
mtext(side=3, line=0, sprintf(Density: %.1f%% [n=%d],
100*rhos[kk], nrow(xy2)))
}
If the mailing list allows it, a PNG is attached.
/Henrik
On Mon, Sep 7, 2009 at 1:50 PM, Philipp Pagelp.pa...@wzw.tum.de wrote:
I have to produce arrangements of 25 simple plots of the type
plot(x,y,pch=.) where there are typically on the order of 2
points.
So, overall, I have about 50 points. When I use the pdf device,
I get file sizes (on a Windows machine) of about 10 MB.
When I then zip the files, I'm down to about 0.5MB, so the original
pdf files were created with a lot of 'air'.
I'm wondering why the pdf files are so large and whether there is an
alternative to produce them smaller. Any ideas?
As far as I know, R cannot directly produce compressed pdf files. You
could postprocess your plots. E.g. pdftk is quite conveniant:
pdftk foo.pdf cat output foo2.pdf compress
Maybe in a loop if many files are involved.
cu
Philipp
--
Dr. Philipp Pagel
Lehrstuhl für Genomorientierte Bioinformatik
Technische Universität München
Wissenschaftszentrum Weihenstephan
85350 Freising, Germany
http://webclu.bio.wzw.tum.de/~pagel/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
attachment: thinScatter,Rex.png__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Using anova(f1, f2) to compare lmer models yields seemingly erroneous Chisq = 0, p = 1
Thank you all for your insight! I am glad to hear, at least, that I am doing something incorrectly (since the results do not make sense), and I am very grateful for your attempts to remedy my very limited (and admittedly self-taught) understanding of multilevel models and R. As I mentioned in the problem statement, predictor.1 explains vastly more variance in outcome than predictor.2 (R2 = 15% vs. 5% in OLS regression, with very large N), and the model estimates are very similar for the multilevel model as for OLS regression. Therefore, I am quite confident that predictor.1 comprises a much better model. I understand that several of you are saying that anova() cannot be used to compare these two multilevel models. Is there *any* way to compare two predictors to see which better predicts the outcome in a multilevel model? f1's lower AIC and BIC, and higher logLik are concordant with the idea that predictor.1 is superior to predictor.2, as best as I understand it, but is there any way to test whether that difference is statistically significant? The only function I can find online is anova() to compare models, but its output is nonsensical and, as you are all saying, it does not apply to my situation anyway. Interestingly, anova() seems to work if I arbitrarily subset my observations, but when I use all the observations anova() generates Chisq = 0. This is probably a red herring but I thought I would mention it in case it is not. Also, I concede that I am confused what you mean that the two models (f1 and f2) are not nested, and therefore anova() cannot be used. What would be an example of a nested model: comparing predictor.1 to a model with both predictor.1 and predictor.2? Surely there must also be a way to compare the predictive power of predictor.1 and predictor.2 to each other in a zero-order sense, but I am at a loss to identify it. Alain Zuur wrote: rapton wrote: When I run two models, the output of each model is generated correctly as far as I can tell (e.g. summary(f1) and summary(f2) for the multilevel model output look perfectly reasonable), and in this case (see below) predictor.1 explains vastly more variance in outcome than predictor.2 (R2 = 15% vs. 5% in OLS regression, with very large N). What I am utterly puzzled by is that when I run an anova comparing the two multilevel model fits, the Chisq comes back as 0, with p = 1. I am pretty sure that fit #1 (f1) is a much better predictor of the outcome than f2, which is reflected in the AIC, BIC , and logLik values. Why might anova be giving me this curious output? How can I fix it? I am sure I am making a dumb error somewhere, but I cannot figure out what it is. Any help or suggestions would be greatly appreciated! -Matt f1 - (lmer(outcome ~ predictor.1 + (1 | person), data=i)) f2 - (lmer(outcome ~ predictor.2 + (1 | person), data=i)) anova(f1, f2) Data: i Models: f1: outcome ~ predictor.1 + (1 | person) f2: outcome ~ predictor.2 + (1 | person) DfAIC BIClogLik Chisq Chi Df Pr(Chisq) f1 6 45443 45489 -22715 f2 25 47317 47511 -23633 0 19 1 Your models are nest nestedit doesn't make sense to do. Alain -- View this message in context: http://www.nabble.com/Using-anova%28f1%2C-f2%29-to-compare-lmer-models-yields-seemingly-erroneous-Chisq-%3D-0%2C-p-%3D-1-tp25297254p25338046.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Very inaccurate circles
Hello,
Please take a look at the attached plot and let me know if this is
normal. The circle has radio I am using R 2.9.2 inside OS X Leopard.
The plot was generated with:
png('bizarre_circle.png')
plot(c(-5,0,0,5), c(0,5,-5,0))
symbols(0,0, circles=c(sqrt(25)), inches=FALSE, add=TRUE)
dev.off()
If I resize the x11's plot window I can manage to make the circle pass
through all the points but it's the first time I have to do this. I
understand this may be happening because of the method to draw
circles. I googled and found this other methods:
http://tolstoy.newcastle.edu.au/R/help/03a/3364.html
But method 2 and method 3 seem to use different units than my points
and working with them is not so straightforward.
1. Is there a way to draw circles in the same unit system than plot
already uses?
2. Is there a way of making method 3 (which uses mm) work with the
same unit system as the plot?
Thank you,
--
slnc
attachment: bizarre_circle.png__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Very inaccurate circles
On 8/09/2009, at 11:58 AM, Juan Alonso wrote:
Hello,
Please take a look at the attached plot and let me know if this is
normal. The circle has radio I am using R 2.9.2 inside OS X Leopard.
The plot was generated with:
png('bizarre_circle.png')
plot(c(-5,0,0,5), c(0,5,-5,0))
symbols(0,0, circles=c(sqrt(25)), inches=FALSE, add=TRUE)
dev.off()
snip
Please don't post implied criticism of software (``very inaccurate'')
when the fault lies not in the software but in your understanding.
There is *nothing* inaccurate about the circles. The symbols() function
is cleverly designed to produce ``circles'' that ***look*** like circles
when plotted. In most cases they ***aren't*** circles, but rather
ellipses,
with eccentricity adjusted to compensate for the aspect ratio of the
plot
on which they are being superimposed.
If you want the symbols()-created circle to pass through points which
really
do lie on a circle you need to make the aspect ratio of the plot
equal to 1:
plot(c(-5,0,0,5), c(0,5,-5,0),asp=1)
symbols(0,0, circles=c(sqrt(25)), inches=FALSE, add=TRUE)
OM!
cheers,
Rolf Turner
##
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] using an array of strings with strsplit, issue when including a space in split criteria
I get a different result: txt - c(sales to 23 August 2008 published 29 August,sales to 6 September 2008 published 11 September) strsplit(txt, 'published ', fixed=TRUE) [[1]] [1] sales to 23 August 2008 29 August [[2]] [1] sales to 6 September 2008 11 September sessionInfo() R version 2.9.0 (2009-04-17) i386-pc-mingw32 locale: LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base On Mon, Sep 7, 2009 at 11:40 AM, Tony Breyaltony.bre...@googlemail.com wrote: Dear all, I'm having a problem understanding why a split does not occur with in the 2nd use of the function strsplit below: # text strings txt - c(sales to 23 August 2008 published 29 August, + sales to 6 September 2008 published 11 September) # first use strsplit(txt, 'published', fixed=TRUE) [[1]] [1] sales to 23 August 2008 29 August [[2]] [1] sales to 6 September 2008 11 September # second use, but with a space ' ' in the split strsplit(txt, 'published ', fixed=TRUE) [[1]] [1] sales to 23 August 2008 29 August [[2]] [1] sales to 6 September 2008 published 11 September Thank you kindly for any help in advance. Tony O/S: Win Vista Ultimate sessionInfo() R version 2.9.2 (2009-08-24) i386-pc-mingw32 locale: LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United Kingdom. 1252;LC_MONETARY=English_United Kingdom. 1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] RODBC_1.3-0 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Omnibus test for main effects in the face of aninteraction containing the main effects.
John, your question is confusing. After reading it twice, I still cannot
figure out what exactly you want to compare.
Your model a is the unrestricted model, and model b is a restricted
version of model a (i.e., b is a hiearchically reduced version of a, or
put differently, all coefficients of b are in a with a having additional
coefficients). Thus, it is appropriate to compare the models (also called
nested models).
Comparing c with a and d with a is also appropriate for the same reason.
However, note that depedent on discipline, it may be highly unconventional
to fit an interaction without all direct effects of the interacted variables
(the reason for this being that you may get biased estimates).
What you might consider is:
1. Run an intercept only model
2. Run a model with group and time
3. Run a model with group, time, and the interaction
Then compare 2 to 1, and 3 to 2. This tells you whether including more
variables (hierarchically) makes your model better.
HTH,
Daniel
On a different note, if lme fits with restricted maximum likelihood, I
think I remember that you cannot compare them. You have to fit them with
maximum likelihood. I am pointing this out because lmer with restricted
maximum likelihood by standard, so lme might too.
-
cuncta stricte discussurus
-
-Ursprüngliche Nachricht-
Von: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] Im
Auftrag von John Sorkin
Gesendet: Monday, September 07, 2009 4:00 PM
An: r-help@r-project.org
Betreff: [R] Omnibus test for main effects in the face of aninteraction
containing the main effects.
R 2.9.1
Windows XP
UPDATE,
Even my first suggestion
anova(fita,fitb) is probably not appropriate as the fixed effects are
different in the two model, so I don't even know how to perform the ombnibus
test for the interaction!
I am fitting a random effects ANOVA with two factors Group which has two
levels and Time which has three levels:
fita-lme(Post~Time+factor(Group)+factor(Group)*Time,
random=~1|SS,data=blah$alldata)
I want to get the omnibus significance tests for each factor and the
interaction. I believe I can get the omnibus test for the interaction by
running the model:
fitb-lme(Post~Time+factor(Group), random=~1|SS,data=blah$alldata) followed
by anova(fita,fitb).
How do I get the omnibus test for the main effects i.e. for Time and
factor(Group)? I could drop each from the model, i.e.
fitc-lme(Post~ factor(Group)+factor(Group)*Time,
random=~1|SS,data=blah$alldata)
fitd-lme(Post~Time+factor(Group)*Time,
random=~1|SS,data=blah$alldata)
and then run
anova(fita,fitc)
anova(fita,fitd)
but I don't like this option as it will have in interaction that contains a
factor that is not included in the model as a main effect. How then do I get
the omnibus test for Time and factor(Group)?
Thanks
John
John David Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
Confidentiality Statement:
This email message, including any attachments, is for th...{{dropped:9}}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] OT - Banker's Algorithum
Hi,
On Mon, Sep 7, 2009 at 5:41 AM, Stephen Liusati...@yahoo.com wrote:
Hi folks,
Can R-Project be used to perform Banker's Algorithum?
Yes.
http://en.wikipedia.org/wiki/Banker%27s_algorithm
You can pretty easily directly transalte the pseudocde listed in the
section below to R:
http://en.wikipedia.org/wiki/Banker%27s_algorithm#Pseudo-Code.5B3.5D
Not really sure what else to say ... for instance, the line:
P = P - {p}
Might be:
P - setdiff(P, p)
Maybe use lists to hold the stuff in P? Don't know what else you're
looking for ...
-steve
--
Steve Lianoglou
Graduate Student: Computational Systems Biology
| Memorial Sloan-Kettering Cancer Center
| Weill Medical College of Cornell University
Contact Info: http://cbio.mskcc.org/~lianos/contact
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Omnibus test for main effects in the face ofaninteraction containing the main effects.
Daniel, When Group is entered as a factor, and the factor has two levels, the ANOVA table gives a p value for each level of the factor. What I am looking for is the omnibus p value for the factor, i.e. the test that the factor (with all its levels) improves the prediction of the outcome. You are correct that normally one could rely on the fact that the model Post-Time+as.factor(Group)+as.factor(Group)*Time contains the model Post-Time+as.factor(Group) and compare the two models using anova(model1,model2). However, my model is has a random effect, the comparison is not so easy. The REML comparions of nested random effects models is not valid when the fixed effects are not the same in the models, which is the essence of the problem in my case. In addition to the REML problem if one wants to perform an omnibus test for Group, one would want to compare nested models, one containing Group, and the other not containing group. This would suggest comparing Post-Time+ as.factor(Group)*Time to Post-Time+Group+as.factor(Group)*Time The quandry here is whether one should or not allow the first model as it is poorly specified - one term of the interaction, as.factor(Group)*Time, as.factor(Group) does not appear as a main effect - a no-no in model building. John John Sorkin M.D., Ph.D. Chief, Biostatistics and Informatics Baltimore VA Medical Center GRECC, University of Maryland School of Medicine Claude D. Pepper OAIC, University of Maryland Clinical Nutrition Research Unit, and Baltimore VA Center Stroke of Excellence University of Maryland School of Medicine Division of Gerontology Baltimore VA Medical Center 10 North Greene Street GRECC (BT/18/GR) Baltimore, MD 21201-1524 (Phone) 410-605-7119 (Fax) 410-605-7913 (Please call phone number above prior to faxing) jsor...@grecc.umaryland.edu Daniel Malter dan...@umd.edu 09/07/09 9:23 PM John, your question is confusing. After reading it twice, I still cannot figure out what exactly you want to compare. Your model a is the unrestricted model, and model b is a restricted version of model a (i.e., b is a hiearchically reduced version of a, or put differently, all coefficients of b are in a with a having additional coefficients). Thus, it is appropriate to compare the models (also called nested models). Comparing c with a and d with a is also appropriate for the same reason. However, note that depedent on discipline, it may be highly unconventional to fit an interaction without all direct effects of the interacted variables (the reason for this being that you may get biased estimates). What you might consider is: 1. Run an intercept only model 2. Run a model with group and time 3. Run a model with group, time, and the interaction Then compare 2 to 1, and 3 to 2. This tells you whether including more variables (hierarchically) makes your model better. HTH, Daniel On a different note, if lme fits with restricted maximum likelihood, I think I remember that you cannot compare them. You have to fit them with maximum likelihood. I am pointing this out because lmer with restricted maximum likelihood by standard, so lme might too. - cuncta stricte discussurus - -Ursprüngliche Nachricht- Von: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] Im Auftrag von John Sorkin Gesendet: Monday, September 07, 2009 4:00 PM An: r-help@r-project.org Betreff: [R] Omnibus test for main effects in the face of aninteraction containing the main effects. R 2.9.1 Windows XP UPDATE, Even my first suggestion anova(fita,fitb) is probably not appropriate as the fixed effects are different in the two model, so I don't even know how to perform the ombnibus test for the interaction! I am fitting a random effects ANOVA with two factors Group which has two levels and Time which has three levels: fita-lme(Post~Time+factor(Group)+factor(Group)*Time, random=~1|SS,data=blah$alldata) I wantinteraction. I believe I can get the omnibus test for the interaction by running the model: fitb-lme(Post~Time+factor(Group), random=~1|SS,data=blah$alldata) followed by anova(fita,fitb). How do I get the omnibus test for the main effects i.e. for Time and factor(Group)? I could drop each from the model, i.e. fitc-lme(Post~ factor(Group)+factor(Group)*Time, random=~1|SS,data=blah$alldata) fitd-lme(Post~Time+factor(Group)*Time, random=~1|SS,data=blah$alldata) and then run anova(fita,fitc) anova(fita,fitd) but I don't like this option as it will have in interaction that contains a factor that is not included in the model as a main effect. How then do I get the omnibus test for Time and factor(Group)? Thanks John John David Sorkin M.D., Ph.D. Chief, Biostatistics and Informatics University of Maryland School of Medicine Division of Gerontology Baltimore VA Medical Center 10 North Greene Street GRECC (BT/18/GR) Baltimore, MD 21201-1524 (Phone) 410-605-7119 (Fax) 410-605-7913 (Please call
Re: [R] Omnibus test for main effects in the face ofaninteraction containing the main effects.
John Sorkin wrote: Daniel, When Group is entered as a factor, and the factor has two levels, the ANOVA table gives a p value for each level of the factor. What I am looking for is the omnibus p value for the factor, i.e. the test that the factor (with all its levels) improves the prediction of the outcome. You are correct that normally one could rely on the fact that the model Post-Time+as.factor(Group)+as.factor(Group)*Time contains the model Post-Time+as.factor(Group) and compare the two models using anova(model1,model2). However, my model is has a random effect, the comparison is not so easy. The REML compari[s]ons of nested random effects models is not valid when the fixed effects are not the same in the models, which is the essence of the problem in my case. So why not use ML comparisons if this is what you want to do ... ? In addition to the REML problem if one wants to perform an omnibus test for Group, one would want to compare nested models, one containing Group, and the other not containing group. This would suggest comparing Post-Time+ as.factor(Group)*Time to Post-Time+Group+as.factor(Group)*Time The quand[a]ry here is whether one should or not allow the first model as it is poorly specified - one term of the interaction, as.factor(Group)*Time, as.factor(Group) does not appear as a main effect - a no-no in model building. John Agreed. I don't think the question makes sense -- the overall effect of group includes both the main effect the interaction. (re) reading Venables, W. N. “Exegeses on Linear Models.” 1998 International S-PLUS User Conference. Washington, DC, 1998. http://www.stats.ox.ac.uk/pub/MASS3/Exegeses.pdf. might be useful ... -- View this message in context: http://www.nabble.com/Re%3A-Omnibus-test-for-main-effects-in-the-face%09ofaninteraction%09containing-the-main-effects.-tp25338728p25338952.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Help with use of rep function in R
Dear List,I am trying to use rep function in the following conditions A = c( 5, 6, 7, 11, 9, 12, 10, 15) B = c(12,15, 21, 31, 25, 27,32, *34*,13,12, 34, 33, 24, 29, 26, *28*,22,14,27,22,21,12,32, 16) I need to repeat each element of A, as many times as each element of B, for the entire length of B. for example, repeat 5, for 12 times, 6 for 15 times,, 15 for 34 times, and then, again, 5 for 13 times, 6 for 12 times,, 15 for 28 times, and so on. I used, the function rep(A, times = B) It didn't work because apparently, the times command , worked for only if the length of A and B was equal. Thank you for your help in advance. -- Subodh Acharya University of Florida Gainesville, FL. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] using an array of strings with strsplit, issue when including a space in split criteria
I am using the exact same version of R as you also on Vista but can't reproduce your result. For me it splits properly. Try starting R like this (modify path if needed) from the Windows cmd line: \Program Files\R\R-2.9.2\bin\Rgui --vanilla and then try it. On Mon, Sep 7, 2009 at 11:40 AM, Tony Breyaltony.bre...@googlemail.com wrote: Dear all, I'm having a problem understanding why a split does not occur with in the 2nd use of the function strsplit below: # text strings txt - c(sales to 23 August 2008 published 29 August, + sales to 6 September 2008 published 11 September) # first use strsplit(txt, 'published', fixed=TRUE) [[1]] [1] sales to 23 August 2008 29 August [[2]] [1] sales to 6 September 2008 11 September # second use, but with a space ' ' in the split strsplit(txt, 'published ', fixed=TRUE) [[1]] [1] sales to 23 August 2008 29 August [[2]] [1] sales to 6 September 2008 published 11 September Thank you kindly for any help in advance. Tony O/S: Win Vista Ultimate sessionInfo() R version 2.9.2 (2009-08-24) i386-pc-mingw32 locale: LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United Kingdom. 1252;LC_MONETARY=English_United Kingdom. 1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] RODBC_1.3-0 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Help with use of rep function in R
On 8/09/2009, at 2:20 PM, Subodh Acharya wrote:
Dear List,I am trying to use rep function in the following conditions
A = c( 5, 6, 7, 11, 9, 12, 10, 15)
B = c(12,15, 21, 31, 25, 27,32, *34*,13,12, 34, 33, 24, 29, 26,
*28*,22,14,27,22,21,12,32,
16)
I need to repeat each element of A, as many times as each element
of B, for
the entire length of B.
for example,
repeat 5, for 12 times, 6 for 15 times,, 15 for 34 times,
and then,
again, 5 for 13 times, 6 for 12 times,, 15 for 28 times,
and so on.
I used, the function
rep(A, times = B)
It didn't work because apparently, the times command , worked for
only if
the length of A and B was equal.
A - c( 5, 6, 7, 11, 9, 12, 10, 15)
B - c(12,15, 21, 31, 25, 27,32, 34,13,12, 34, 33, 24, 29, 26,
28,22,14,27,22,21,12,32,16)
AA - rep(A,length.out=length(B))
RA - rep(AA,times=B)
cheers,
Rolf Turner
##
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Very inaccurate circles
Try this:
library(plotrix)
plot(c(-5,0,0,5), c(0,5,-5,0), asp = 1)
draw.circle(0, 0, 5)
On Mon, Sep 7, 2009 at 7:58 PM, Juan Alonsos...@slnc.me wrote:
Hello,
Please take a look at the attached plot and let me know if this is
normal. The circle has radio I am using R 2.9.2 inside OS X Leopard.
The plot was generated with:
png('bizarre_circle.png')
plot(c(-5,0,0,5), c(0,5,-5,0))
symbols(0,0, circles=c(sqrt(25)), inches=FALSE, add=TRUE)
dev.off()
If I resize the x11's plot window I can manage to make the circle pass
through all the points but it's the first time I have to do this. I
understand this may be happening because of the method to draw
circles. I googled and found this other methods:
http://tolstoy.newcastle.edu.au/R/help/03a/3364.html
But method 2 and method 3 seem to use different units than my points
and working with them is not so straightforward.
1. Is there a way to draw circles in the same unit system than plot
already uses?
2. Is there a way of making method 3 (which uses mm) work with the
same unit system as the plot?
Thank you,
--
slnc
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Help with use of rep function in R
Dear List,I am trying to use rep function in the following conditions A = c( 5, 6, 7, 11, 9, 12, 10, 15) B = c(12,15, 21, 31, 25, 27,32, *34*,13,12, 34, 33, 24, 29, 26, *28*,22,14,27,22,21,12,32, 16) I need to repeat each element of A, as many times as each element of B, for the entire length of B. for example, repeat 5, for 12 times, 6 for 15 times,, 15 for 34 times, and then, again, 5 for 13 times, 6 for 12 times,, 15 for 28 times, and so on. I used, the function rep(A, times = B) It didn't work because apparently, the times command , worked for only if the length of A and B was equal. Thank you for your help in advance. -- Subodh Acharya University of Florida Gainesville, FL. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
