Re: [R] finding x values to meet a y
Antonio Olinto wrote: Hi, I'd like to find which values of x will give me a y. In other words, in the example of a gaussian curve, I want to find the values of x that will give me a density, let's say, of 0.02. curve(((1/(sqrt(2*pi)*10))*exp(-((x-50)^2)/(2*10^2))),xlim=c(0,100)) Numerical optimization might help: optimize( function(x) (1/(sqrt(2*pi)*10) * exp(-((x-50)^2)/(2*10^2)) - 0.02)^2, interval = c(0, 50)) Uwe Ligges Thanks for any help, Antonio Olinto - WebMail Bignet - O seu provedor do litoral www.bignet.com.br __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] geodesic distance
I don't know exactly, but you can have a look at the sna package,
there should be a function for geodesic computation in social network
analysis...
regards,
Simone
Il giorno 03/ago/06, alle ore 11:14, stefano iacus ha scritto:
Hi,
has anyone ever seen implemented in R the following geodesic
distance between positive definite pxp matrices A and B?
d(A,B) = \sum_{i=1}^p (\log \lambda_i)^2
were \lambda is the solution of det(A -\lambda B) = 0
thanks
stefano
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-
guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] integrate() problem {was mathematica - r ...}
Leo == Leo G?rtler [EMAIL PROTECTED]
on Tue, 08 Aug 2006 00:13:19 +0200 writes:
Leo Dear R-list,
Leo I try to transform a mathematica script to R.
Leo ###relevant part of the Mathematica script
Leo (* p_sv *)
Leo dd = NN (DsD - DD^2);
Leo lownum = NN (L-DD)^2;
Leo upnum = NN (H-DD)^2;
Leo low = lownum/(2s^2);
Leo up = upnum/(2s^2);
Leo psv = NIntegrate[1/(s^NN) Exp[-dd/(2s^2)]
Leo (Gamma[1/2,0,up] + Gamma[1/2,0,low]),{s,sL,sH}, MinRecursion- 3];
Leo PSV = psv/Sqrt[2NN];
Leo Print[- Results ];
Leo Print[ ];
Leo Print[p(sv|D_1D_2I) = const. ,N[PSV,6]];
Leo
Leo # R part
Leo library(fOptions)
Leo ###raw values for reproduction
Leo NN - 58
Leo dd - 0.411769
Leo lownum - 20.81512
Leo upnum - 6.741643
Leo sL - 0.029
Leo sH - 0.092
Leo ###
Leo integpsv - function(s) { 1 / (s^NN) * exp(-dd / (2 * s^2)) *
Leo( (igamma((upnum/(2*s^2)),1/2) - igamma(0,1/2) ) +
Leo(igamma((lownum/(2*s^2)),1/2) - igamma(0,1/2) ) )
Leo }
Leo psv - integrate(integpsv, lower=sL, upper=sH)
Leo PSV - psv$value / sqrt(2*NN)
Leo print(- Results \n)
Leo print(paste(p(sv|D_1D_2I) = const. ,PSV, sep=))
Leo The results of variable PSV are not the same.
Leo In mathematica - PSV ~ 2.67223e+47
Leo with rounding errors due to the initial values, in R - PSV ~ 1.5e+47
Leo I am not that familiar with gamma functions and integration, thus I
Leo assume there the source of the problem can be located.
Yes.
A few remarks
1) No need to use package fOptions and igamma();
standard R's pgamma() is all you need
{igamma() has added value only for *complex* arguments!}
2) igamma(0, 1/2) == pgamma(0, 1/2) == 0 , so you can really
drop them from your integrand.
integpsv - function(s) {
1 / (s^NN) * exp(-dd / (2 * s^2)) *
( pgamma(upnum/(2*s^2), 1/2) + pgamma(lownum/(2*s^2), 1/2) )
}
3) But then the problem could really be with the algorithm used in
integrate(), and indeed if you plot your integrand
plot(integpsv, from= sL, to = sH)
you see that indeed your integrand looks ``almost
constant'' in the left half --- whereas that is actually not
true but the range of the integrand varies so dramatically
that it ``looks like'' constant 0 upto about x= .06.
Something which help(integrate) warns about.
However, if I experiment, using integrate() in two parts, or using many other
numerical integration approximators,
all methods give ( your 'psv', not PSV )
integrate(integpsv, lower=sL, upper=sH)
a value of 1.623779e+48 (which leads to your PSV of 1.5076e+47)
Could it be that you are not using the same definition of
incomplete gamma in Mathematica and R ?
Martin Maechler, ETH Zurich
Leo Thanks for helping me to adjust the sript.
Leo best wishes
Leo leo
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] finding x values to meet a y
AOlinto == Antonio Olinto [EMAIL PROTECTED] on Mon, 7 Aug 2006 19:49:42 -0300 writes: AOlinto Hi, AOlinto I'd like to find which values of x will give me a y. AOlinto In other words, in the example of a gaussian curve, I want to find the values of AOlinto x that will give me a density, let's say, of 0.02. AOlinto curve(((1/(sqrt(2*pi)*10))*exp(-((x-50)^2)/(2*10^2))),xlim=c(0,100)) In other words, you are trying to *invert a function* numerically, i.e., find x such that f(x) = y0 Now, this is trivially the same as finding the root (or zero) of the function f~(x) = f(x) - y0 == Use *the* R root finding function : uniroot() [or polyroot() if f(.) is a polynomial] Martin Maechler, ETH Zurich __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Constrain coefs. in linear model to sum to 0
BillV == [EMAIL PROTECTED] on Tue, 8 Aug 2006 15:52:20 +1000 writes: BillV Gorjanc Gregor asks: Hello! I would like to use constrain to sum coeficients of a factor to 0 instead of classical corner contraint i.e. I would like to fit a model like lm(y ~ 1 + effectA + effectB) and say get parameters intercept effectA_1 effectA_2 effectB_1 effectB_2 effectB_3 where effectA_1 represents deviation of level A_1 from intercept and sum(effectA_1, effectA_2) = 0 and the same for factor B. Is this possible to do? BillV Yes, but not quite as simply as you would like. If you set BillV options(contrasts = c(contr.sum, contr.poly)) BillV for example, then factor models are parametrised as BillV you wish above, BUT you don't get all the effects BillV directly BillV In your case above, for example, if fm is the fitted BillV model object, then BillV coef(fm) BillV Will give you intercept, effectA_2, effectB_2, BillV effectB_3. The remaining effects*_1 you will need to BillV calculate as the negative of the sum of all the BillV others. BillV This gets a bit more complicated when you have BillV crossed terms, a*b, but the same principle applies. Further, note that there have been functions dummy.coef() and model.tables() that have been intended to do exactly this. Unfortunately, these two functions only work correctly in simpler cases (e.g., complete balanced designs, IIRC) E.g. help(model.tables) has Warning: The implementation is incomplete, and only the simpler cases have been tested thoroughly. and help(dummy.coef) Details: [] The method used has some limitations, and will give incomplete results for terms such as 'poly(x, 2))'. However, it is adequate for its main purpose, 'aov' models. [] Warning: This function is intended for human inspection of the output: it should not be used for calculations. Use coded variables for all calculations. The results differ from S for singular values, where S can be incorrect. When we last discussed this, IIRC, we (some within R-core) were waiting for interested (and knowledgable) users to provide improved code for these functions. Hey, if you want to become immortal in the R annals, now is your chance! ;-) BillV Bill Venables. Lep pozdrav / With regards, Gregor Gorjanc [.] -- One must learn by doing the thing; for though you think you know it, you have no certainty until you try. Sophocles ~ 450 B.C. BillV Well, now's your chance! indeed! Martin Maechler, ETH Zurich __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] OT: RBanking
Hi all, This might be a little off-topic, but I wanted to let you all know about a project for which I used R. I started an (free for all, open source) offline banking tool, programmed in R and using the tcltk library for the user interface. It has several nice features already :) It can use regular expressions to match/group bank transactions and can obviously use R's power to plot nice graphs. Personally, I found it too difficult to make nice graphs with a program like gnucash (which seems more for experts instead of home-users??) There is support for several Dutch banks at the moment (I don't know about internationalization). The program is not yet user-friendly and not trivial to install, but if you are interested I invite you to try it. It can be checked out with the command: svn co https://svn.sourceforge.net/svnroot/rbanking/trunk rbanking And the website is at http://rbanking.sourceforge.net (thereafter you would need to install several dependencies like BWidgets and tablelist. It would be great if there's a group of people who would slowly like to further develop this, like myself. Any help is appreciated :) Jonne. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Source installation error: gfortran and gcc disagree on int and double ...
First, you replied to the list and not to me, which was discourteous. Your error does indeed appear to be an incorrectly installed compiler. conftestf.o: In function `cftest_': /home/mike/Desktop/R-2.3.1/conftestf.f:9: undefined reference to `_gfortran_pow_r8_i4' this is in -lgfortran, so that is not being found or is broken. One possibility is a missing symlink that is in a -devel RPM. This is not an R issue. On Mon, 7 Aug 2006, Mike wrote: Prof Brian Ripley wrote: First, you do not need -fPIC: it is the same as -fpic which R selects on your platform. Right, I commented the flags and configure goes just as far Second, please look at config.log to find the exact problem: it well be that your compilers are not properly installed (as was the case in the reference you quote below). from config.log it seem that there are a few missing inclides and syntax errors in confdefs.h. (I can quote the specifics, if necessary). The last lines before the error message are: configure:27770: checking whether gfortran and gcc agree on int and double conftest.c: In function 'main': conftest.c:28: warning: implicit declaration of function 'printf' conftest.c:28: warning: incompatible implicit declaration of built-in function 'printf' conftest.c:29: warning: implicit declaration of function 'exit' conftest.c:29: warning: incompatible implicit declaration of built-in function 'exit' conftestf.o: In function `cftest_': /home/mike/Desktop/R-2.3.1/conftestf.f:9: undefined reference to `_gfortran_pow_r8_i4' collect2: ld returned 1 exit status configure:27848: WARNING: gfortran and gcc disagree on int and double configure:27850: error: Maybe change CFLAGS or FFLAGS? Finally, your compilers are pretty obselete (there have been 4.0.2, 4.0.3, 4.1.0 and 4.1.1), so you should be updating them. Changing compiler is not an option for me. Could you advise me how to configure the flags? -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Netiquette, was Re: ... gfortran and gcc...
Prof Brian Ripley [EMAIL PROTECTED] writes: First, you replied to the list and not to me, which was discourteous. You mean that he replied to the list *only*, I hope. I usually consider it offensive when people reply to me and not the list (reasons including: It feels like being grabbed by the sleeve, I might not actually be the best source for the answer, and it's withholding the answer from the rest of the subscribers.) -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Netiquette, was Re: ... gfortran and gcc...
On Tue, 8 Aug 2006, Peter Dalgaard wrote: Prof Brian Ripley [EMAIL PROTECTED] writes: First, you replied to the list and not to me, which was discourteous. You mean that he replied to the list *only*, I hope. Yes, and it was written as if to me, and was a reply to an email from me. I usually consider it offensive when people reply to me and not the list (reasons including: It feels like being grabbed by the sleeve, I might not actually be the best source for the answer, and it's withholding the answer from the rest of the subscribers.) We do ask people to copy to the list. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] fixed effects following lmer and mcmcsamp - which to present?
Dear all, I am running a mixed model using lmer. In order to obtain CI of individual coefficients I use mcmcsamp. However, I need advice which values that are most appropriate to present in result section of a paper. I have not used mixed models and lmer so much before so my question is probably very naive. However, to avoid to much problems with journal editors and referees addicted to p-values, I would appreciate advice of which values of the output for the fixed factor that would be most appropriate to present in a result section, in order to convince them of the p-value free 'lmer-mcmcsamp'-approach! Using the example from the help page on lmer: fm1 - lmer(Reaction ~ Days + (Days|Subject), sleepstudy) ...I obtain the following for 'Days': summary(fm1) ... Estimate Std. Error t value Days 10.4673 1.54586.771 ...and from mcmcsamp: summary(mcmcsamp(fm1 , n = 1)) 1. Empirical mean and standard deviation for each variable, plus standard error of the mean: Mean SD Naive SE Time-series SE Days 10.4695 1.7354 0.017354 0.015921 2. Quantiles for each variable: 2.5% 25% 50% 75%97.5% Days 7.02279.3395 10.4712 11.5719 13.957 The standard way of presenting coefficients following a 'non-lmer' output is often (beta=..., SE=..., statistic=..., P=...). What would be the best equivalent in a 'lmer-mcmcsamp-context'? (beta=..., CI=...) is a good start I believe. But which beta? And what else? I assume that the a 95% CI in this case would be 7.0227-13.957 (please, do correct me I have completely misunderstood!). But which would be the corresponding beta? 10.4673?, 10.4695? 10.4712? Is the t-value worth presenting or is it 'useless' without corresponding degrees of freedom and P-value? If I present the mcmcsamp-CI, does it make sense to present any of the three SE obtained in the output above? BTW, I have no idea what Naive SE, Time-series SE means. Could not find much in help and pdfs to coda or Matrix, or in Google. Thanks in advance for any advice and hints to help-texts I have missed! Best regards, Henrik __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] gamm question
Hello, I have two gamm question (I am using gamm in mgcv). 1. In have, say 5 time series. Monthly data, 20 year. The 5 time series are from 5 stations. The data are in vectors, so I have fitted something along the lines of: tmp-gamm(Y ~ s(Year,by=station1)+s(Year,by=station2)+ s(Year,by=station3)+s(Year,by=station4)+ s(Year,by=station5)+ factor(station)*factor(month), correlation=corAR1(form=~MyTime|station), famliy=gaussian) station is just a long vector with ones, twos, threes, fours and fives. MyTime defines the order of time (and has values 1 to 240) This model fits a Year smoother on each station and, and has one auto-regressive parameter (whihc is about 0.3). How woul I allow for 5 different AR1 parameters (one per station)? So far so good. It runs..and get the output. The problem is that the errors from station 1 are correlated with those of station 2 (as they are close in space). Same holds for other stations. The cross-correlation is about 0.5. How do I build in this between-station correlation? So..I have within station auto-correlation (dealt by the AR1 parameter), and between station correlation. Question 2: Why is Simon Wood using in his 2006 book only ML estimation and not REML? I thought that REML was used to compare models with different random components and ML estimation models with different fixed components? Kind regards, Piet Bell - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] fixed effects constant in mcmcsamp
I'm fitting a GLMM to some questionnaire data. The structure is J individuals,
nested within I areas, all of whom answer the same K (ordinal) questions. The
model I'm using is based on so-called continuation ratios, so that it can be
fitted using the lme4 package.
The lmer function fits the model just fine, but using mcmcsamp to judge the
variability of the parameter estimates produces some strange results. The
posterior sample is constant for the fixed effects, and the estimates of the
variance components are way out in the tails of their posterior samples.
The model I'm using says (for l = 1, ..., L - 1)
logit P(X[ijk] = l | X[ijk] = l, U[i], V[j], W[k]) = U[i] + V[j] + W[k] + a[l]
where X[ijk] is the ordinal response to question k for individual j in area i.
The U, V, and W are random effects and the a's are fixed effects. Here's a
function to simulate data which mimics this setup (with a sequence of binary
responses Y[ijkl] = 1 iff X[ijk] = l):
sim - function(n = c(10, 10, 10), sd = c(0.5, 2, 0.5), a = seq(-5, 1, 2)) {
u - rnorm(n[1], 0, sd[1])
v - rnorm(prod(n[1:2]), 0, sd[2])
w - rnorm(n[3], 0, sd[3])
i - factor(rep(1:n[1], each = prod(n[2:3])))
j - factor(rep(1:prod(n[1:2]), each = n[3]))
k - factor(rep(1:n[3], prod(n[1:2])))
df - NULL
for(l in 1:length(a)) {
y - rbinom(length(i), 1, plogis(u[i] + v[j] + w[k] + a[l]))
df - rbind(df, data.frame(i, j, k, l, y))
i - i[!y]
j - j[!y]
k - k[!y]
}
df$l - factor(df$l)
return(df)
}
And here's a function which shows the difficulties I've been having:
test - function(seed = 10, ...) {
require(lme4)
set.seed(seed)
df - sim(...)
df.lmer - lmer(y ~ 0 + l + (1 | i) + (1 | j) + (1 | k), family = binomial,
data = df)
df.mcmc - mcmcsamp(df.lmer, 1000, trans = FALSE)
print(summary(df.lmer))
print(summary(df.mcmc))
densityplot(~ as.numeric(df.mcmc) | factor(rep(colnames(df.mcmc), each =
1000)), scales = list(relation = free))
}
Running
test()
gives the following:
Generalized linear mixed model fit using PQL
Formula: y ~ 0 + l + (1 | i) + (1 | j) + (1 | k)
Data: df
Family: binomial(logit link)
AIC BIClogLik deviance
2316.133 2359.166 -1151.066 2302.133
Random effects:
Groups NameVariance Std.Dev.
j (Intercept) 4.15914 2.03940
k (Intercept) 0.25587 0.50584
i (Intercept) 0.56962 0.75473
number of obs: 3455, groups: j, 100; k, 10; i, 10
Estimated scale (compare to 1) 0.8747598
Fixed effects:
Estimate Std. Error z value Pr(|z|)
l1 -4.502340.40697 -11.0632 2.2e-16 ***
l2 -3.276430.38177 -8.5821 2.2e-16 ***
l3 -1.052770.36566 -2.8791 0.003988 **
l4 0.765380.36832 2.0780 0.037706 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
l1l2l3
l2 0.843
l3 0.841 0.900
l4 0.805 0.865 0.921
l1 l2 l3 l4
Min. :-4.502 Min. :-3.276 Min. :-1.053 Min. :0.7654
1st Qu.:-4.502 1st Qu.:-3.276 1st Qu.:-1.053 1st Qu.:0.7654
Median :-4.502 Median :-3.276 Median :-1.053 Median :0.7654
Mean :-4.502 Mean :-3.276 Mean :-1.053 Mean :0.7654
3rd Qu.:-4.502 3rd Qu.:-3.276 3rd Qu.:-1.053 3rd Qu.:0.7654
Max. :-4.502 Max. :-3.276 Max. :-1.053 Max. :0.7654
j.(In) k.(In) i.(In) deviance
Min. :1.911 Min. :0.0509 Min. :0.06223 Min. :2011
1st Qu.:2.549 1st Qu.:0.1310 1st Qu.:0.19550 1st Qu.:2044
Median :2.789 Median :0.1756 Median :0.25581 Median :2054
Mean :2.824 Mean :0.2085 Mean :0.29948 Mean :2054
3rd Qu.:3.070 3rd Qu.:0.2463 3rd Qu.:0.34640 3rd Qu.:2064
Max. :4.615 Max. :0.8804 Max. :3.62168 Max. :2107
As you can see, the posterior samples from the fixed effects are constant (at
the inital estimates) and the estimates of the variance components aren't within
the IQ ranges of their posterior samples.
I understand from various posts that mcmcsamp is still a work in progress, and
may not work on every model. Is this one of those cases? I'm using R 2.3.1 and
lme4 0.995-2 on Windows XP.
Daniel Farewell
Cardiff University
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] How to convert list elements to data.frames or vectors?
Dear R mailing-list comunity! I'm currently trying to implement an R method. I have two sets of data that I convert into a data.frame each. These data.frames I'd like to append to a list: # generate a list listTable-list() # add one set of data x-1000 ;y-1 ;listTable[[length(listTable) + 1]] - data.frame(matrix(rnorm(x*y), nrow=y)); rm(x); rm(y) # add another set of data (same command) x-1000 ;y-1 ;listTable[[length(listTable) + 1]] - data.frame(matrix(rnorm(x*y), nrow=y)); rm(x); rm(y) My objective is to performed some hypothesis tests on the data. To test if that works out correctly, I tried first using an unpaired t-test (therfore the data.frames consist only of one row each in the example). # alternative, var.equal and conf.level shall # be arguments of my method as well (alternative=two.sided, # var.equal=TRUE, conf.level=0.95) t.test(listTable[[1]][1,], listTable[[2]][1,], alternative=alternative, paired=FALSE, var.equal=var.equal, conf.level=conf.level) And an F-test, throwing an error, like there were not enough observations in the x vector of the test's input. # The F-test (ratio=1, alternative=two.sided, conf.level=0.95) var.test(listTable[[1]][1,], listTable[[2]][1,], ratio=ratio, alternative=alternative, conf.level=conf.level) I figured out, those tests work perfectly, using vectors instead of my list elements with the same argument values, hence there should be no problem with the parameters, I guess. So, my problems would be the list elements like listTable[[1]][1,]. They are no vectors but lists themselves containing each only one element?! I tried several things without any success to change that. I need to have a list like structure and couldn't find a way how to convert the list elements back to data.frames or vectors. Thus I now have a bunch of basic questions on R: 1. If I put a data.frame into a list, how can I get it back as data.frame? 2. How can I get a single row of a data.frame, stored in a list, as vector and not as list of elements? 3. Is a list at all the correct structure for my deeds? 4. Why is this only a problem for the F-test and it seems to be no problem for the t-test? Regards and TIA, Lothar Rubusch __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Spline Extrapolation(NURBS)
Hi R-help
Does R supports Non Uniform Rational B Splines Extrapolation. If yes then
please help me in knowing the way. Further, Is there any reference
regarding Spline Extrapolation (pls note I am not refering to Spline
Interpolation).
thanks a lot in advance
-gaurav
DISCLAIMER AND CONFIDENTIALITY CAUTION:\ \ This message and ...{{dropped}}
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to convert list elements to data.frames or vectors?
-BEGIN PGP SIGNED MESSAGE- Hash: SHA1 Patrick Burns wrote: Chapter 1 of S Poetry might help you come to grips with data structures in R. Patrick Burns [EMAIL PROTECTED] +44 (0)20 8525 0696 http://www.burns-stat.com (home of S Poetry and A Guide for the Unwilling S User) Lothar Botelho-Machado wrote: Dear R mailing-list comunity! I'm currently trying to implement an R method. I have two sets of data that I convert into a data.frame each. These data.frames I'd like to append to a list: # generate a list listTable-list() # add one set of data x-1000 ;y-1 ;listTable[[length(listTable) + 1]] - data.frame(matrix(rnorm(x*y), nrow=y)); rm(x); rm(y) # add another set of data (same command) x-1000 ;y-1 ;listTable[[length(listTable) + 1]] - data.frame(matrix(rnorm(x*y), nrow=y)); rm(x); rm(y) My objective is to performed some hypothesis tests on the data. To test if that works out correctly, I tried first using an unpaired t-test (therfore the data.frames consist only of one row each in the example). # alternative, var.equal and conf.level shall # be arguments of my method as well (alternative=two.sided, # var.equal=TRUE, conf.level=0.95) t.test(listTable[[1]][1,], listTable[[2]][1,], alternative=alternative, paired=FALSE, var.equal=var.equal, conf.level=conf.level) And an F-test, throwing an error, like there were not enough observations in the x vector of the test's input. # The F-test (ratio=1, alternative=two.sided, conf.level=0.95) var.test(listTable[[1]][1,], listTable[[2]][1,], ratio=ratio, alternative=alternative, conf.level=conf.level) I figured out, those tests work perfectly, using vectors instead of my list elements with the same argument values, hence there should be no problem with the parameters, I guess. So, my problems would be the list elements like listTable[[1]][1,]. They are no vectors but lists themselves containing each only one element?! I tried several things without any success to change that. I need to have a list like structure and couldn't find a way how to convert the list elements back to data.frames or vectors. Thus I now have a bunch of basic questions on R: 1. If I put a data.frame into a list, how can I get it back as data.frame? 2. How can I get a single row of a data.frame, stored in a list, as vector and not as list of elements? 3. Is a list at all the correct structure for my deeds? 4. Why is this only a problem for the F-test and it seems to be no problem for the t-test? Regards and TIA, Lothar Rubusch __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Thank you! I was rather concentrated in finding something especially for R and didn't search for S as well. I'll have a look in S Poetry now, It seems like that's exactly the thing I was looking for, the last days! Great! Lothar Rubusch -BEGIN PGP SIGNATURE- Version: GnuPG v1.4.3 (GNU/Linux) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org iD8DBQFE2HrUHRf7N9c+X7sRAvEeAJwPQTjDm0qmtz15mWJIPEmF5atIzQCfa0kC YyzgDhUMoeXdeJ4LzRHWTc4= =MIq/ -END PGP SIGNATURE- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Incomplete gamma function (was integrate(), {was mathematica - r ...})
MM == Martin Maechler [EMAIL PROTECTED]
on Tue, 8 Aug 2006 09:55:50 +0200 writes:
Leo == Leo Gürtler [EMAIL PROTECTED]
on Tue, 08 Aug 2006 00:13:19 +0200 writes:
Leo Dear R-list,
Leo I try to transform a mathematica script to R.
Leo ###relevant part of the Mathematica script
Leo (* p_sv *)
Leo dd = NN (DsD - DD^2);
Leo lownum = NN (L-DD)^2;
Leo upnum = NN (H-DD)^2;
Leo low = lownum/(2s^2);
Leo up = upnum/(2s^2);
Leo psv = NIntegrate[1/(s^NN) Exp[-dd/(2s^2)]
Leo(Gamma[1/2,0,up] + Gamma[1/2,0,low]),{s,sL,sH}, MinRecursion- 3];
LeoPSV = psv/Sqrt[2NN];
Leo Print[- Results ];
Leo Print[ ];
Leo Print[p(sv|D_1D_2I) = const. ,N[PSV,6]];
Leo
Leo # R part
Leo library(fOptions)
Leo ###raw values for reproduction
Leo NN - 58
Leo dd - 0.411769
Leo lownum - 20.81512
Leo upnum - 6.741643
Leo sL - 0.029
Leo sH - 0.092
Leo ###
Leo integpsv - function(s) { 1 / (s^NN) * exp(-dd / (2 * s^2)) *
Leo( (igamma((upnum/(2*s^2)),1/2) - igamma(0,1/2) ) +
Leo(igamma((lownum/(2*s^2)),1/2) - igamma(0,1/2) ) )
Leo }
Leo psv - integrate(integpsv, lower=sL, upper=sH)
Leo PSV - psv$value / sqrt(2*NN)
Leo print(- Results \n)
Leo print(paste(p(sv|D_1D_2I) = const. ,PSV, sep=))
Leo The results of variable PSV are not the same.
Leo In mathematica - PSV ~ 2.67223e+47
Leo with rounding errors due to the initial values, in R - PSV ~ 1.5e+47
Leo I am not that familiar with gamma functions and integration, thus I
Leo assume there the source of the problem can be located.
MM Yes.
MM A few remarks
MM 1) No need to use package fOptions and igamma();
MM standard R's pgamma() is all you need
MM {igamma() has added value only for *complex* arguments!}
MM 2) igamma(0, 1/2) == pgamma(0, 1/2) == 0 , so you can really
MM drop them from your integrand.
MM integpsv - function(s) {
MM1 / (s^NN) * exp(-dd / (2 * s^2)) *
MM( pgamma(upnum/(2*s^2), 1/2) + pgamma(lownum/(2*s^2), 1/2) )
MM }
[]
MM However, if I experiment, using integrate() in two parts, or using many
other
MM numerical integration approximators,
MM all methods give ( your 'psv', not PSV )
MM integrate(integpsv, lower=sL, upper=sH)
MM a value of 1.623779e+48 (which leads to your PSV of 1.5076e+47)
MM Could it be that you are not using the same definition of
MM incomplete gamma in Mathematica and R ?
Offlist, Leo sent me Mathematica's definition
of
Gamma[a, z0, z1] := integral_z0^z1 t^(a-1) exp(-t) dt
Now if you compare this with what ?pgamma (not ?gamma !) tells you,
namely that R uses (Abramowitz and Stegun's definition of the
incomplete gamma function)
pgamma(x, a) = 1/ Gamma(a) * integral_0^x t^(a-1) exp(-t) dt
which has a normalizing factor: In your case above, it is
Gamma(1/2) with its well-known value of sqrt(pi).
And indeed, if you multiply the current result by sqrt(pi), you
get what you want -- and did get from Mathematica:
(1.623779e+48 / sqrt(2*NN)) * sqrt(pi)
[1] 2.672224e+47
Regards,
Martin Maechler, ETH Zurich
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Netiquette, was Re: ... gfortran and gcc...
Thank you both. I would prefer to communicate through the list only. Mike. On Tue August 8 2006 04:47, Prof Brian Ripley wrote: On Tue, 8 Aug 2006, Peter Dalgaard wrote: Prof Brian Ripley [EMAIL PROTECTED] writes: First, you replied to the list and not to me, which was discourteous. You mean that he replied to the list *only*, I hope. Yes, and it was written as if to me, and was a reply to an email from me. I usually consider it offensive when people reply to me and not the list (reasons including: It feels like being grabbed by the sleeve, I might not actually be the best source for the answer, and it's withholding the answer from the rest of the subscribers.) We do ask people to copy to the list. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Frequency Distribution
Hi, Could someone please suggest where I might find some instructions / tutorials / FAQs that describe how to create a frequency distribution and cumulative frequency distribution in R using different class withs. I have about a 2-million observations (distances between points ranging from sub-millimetre to about 400km, and I want to get a feel for how they are distributed). I'd like the output as a table / data rather than an graph. I've searched Google and R's help for obvious terms, and while I've found much information on graphing/plotting, I haven't hit on anything for this. (I only downloaded R about 2 hours ago, apologies if this is obviously documented somewhere I missed.) Regards Michael. Send instant messages to your online friends http://au.messenger.yahoo.com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to export data to Excel Spreadsheet?
Hi On 7 Aug 2006 at 8:00, Berton Gunter wrote: From: Berton Gunter [EMAIL PROTECTED] To: 'Paul Smith' [EMAIL PROTECTED], 'R-Help' R-help@stat.math.ethz.ch Date sent: Mon, 7 Aug 2006 08:00:00 -0700 Organization: Genentech Inc. Subject:Re: [R] How to export data to Excel Spreadsheet? You can also usually copy and paste to/from the Windows clipboard by using file='clipboard' in file i/o or via description = 'clipboard' using connections. I haven't checked all details of this, so there may be some glitches. No problem write.excel-function(tab, ...) write.table( tab, clipboard, sep=\t, row.names=F) works, at least with my version of Excel. Of course after Ctrl-Ving in Excel HTH Petr -- Bert Gunter -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Paul Smith Sent: Monday, August 07, 2006 5:23 AM To: R-Help Subject: Re: [R] How to export data to Excel Spreadsheet? On 8/7/06, Xin [EMAIL PROTECTED] wrote: I try to export my output's data to Excel spreadsheet. My outputs are: comb3 [,1] [,2] [,3] [1,] a b c [2,] a b d [3,] a b e [4,] a b f [5,] a b g See ? write.table ? write.csv Paul __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Petr Pikal [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] (Fwd) Re: paired t-test. Need to rearrange data?
--- Forwarded message follows --- From: Petr Pikal [EMAIL PROTECTED] To: Henrik Parn [EMAIL PROTECTED] Subject:Re: [R] paired t-test. Need to rearrange data? Date sent: Tue, 08 Aug 2006 16:13:47 +0200 Hi Uff, it takes me a bit headache but this shall do it ?unstack ?table ?t new.list-unstack(test.data,y~id) new.test.data-t(as.data.frame(new.list[table(test.data$id)1])) # if you have more than 2 catches you need to modify this^^^ t.test(new.test.data[,1], new.test.data[,2], paired=T) HTH Petr On 7 Aug 2006 at 16:43, Henrik Parn wrote: Date sent: Mon, 07 Aug 2006 16:43:45 +0200 From: Henrik Parn [EMAIL PROTECTED] Send reply to: [EMAIL PROTECTED] Organization: NTNU To: Petr Pikal [EMAIL PROTECTED] Subject:Re: [R] paired t-test. Need to rearrange data? Dear Peter, Thanks a lot for your rapid answer! I am afraid that I presented an example data set in my previous question that was to nicely organised. In fact it more looked like the data I want to acheive rather than the data I recieved... So,may I bother you with a follow-up question: I have data on morphology on birds from several different years. Each year about 50-100 birds are captured, marked with a unique id and measured. The recapture rate is very low and on average 1-2 birds are recaptured the subsequent year. Again a small example data set, which hopefully is more realistic: year - as.factor(rep(1:4,each= 5)) # the four study years id - c(a, b, c, d, e, b, f, g, h, a, i, g, j, k, l, j, m, n, l, o) # id of the birds captured y - rnorm(20) # the measure taken test.data - data.frame(year, id, y) So, in this small example, bird a and b is captured on year 1 and recaptured year 2, bird g captured on year 2 and 3, and so on. On birds that are captured more than once, I would like to do a paired t.test where I compare the measure taken at first capture with that taken on the second time. I suppose I need to add a vector indicating whether it is the first or second time a bird is captured: test.data$time - ifelse(duplicated(test.data$id), 2, 1) But then I am stuck...Do you have a hint of how to proceed? How to select the y-values for which I have both a capture and a recapture? Thanks in advance for taking your time! Best regards, Henrik Petr Pikal wrote: Hi one possibility is t(unstack(test.data, y~id)) HTH Petr On 6 Aug 2006 at 16:13, Henrik Parn wrote: Date sent: Sun, 06 Aug 2006 16:13:29 +0200 From:Henrik Parn [EMAIL PROTECTED] Organization:NTNU To: R-help r-help@stat.math.ethz.ch Subject: [R] paired t-test. Need to rearrange data? Send reply to: [EMAIL PROTECTED] mailto:[EMAIL PROTECTED] mailto:[EMAIL PROTECTED] Dear all, I have received some data arranged like this: # original data id - rep(letters[1:6], each=2) time - as.factor(rep(1:2, 6)) y - as.vector(replicate(6, c(rnorm(n=1, mean=1), rnorm(n=1, mean=2 test.data - data.frame(id, time, y) test.data I would like to perform a paired t-test of the y-values at time=1 against those at time=2, with samples paired by their id. Is it necessary to arrange the data in a format like this: # rearranged data id - letters[1:6] y1 - replicate(6, rnorm(n=1, mean=1)) # y-value at time = 1 y2 - replicate(6, rnorm(n=1, mean=2)) # y-value at time = 2 test.data2 - data.frame(id, y1, y2) test.data2 ...in order to perform a paired t-test? t.test(y1, y2, paired=T) If yes, which is the most convenient way to rearrange the data? Or is it possible to apply the paired t-test function to the original data set? And a side question: In my examples, I suppose can I use set.seed to reproduce the 'rnorm-values' created in the 'original data' also in my the 'rearranged data'. Can someone give me a hint of how to apply the same 'seed' to all the rnorms? Thanks a lot in advance! Henrik __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Petr Pikal [EMAIL PROTECTED] -- Henrik Pärn Department of Biology NTNU 7491 Trondheim Norway +47 735 96282 (office) +47 909 89 255 (mobile) +47 735 96100 (fax) --- End of forwarded message ---Petr Pikal [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal,
[R] Call for Beta Testers: R+ (read R plus) for Solaris and Linux:
Code-named R+Team, our commercially supported R group is looking for beta testers for R+ on Solaris and R+ on Linux. We'd love to get your feedback for our R+Professional version that has additional funcationalities for enterprise support. The sooner we get your feedback and inputs, the faster we'll make changes for the final release! To participate in our beta test program, please email [EMAIL PROTECTED] or go online www.xlsolutions-corp.com/contact.htm The windows (Vista) version of R+ is also into development with a great graphical user interface, and we'll be calling for beta testers this fall or after Vista release. If you want to beta test the vista version, please email Jennifer McDonald, [EMAIL PROTECTED] We look forward to hearing your comments and inputs on R+ ... please feel free to suggest a final name for our commercially supported R. Drew Smith R+ Group Manager 206 686 1578 [EMAIL PROTECTED] www.xlsolutions-corp.com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] More Plots
Hi, How can we plot two graphs ex. lets say correlation ratio in the same window? I mean in the window I have : 1. Graph of correlation having X Y axes 2. Graph of ratio having A B axes one above the other. Thanks, Sonal __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] fixed effects following lmer and mcmcsamp - which to present?
On 8/8/06, Henrik Parn [EMAIL PROTECTED] wrote: Dear all, I am running a mixed model using lmer. In order to obtain CI of individual coefficients I use mcmcsamp. However, I need advice which values that are most appropriate to present in result section of a paper. I have not used mixed models and lmer so much before so my question is probably very naive. However, to avoid to much problems with journal editors and referees addicted to p-values, I would appreciate advice of which values of the output for the fixed factor that would be most appropriate to present in a result section, in order to convince them of the p-value free 'lmer-mcmcsamp'-approach! Using the example from the help page on lmer: fm1 - lmer(Reaction ~ Days + (Days|Subject), sleepstudy) ...I obtain the following for 'Days': summary(fm1) ... Estimate Std. Error t value Days 10.4673 1.54586.771 ...and from mcmcsamp: summary(mcmcsamp(fm1 , n = 1)) 1. Empirical mean and standard deviation for each variable, plus standard error of the mean: Mean SD Naive SE Time-series SE Days 10.4695 1.7354 0.017354 0.015921 2. Quantiles for each variable: 2.5% 25% 50% 75%97.5% Days 7.02279.3395 10.4712 11.5719 13.957 The standard way of presenting coefficients following a 'non-lmer' output is often (beta=..., SE=..., statistic=..., P=...). What would be the best equivalent in a 'lmer-mcmcsamp-context'? (beta=..., CI=...) is a good start I believe. But which beta? And what else? I assume that the a 95% CI in this case would be 7.0227-13.957 (please, do correct me I have completely misunderstood!). But which would be the corresponding beta? 10.4673?, 10.4695? 10.4712? Is the t-value worth presenting or is it 'useless' without corresponding degrees of freedom and P-value? If I present the mcmcsamp-CI, does it make sense to present any of the three SE obtained in the output above? BTW, I have no idea what Naive SE, Time-series SE means. Could not find much in help and pdfs to coda or Matrix, or in Google. I would definitely report the REML value 10.4673 as the estimate. This is a reproducible estimate - the other two are stochastic in that you will get different values from different mcmcsamp runs. However, the reproducibility is high. Notice that they all round to 10.47 and two decimal places is quite enough precision for reporting this value when you consider that a 95% confidence interval is [7.02,13.96]. Another way of evaluating a confidence interval using the coda package is the HPDinterval function. (HPD stands for Highest Posterior Density) It returns the shortest interval with a 95% probability content in the empirical distribution. Here are values from a sample that I evaluated data(sleepstudy) fm1 - lmer(Reaction ~ Days + (Days|Subject), sleepstudy) fs1 - mcmcsamp(fm1, 1) library(coda) summary(fs1) ... 1. Empirical mean and standard deviation for each variable, plus standard error of the mean: Mean SD Naive SE Time-series SE Days 10.4810 1.6910 0.016910 0.015762 2. Quantiles for each variable: 2.5% 25% 50% 75%97.5% Days 7.09759.3971 10.4851 11.5751 13.810 fixef(fm1) (Intercept)Days 251.4051010.46729 HPDinterval(fs1) lower upper Days 7.038076 13.734493 attr(,Probability) [1] 0.95 If you would like to use a t-distribution to calculate a confidence interval, I would argue that the degrees of freedom for the t should be somewhere between n - p = 178 in this case (n = number of observations, p = number of coefficients in the fixed effects or rank(X) where X is the model matrix for the fixed effects) and n - rank([Z:X]) = 144. (there are 36 random effects and 2 fixed effects but the rank of [Z:X] = 36) If we use the lower value of 144 we produce a confidence interval of 10.4673 + c(-1,1) * qt(0.975, df = 144) * 1.5458 [1] 7.41191 13.52269 Notice that this interval is contained in the HPD interval and in the interval obtained from the 2.5% and 97.5% quantiles of the empirical distribution. I attribute this to the fact that the standard errors are calculated conditional on the variance of the random effects. Thus the t-based confidence interval takes into account imprecision of the estimate of \sigma^2 but it assumes that the variance of the random effects is fixed at the estimated values. (That's not quite true - it is the relative variance that is assumed fixed - but the effect is the same.) The MCMC sample allows all the parameters to vary and I would claim is therefore a better measure of the marginal variability of this parameter. However, as stated above, the HPD interval is stochastic. I would create more than one sample
[R] prefixing list names in print
With print(list(A=a,B=b)) it displays $A [1] a $B [1] b I would like to add a common prefix to all the list tags after the $. Pasting the prefix to the names does not work (appear after the $). For example if the prefix would be P, it should display: P$A [1] a P$B [1] b I tried to add a name attribute to the list or to add a prefix=P to print but nothing works. Any hint? Thanks, Laurent. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Netiquette, was Re: ... gfortran and gcc...
What could be the reason, to respond not only to the list? I did not see an advantage, to receive a response twice, once directly, once by the list. Is it wrong, to assume that someone who writes to the list, does also receive all the postings on the list? Heinz At 08:09 08.08.2006 -0500, Mike wrote: Thank you both. I would prefer to communicate through the list only. Mike. On Tue August 8 2006 04:47, Prof Brian Ripley wrote: On Tue, 8 Aug 2006, Peter Dalgaard wrote: Prof Brian Ripley [EMAIL PROTECTED] writes: First, you replied to the list and not to me, which was discourteous. You mean that he replied to the list *only*, I hope. Yes, and it was written as if to me, and was a reply to an email from me. I usually consider it offensive when people reply to me and not the list (reasons including: It feels like being grabbed by the sleeve, I might not actually be the best source for the answer, and it's withholding the answer from the rest of the subscribers.) We do ask people to copy to the list. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Take random sample from class variable
There may be better ways, but this should work: R p.yes - 0.7 R n.yes - rbinom(1, nof.sample, p.yes) R n.no - nof.sample - n.yes R dat.yes - mydat[sample(which(mydat$Class == yes), n.yes, replace=TRUE),] R dat.no - mydat[sample(which(mydat$Class == no), n.no, replace=TRUE),] You can rbind() them, and shuffle the rows if you wish. Andy From: Muhammad Subianto Dear all, Suppose I have a dataset like below, then I take for example, 100 random sample class variable where contains yes and no respectively, 70% and 30%. I need a new 100 random sample from mydat dataset, but I can't get the result. Thanks you very much for any helps. Best, Muhammad Subianto mydat - data.frame(size=c(30,12,15,10,12,12,25,30,20,14), A=c(0,1,0,1,0,1,1,1,0,0), B=c(1,1,0,1,0,1,1,0,0,1), C=c(0,0,1,1,0,0,1,1,0,0), D=c(1,1,1,1,0,1,0,0,1,1), E=c(1,1,0,1,1,1,1,1,1,0), Class=c(yes,yes,no,yes,yes,no,yes,no,yes,yes)) mydat # Maximal data from dataset max.size - sum(mydat$size);max.size # I need sample random nof.sample - 100 set.seed(123) sample.class - sample(c(yes,no), nof.sample, prob=c(.7, .3), replace=TRUE) sample.class sampledat.class - mydat[sample.class,] sampledat.class __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] locating intervals
Hi , I have two sorted vectors X and Xi, where the range of Xi lies within the range of X. For an element in Xi, I want to find the neigbouring data in X, e.g. find an index ix so that for element number k, then X[ix[k]] X[k] X[ix[k] +1] # also OK with = on either one, but not both This is easy to code by looping over the data in X,Xi, but I suspect there may be a faster and more elegant way to do this in R. In Python (Numeric) the same can be achieved with ix=Numeric.searchsorted(X[1:-1],Xi), which is quite compact. So, does anyone know of a corresponding R call that can achive the same? Sincerely, Halldór -- Halldor Bjornsson Weatherservice R D Icelandic Met. Office [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] prefixing list names in print
Laurent Deniau wrote: With print(list(A=a,B=b)) it displays $A [1] a $B [1] b I would like to add a common prefix to all the list tags after the $. Pasting the prefix to the names does not work (appear after the $). For example if the prefix would be P, it should display: P$A [1] a P$B [1] b I tried to add a name attribute to the list or to add a prefix=P to print but nothing works. Any hint? This is a very internal feature of print(). At a first quick look, I think you will have to change the R sources and recompile. Conclusion: Don't do it. Uwe Ligges Thanks, Laurent. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] More Plots
[EMAIL PROTECTED] a écrit : Hi, How can we plot two graphs ex. lets say correlation ratio in the same window? I mean in the window I have : 1. Graph of correlation having X Y axes 2. Graph of ratio having A B axes one above the other. Thanks, Sonal __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Hi, you can simply do that using the command par(new=TRUE) between your two graphs. example : x-rnorm(100,0,2) y-2*x+rnorm(100,0,2) plot(y~x) par(new=TRUE) z-rnorm(100,10,2) plot(z,col='red') As you can notice the second plot is roughly added to the new one, even if the scales are different. So you basically have to specify the axes of the second plot, using the axis() command (after specifying xaxt='n' and yaxt='n' in the second plot for removing its axes). Tristan __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] locating intervals (corrected version)
I have corrected a typo in my previous posting. In what follows the line with the inequality is correct Hi , I have two sorted vectors X and Xi, where the range of Xi lies within the range of X. For an element in Xi, I want to find the neigbouring data in X, e.g. find an index ix so that for element number k, then X[ix[k]] Xi[k] X[ix[k] +1] # also OK with = on either one, but not both This is easy to code by looping over the data in X,Xi, but I suspect there may be a faster and more elegant way to do this in R. In Python (Numeric) the same can be achieved with ix=Numeric.searchsorted(X[1:-1],Xi), which is quite compact. So, does anyone know of a corresponding R call that can achive the same? Sincerely, Halldór -- Halldór Björnsson Deildarstj. Ranns. Þróun Veðursvið Veðurstofu Íslands Halldór Bjornsson Weatherservice R D Icelandic Met. Office [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] prefixing list names in print
On Tue, 8 Aug 2006, Laurent Deniau wrote: With print(list(A=a,B=b)) it displays $A [1] a $B [1] b I would like to add a common prefix to all the list tags after the $. `prefix' ... `after'? You seem to want to prefix component names: why? What do you want for component $a$b$c? For unnamed components? Pasting the prefix to the names does not work (appear after the $). For example if the prefix would be P, it should display: P$A [1] a P$B [1] b I tried to add a name attribute to the list or to add a prefix=P to print but nothing works. Any hint? You will need to alter the C code to do this. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Netiquette, was Re: ... gfortran and gcc...
[Re-sending to the list only for archiving, as my original reply had too many recipients and I cancelled it.] 1. One need not be subscribed to the list to be able to post. Thus, indeed, a poster may not see all postings. 2. On the relatively rare occasion (thanks to Martin) where the server seems to incur delays in sending out posts and replies, copying the original poster on your reply ensures that they will get the reply in a timely fashion. HTH, Marc Schwartz On Tue, 2006-08-08 at 17:41 +0100, Heinz Tuechler wrote: What could be the reason, to respond not only to the list? I did not see an advantage, to receive a response twice, once directly, once by the list. Is it wrong, to assume that someone who writes to the list, does also receive all the postings on the list? Heinz At 08:09 08.08.2006 -0500, Mike wrote: Thank you both. I would prefer to communicate through the list only. Mike. On Tue August 8 2006 04:47, Prof Brian Ripley wrote: On Tue, 8 Aug 2006, Peter Dalgaard wrote: Prof Brian Ripley [EMAIL PROTECTED] writes: First, you replied to the list and not to me, which was discourteous. You mean that he replied to the list *only*, I hope. Yes, and it was written as if to me, and was a reply to an email from me. I usually consider it offensive when people reply to me and not the list (reasons including: It feels like being grabbed by the sleeve, I might not actually be the best source for the answer, and it's withholding the answer from the rest of the subscribers.) We do ask people to copy to the list. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] multinom details
Hello, the multinom() procedure prints a lot of information during the fitting process. Is there a way to disable this verbose mode? Thanks, Bruno __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] multinom details
Try adding trace=FALSE to the call. ?multinom says ... is passed to nnet(), and you'd find trace documented in ?nnet. Please remember to mention the add-on package(s) you're using. Andy From: Bruno L. Giordano Hello, the multinom() procedure prints a lot of information during the fitting process. Is there a way to disable this verbose mode? Thanks, Bruno __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] fixed effects constant in mcmcsamp
Thank you for the thorough report.
On 8/8/06, Daniel Farewell [EMAIL PROTECTED] wrote:
I'm fitting a GLMM to some questionnaire data. The structure is J individuals,
nested within I areas, all of whom answer the same K (ordinal) questions. The
model I'm using is based on so-called continuation ratios, so that it can be
fitted using the lme4 package.
The lmer function fits the model just fine, but using mcmcsamp to judge the
variability of the parameter estimates produces some strange results. The
posterior sample is constant for the fixed effects, and the estimates of the
variance components are way out in the tails of their posterior samples.
The model I'm using says (for l = 1, ..., L - 1)
logit P(X[ijk] = l | X[ijk] = l, U[i], V[j], W[k]) = U[i] + V[j] + W[k] +
a[l]
where X[ijk] is the ordinal response to question k for individual j in area i.
The U, V, and W are random effects and the a's are fixed effects. Here's a
function to simulate data which mimics this setup (with a sequence of binary
responses Y[ijkl] = 1 iff X[ijk] = l):
sim - function(n = c(10, 10, 10), sd = c(0.5, 2, 0.5), a = seq(-5, 1, 2)) {
u - rnorm(n[1], 0, sd[1])
v - rnorm(prod(n[1:2]), 0, sd[2])
w - rnorm(n[3], 0, sd[3])
i - factor(rep(1:n[1], each = prod(n[2:3])))
j - factor(rep(1:prod(n[1:2]), each = n[3]))
k - factor(rep(1:n[3], prod(n[1:2])))
df - NULL
for(l in 1:length(a)) {
y - rbinom(length(i), 1, plogis(u[i] + v[j] + w[k] + a[l]))
df - rbind(df, data.frame(i, j, k, l, y))
i - i[!y]
j - j[!y]
k - k[!y]
}
df$l - factor(df$l)
return(df)
}
And here's a function which shows the difficulties I've been having:
test - function(seed = 10, ...) {
require(lme4)
set.seed(seed)
df - sim(...)
df.lmer - lmer(y ~ 0 + l + (1 | i) + (1 | j) + (1 | k), family = binomial,
data = df)
df.mcmc - mcmcsamp(df.lmer, 1000, trans = FALSE)
print(summary(df.lmer))
print(summary(df.mcmc))
densityplot(~ as.numeric(df.mcmc) | factor(rep(colnames(df.mcmc), each =
1000)), scales = list(relation = free))
}
Running
test()
gives the following:
Generalized linear mixed model fit using PQL
Formula: y ~ 0 + l + (1 | i) + (1 | j) + (1 | k)
Data: df
Family: binomial(logit link)
AIC BIClogLik deviance
2316.133 2359.166 -1151.066 2302.133
Random effects:
Groups NameVariance Std.Dev.
j (Intercept) 4.15914 2.03940
k (Intercept) 0.25587 0.50584
i (Intercept) 0.56962 0.75473
number of obs: 3455, groups: j, 100; k, 10; i, 10
Estimated scale (compare to 1) 0.8747598
Fixed effects:
Estimate Std. Error z value Pr(|z|)
l1 -4.502340.40697 -11.0632 2.2e-16 ***
l2 -3.276430.38177 -8.5821 2.2e-16 ***
l3 -1.052770.36566 -2.8791 0.003988 **
l4 0.765380.36832 2.0780 0.037706 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
l1l2l3
l2 0.843
l3 0.841 0.900
l4 0.805 0.865 0.921
l1 l2 l3 l4
Min. :-4.502 Min. :-3.276 Min. :-1.053 Min. :0.7654
1st Qu.:-4.502 1st Qu.:-3.276 1st Qu.:-1.053 1st Qu.:0.7654
Median :-4.502 Median :-3.276 Median :-1.053 Median :0.7654
Mean :-4.502 Mean :-3.276 Mean :-1.053 Mean :0.7654
3rd Qu.:-4.502 3rd Qu.:-3.276 3rd Qu.:-1.053 3rd Qu.:0.7654
Max. :-4.502 Max. :-3.276 Max. :-1.053 Max. :0.7654
j.(In) k.(In) i.(In) deviance
Min. :1.911 Min. :0.0509 Min. :0.06223 Min. :2011
1st Qu.:2.549 1st Qu.:0.1310 1st Qu.:0.19550 1st Qu.:2044
Median :2.789 Median :0.1756 Median :0.25581 Median :2054
Mean :2.824 Mean :0.2085 Mean :0.29948 Mean :2054
3rd Qu.:3.070 3rd Qu.:0.2463 3rd Qu.:0.34640 3rd Qu.:2064
Max. :4.615 Max. :0.8804 Max. :3.62168 Max. :2107
As you can see, the posterior samples from the fixed effects are constant (at
the inital estimates) and the estimates of the variance components aren't
within
the IQ ranges of their posterior samples.
I understand from various posts that mcmcsamp is still a work in progress, and
may not work on every model. Is this one of those cases? I'm using R 2.3.1 and
lme4 0.995-2 on Windows XP.
In recent versions of the lme4/Matrix packages setting verbose = TRUE
in a call to mcmcsamp for a generalized linear mixed model causes the
Metropolis-Hastings ratio for the proposed change in the fixed effects
and random effects to be printed. When I ran your example those
values were always 0 which is either extremely bad luck or a bug. My
guess is that it's a bug.
Thanks for the report with the reproducible example.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide
Re: [R] Sampling from a Matrix
From: Marc Schwartz
On Fri, 2006-08-04 at 12:46 -0400, Daniel Gerlanc wrote:
Hello all,
Consider the following problem:
There is a matrix of probabilities:
set.seed(1)
probs - array(abs(rnorm(25, sd = 0.33)), dim = c(5,5),
dimnames =
list(1:5, letters[1:5])) probs
a b c de
1 0.21 0.27 0.50 0.0148 0.303
2 0.06 0.16 0.13 0.0053 0.258
3 0.28 0.24 0.21 0.3115 0.025
4 0.53 0.19 0.73 0.2710 0.656
5 0.11 0.10 0.37 0.1960 0.205
I want to sample 3 values from each row.
One way to do this follows:
index - 1:ncol(probs)
for(i in 1:nrow(probs)){
## gets the indexes of the values chosen
sample(index, size = 3, replace = TRUE, prob = probs[i, ])
}
Is there a another way to do this?
Thanks!
t(apply(probs, 1, function(x) sample(x, 3)))
[,1] [,2] [,3]
1 0.210 0.5000 0.0148
2 0.258 0.0053 0.1300
3 0.025 0.2800 0.3115
4 0.190 0.5300 0.2710
5 0.196 0.1000 0.1100
Hmm... If I read Daniel's code (which is different from his description)
correctly, that doesn't seem to be what he wanted. Perhaps something like
this:
apply(probs, 1, function(p) sample(1:ncol(probs), 3, replace=TRUE, prob=p))
Andy
See ?apply and ?t
HTH,
Marc Schwartz
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Fitting data with optim or nls--different time scales
Hi,
I have a system of ODE's I can solve with lsoda.
Model=function(t,x,parms)
{
#parameter definitions
lambda=parms[1]; beta=parms[2];
d = parms[3]; delta = parms[4];
p=parms[5];c=parms[6]
xdot[1] = lambda - (d*x[1])- (beta*x[3]*x[1])
xdot[2] = (beta*x[3]*x[1]) - (delta*x[2])
xdot[3] = (p*x[2]) - (c*x[3])
return(list(xdot))
}
I want to fit the output out[,4] to experimental data that is only
available on days 0, 7, 12, 14, 17, and 20. I don't know how to set up
optim or nls so that it takes out[,4] on the appropriate day, but still
runs lsoda on a time scale of 0.01 day.
Below is the function I've been using to run 'optim', at the
course-grained time scale:
Modelfit=function(s) {
parms[1:4]=s[1:4]; times=c(0,7,12,14,17,20,25)
out=lsoda(x0,times,Model,parms)
mse=mean((log10(out[,4])-log10(i(times)))^2)
# cat(times)
return(mse)
}
#x0=c(T0,I0,V0)
x0=c(2249,0,1)
#parms(lambda, beta, d, delta, p, c)
parms[5:6]=c(1.0,23)
s0=c(49994,8456,6.16E-8,0.012) #initial values
fit=optim(s0,Modelfit)
Right now, lsoda is being run on too course-grained a time scale in the
function Modelfit. Most examples of optim and nls I have found compare
two data sets at the same times, and run lsoda on the time scale the
data is available at, but I would like to run lsoda at a finer scale, and
only compare the appropriate time points with the experiment. I have also
tried using nls, but I have the same problem. Does anyone have
suggestions?
Thank you very much,
Leslie
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] locating intervals (corrected version)
?findInterval() could be of help in this case. Best, Dimitris Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/(0)16/336899 Fax: +32/(0)16/337015 Web: http://med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm Quoting halldor bjornsson [EMAIL PROTECTED]: I have corrected a typo in my previous posting. In what follows the line with the inequality is correct Hi , I have two sorted vectors X and Xi, where the range of Xi lies within the range of X. For an element in Xi, I want to find the neigbouring data in X, e.g. find an index ix so that for element number k, then X[ix[k]] Xi[k] X[ix[k] +1] # also OK with = on either one, but not both This is easy to code by looping over the data in X,Xi, but I suspect there may be a faster and more elegant way to do this in R. In Python (Numeric) the same can be achieved with ix=Numeric.searchsorted(X[1:-1],Xi), which is quite compact. So, does anyone know of a corresponding R call that can achive the same? Sincerely, Halldór -- Halldór Björnsson Deildarstj. Ranns. Þróun Veðursvið Veðurstofu Íslands Halldór Bjornsson Weatherservice R D Icelandic Met. Office [[alternative HTML version deleted]] Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] More Plots
--- [EMAIL PROTECTED] wrote: Hi, How can we plot two graphs ex. lets say correlation ratio in the same window? I mean in the window I have : 1. Graph of correlation having X Y axes 2. Graph of ratio having A B axes one above the other. Thanks, Sonal ?par and have a look at mfcol or mfrow. Example a - c(1:5) b - c(11:15) x - c(21:30) y - c(31:40) par(mfcol=c(2,1)) plot(a,b) plot (x,y) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] prefixing list names in print
A simple function will do what you want, customize this as needed:
lprint - function(lst,prefix)
{
for (i in 1:length(lst)) {
cat(paste(prefix,$,names(lst)[i],sep=),\n)
print(lst[[i]])
cat(\n)
}
}
P - list(A=a,B=b)
lprint(P,Prefix)
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Laurent Deniau
Sent: Tuesday, August 08, 2006 12:25 PM
To: R-help
Subject: [R] prefixing list names in print
With
print(list(A=a,B=b))
it displays
$A
[1] a
$B
[1] b
I would like to add a common prefix to all the list tags after the $.
Pasting the prefix to the names does not work (appear after the $).
For example if the prefix would be P, it should display:
P$A
[1] a
P$B
[1] b
I tried to add a name attribute to the list or to add a
prefix=P to
print but nothing works. Any hint?
Thanks,
Laurent.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] trellis in black and white
I am using this code but the output is not totally black and white. The dots and lines are green and the panel labels are light green or tan. I need to adapt this code to get just black on a wight background. trellis.device(device=jpeg, file=C:/Documents and Settings/tables/figure2_March27b.jpg,theme=col.whitebg) print( xyplot(AWGT ~ log(pcb_80) | malex*romanix, data=pcb_graph3a, auto.key = list(lines = TRUE, points = TRUE), ylab=Birth Weight, xlab=log PCB, type=c(p, smooth), span=.8) ) dev.off() thanks -- Dean Sonneborn, MS Programmer Analyst Department of Public Health Sciences University of California, Davis (530) 754-9516 [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Frequency Distribution
Thankyou William. I found the package and read through the documentation. I'm not a statistican, so it was largely over my head. I was looking for a command/function that described itself as performing a frequency distribution, and could not find anything obvious enough. What did you have in mind in the package that you thought may help? All I'm looking to do is ask it to give me frequencies and cumulative frequencies for the whole dataset, using intervale widths of 100 or 1000 (in much the same way the data would have to have been binned before producing a histogram. Regards Michael. --- William Asquith [EMAIL PROTECTED] wrote: You might be interested in the lmomco package that supports many nontraditional and traditional distributions. William A. On Aug 8, 2006, at 9:00 AM, Michael Zatorsky wrote: Hi, Could someone please suggest where I might find some instructions / tutorials / FAQs that describe how to create a frequency distribution and cumulative frequency distribution in R using different class withs. I have about a 2-million observations (distances between points ranging from sub-millimetre to about 400km, and I want to get a feel for how they are distributed). I'd like the output as a table / data rather than an graph. I've searched Google and R's help for obvious terms, and while I've found much information on graphing/plotting, I haven't hit on anything for this. (I only downloaded R about 2 hours ago, apologies if this is obviously documented somewhere I missed.) Regards Michael. Send instant messages to your online friends http:// au.messenger.yahoo.com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting- guide.html and provide commented, minimal, self-contained, reproducible code. Send instant messages to your online friends http://au.messenger.yahoo.com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Help with short time series
Thanks for the help. I provide below the dataset I'm using, it's a little bit different from what I was describing (sorry for that). The streams are 3 and I have an unequal number of years for each stream. Stream Density Year 1 Zak0.20 2000 2 Zak0.36 2001 3 Zak0.41 2002 4 Zak0.34 2003 5 Zak0.28 2004 6 Gor0.08 1999 7 Gor0.05 2000 8 Gor0.14 2001 9 Gor0.16 2002 10Gor0.13 2003 11Gat0.18 2004 12Gat0.10 2001 13Gat0.37 2002 14Gat0.57 2003 15Gat0.47 2004 I tried to follow the suggestions of Dieter, but the model does not fit. Any suggestion will be appreciated Dear R-list, I have a statistical problem with the comparison of two short time-series of density data in an ecological framework. I have to compare two short time series (5 years, one value for each year) of species density data (it is the density of fish in two different streams) to test if the two means of the five densities are significantly different, so basically if the two mean stream-specific fish densities are significantly different. I don't think I can use a straight t-test due to the problem of autocorrelation and I don't think I can use a repeated measure ANOVA as I don't have any replicates. Any help would be greatly appreciated. try something like library(nlme) summary(lme(dens~stream+year,data=mystreamdata,random=~year|stream)) This should also give you an estimate if the slopes are different if you test against the simplified model summary(lme(dens~stream+year,data=mystreamdata,random=~1|stream)) Since you did not provide a short example data set, this is only approximatively right. Dieter -- Universita' degli Studi di Parma (http://www.unipr.it) _ Simone Vincenzi, PhD Student Department of Environmental Sciences University of Parma Parco Area delle Scienze, 33/A, 43100 Parma, Italy Phone: +39 0521 905696 Fax: +39 0521 906611 e.mail: [EMAIL PROTECTED] -- -- Universita' degli Studi di Parma (http://www.unipr.it) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] parameter yaxs / function hist (graphics)
Dear Paulo, Thank you for your reply. But I doubt yours is a proper solution to my request, for some reasons: 1. The position of the axis graphed with the command axis(1, line=-1) depends on the size of the graphics window. 2. After your graph is on the screen, in case one may want a boxed graph, a box() command will produce a histogram floating in the air, not touching the horizontal axis. Of course, one could build a proper box (with labels, etc.) by means of more primitive graphics functions like lines (package graphics) and others, but I think that would mean a lot of work. Thank you again. Paulo Barata --- Paulo Barata Fundacao Oswaldo Cruz / Oswaldo Cruz Foundation Rua Leopoldo Bulhoes 1480 - 8A 21041-210 Rio de Janeiro - RJ Brasil E-mail: [EMAIL PROTECTED] --- Paulo Justiniano Ribeiro Jr wrote: Paulo One possibility is to draw the histogram without axes and then add them wherever you want. For instance with something along the lines: x - rnorm(500) hist(x, axes=F) axis(1, line=-1) For more details: ?axis best P.J. Paulo Justiniano Ribeiro Jr LEG (Laboratório de Estatística e Geoinformação) Departamento de Estatística Universidade Federal do Paraná Caixa Postal 19.081 CEP 81.531-990 Curitiba, PR - Brasil Tel: (+55) 41 3361 3573 Fax: (+55) 41 3361 3141 e-mail: [EMAIL PROTECTED] http://www.est.ufpr.br/~paulojus On Mon, 7 Aug 2006, Paulo Barata wrote: Dear R users, The parameters xaxs and yaxs (function par, package graphics) seem not to work with the function hist (package graphics), even when the parameters xlim and ylim are defined. Is there any way to make yaxs=i and xaxs=i work properly with the function hist, mainly to produce histograms that touch the horizontal axis? The R documentation and the R mailing lists archive don't seem to be of help here. I am using R 2.3.1, running under Windows XP. ## Example: x - rnorm(100) hist(x,breaks=seq(-4,4,0.5),ylim=c(0,40),yaxs=i, xlim=c(-4,4),xaxs=i) box() Thank you very much. Paulo Barata -- Paulo Barata Fundacao Oswaldo Cruz / Oswaldo Cruz Foundation Rua Leopoldo Bulhoes 1480 - 8A 21041-210 Rio de Janeiro - RJ Brasil E-mail: [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Getting data out of a loop
A stupid question but I just cannot see how to do
this.
I have a loop that does some calculations and puts the
results in a vector for each iteration, but I cannot
see how to get the data out of the loop in such a way
that I can use it. I can print it but how do I get it
into a set of vectors or what ever.
Any help gratefully received. Thanks
Example
cata - c( 3,5,6,8,0, NA)
catb - c( 1,2,3,4,5,6)
doga - c(3,5,3,6,4, 0)
dogb - c(2,4,6,8,10, 12)
rata - c (NA, 5, 5, 4, 9, 0)
ratb - c( 1,2,3,4,5,6)
bata - c( 12, 42,NA, 45, 32, 54)
batb - c( 13, 15, 17,19,21,23)
id - Cs(a,b,b,c,a,b)
site - c(1,1,4,4,1,4)
mat1 - cbind(cata, catb, doga, dogb, rata, ratb,
bata, batb)
Df - data.frame(site, id, mat1)
nn - levels(Df$id)
Df
nn
rate - c(2,3,4)
for (i in 1: length(nn)) {
dd- subset(Df, id==nn[i])
scat - sum(c(dd$cata,dd$catb), na.rm=T)
sdog - sum(c(dd$doga,dd$dogb), na.rm=T)
srat - sum(c(dd$rata, dd$ratb), na.rm=T)
sbat - sum(c(dd$bata,dd$batb), na.rm=T)
sss - c(scat,sdog, srat,sbat) * rate[i]
print(sss)
}
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Getting data out of a loop
Store it in a 'list'; see modified script below:
cata - c( 3,5,6,8,0, NA)
catb - c( 1,2,3,4,5,6)
doga - c(3,5,3,6,4, 0)
dogb - c(2,4,6,8,10, 12)
rata - c (NA, 5, 5, 4, 9, 0)
ratb - c( 1,2,3,4,5,6)
bata - c( 12, 42,NA, 45, 32, 54)
batb - c( 13, 15, 17,19,21,23)
id - c('a', 'b', 'b', 'c', 'a', 'b')
site - c(1,1,4,4,1,4)
mat1 - cbind(cata, catb, doga, dogb, rata, ratb,
bata, batb)
Df - data.frame(site, id, mat1)
nn - levels(Df$id)
Df
nn
rate - c(2,3,4)
Result - list()
for (i in 1: length(nn)) {
dd- subset(Df, id==nn[i])
scat - sum(c(dd$cata,dd$catb), na.rm=T)
sdog - sum(c(dd$doga,dd$dogb), na.rm=T)
srat - sum(c(dd$rata, dd$ratb), na.rm=T)
sbat - sum(c(dd$bata,dd$batb), na.rm=T)
sss - c(scat,sdog, srat,sbat) * rate[i]
Result[[i]] - sss
print(sss)
}
Result
On 8/8/06, John Kane [EMAIL PROTECTED] wrote:
A stupid question but I just cannot see how to do
this.
I have a loop that does some calculations and puts the
results in a vector for each iteration, but I cannot
see how to get the data out of the loop in such a way
that I can use it. I can print it but how do I get it
into a set of vectors or what ever.
Any help gratefully received. Thanks
Example
cata - c( 3,5,6,8,0, NA)
catb - c( 1,2,3,4,5,6)
doga - c(3,5,3,6,4, 0)
dogb - c(2,4,6,8,10, 12)
rata - c (NA, 5, 5, 4, 9, 0)
ratb - c( 1,2,3,4,5,6)
bata - c( 12, 42,NA, 45, 32, 54)
batb - c( 13, 15, 17,19,21,23)
id - Cs(a,b,b,c,a,b)
site - c(1,1,4,4,1,4)
mat1 - cbind(cata, catb, doga, dogb, rata, ratb,
bata, batb)
Df - data.frame(site, id, mat1)
nn - levels(Df$id)
Df
nn
rate - c(2,3,4)
for (i in 1: length(nn)) {
dd- subset(Df, id==nn[i])
scat - sum(c(dd$cata,dd$catb), na.rm=T)
sdog - sum(c(dd$doga,dd$dogb), na.rm=T)
srat - sum(c(dd$rata, dd$ratb), na.rm=T)
sbat - sum(c(dd$bata,dd$batb), na.rm=T)
sss - c(scat,sdog, srat,sbat) * rate[i]
print(sss)
}
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Jim Holtman
Cincinnati, OH
+1 513 646 9390
What is the problem you are trying to solve?
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Split-plot model
How do I set up my model equation in aov to analyze a split-plot design? I have two factors (CO2 and NITROGEN), each with two levels (high and ambient). CO2 is my whole-plot factor with three replicates for each level (i.e., 6 rooms total). Is this syntax below correct? summary(aov(response ~ ROOM + CO2*NITROGEN + Error(ROOM/CO2))) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Netiquette, was Re: ... gfortran and gcc...
I agree. Also, sending a copy to the poster means that they are likely to get it first which seems like a desirable courtesy. On 8/8/06, Marc Schwartz (via MN) [EMAIL PROTECTED] wrote: [Re-sending to the list only for archiving, as my original reply had too many recipients and I cancelled it.] 1. One need not be subscribed to the list to be able to post. Thus, indeed, a poster may not see all postings. 2. On the relatively rare occasion (thanks to Martin) where the server seems to incur delays in sending out posts and replies, copying the original poster on your reply ensures that they will get the reply in a timely fashion. HTH, Marc Schwartz On Tue, 2006-08-08 at 17:41 +0100, Heinz Tuechler wrote: What could be the reason, to respond not only to the list? I did not see an advantage, to receive a response twice, once directly, once by the list. Is it wrong, to assume that someone who writes to the list, does also receive all the postings on the list? Heinz At 08:09 08.08.2006 -0500, Mike wrote: Thank you both. I would prefer to communicate through the list only. Mike. On Tue August 8 2006 04:47, Prof Brian Ripley wrote: On Tue, 8 Aug 2006, Peter Dalgaard wrote: Prof Brian Ripley [EMAIL PROTECTED] writes: First, you replied to the list and not to me, which was discourteous. You mean that he replied to the list *only*, I hope. Yes, and it was written as if to me, and was a reply to an email from me. I usually consider it offensive when people reply to me and not the list (reasons including: It feels like being grabbed by the sleeve, I might not actually be the best source for the answer, and it's withholding the answer from the rest of the subscribers.) We do ask people to copy to the list. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] More Plots
How can we plot two graphs ex. lets say correlation ratio in the same window? I mean in the window I have : 1. Graph of correlation having X Y axes 2. Graph of ratio having A B axes one above the other. Why do you want to do this? It is not a good idea unless you are trying to confusing or mislead people. Hadey __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Getting data out of a loop
I'd create an empty dataframe prior to the loop.
cata - c( 3,5,6,8,0, NA)
catb - c( 1,2,3,4,5,6)
doga - c(3,5,3,6,4, 0)
dogb - c(2,4,6,8,10, 12)
rata - c (NA, 5, 5, 4, 9, 0)
ratb - c( 1,2,3,4,5,6)
bata - c( 12, 42,NA, 45, 32, 54)
batb - c( 13, 15, 17,19,21,23)
id - c('a', 'b', 'b', 'c', 'a', 'b')
site - c(1,1,4,4,1,4)
mat1 - cbind(cata, catb, doga, dogb, rata, ratb,
bata, batb)
Df - data.frame(site, id, mat1)
nn - levels(Df$id)
Df
nn
rate - c(2,3,4)
Result - data.frame(matrix(NA,length(nn),4))
for (i in 1: length(nn)) {
dd- subset(Df, id==nn[i])
scat - sum(c(dd$cata,dd$catb), na.rm=T)
sdog - sum(c(dd$doga,dd$dogb), na.rm=T)
srat - sum(c(dd$rata, dd$ratb), na.rm=T)
sbat - sum(c(dd$bata,dd$batb), na.rm=T)
sss - c(scat,sdog, srat,sbat) * rate[i]
Result[i,] - sss
print(sss)
}
Result
From: John Kane [EMAIL PROTECTED]
To: R R-help r-help@stat.math.ethz.ch
Subject: [R] Getting data out of a loop
Date: Tue, 8 Aug 2006 18:00:23 -0400 (EDT)
A stupid question but I just cannot see how to do
this.
I have a loop that does some calculations and puts the
results in a vector for each iteration, but I cannot
see how to get the data out of the loop in such a way
that I can use it. I can print it but how do I get it
into a set of vectors or what ever.
Any help gratefully received. Thanks
Example
cata - c( 3,5,6,8,0, NA)
catb - c( 1,2,3,4,5,6)
doga - c(3,5,3,6,4, 0)
dogb - c(2,4,6,8,10, 12)
rata - c (NA, 5, 5, 4, 9, 0)
ratb - c( 1,2,3,4,5,6)
bata - c( 12, 42,NA, 45, 32, 54)
batb - c( 13, 15, 17,19,21,23)
id - Cs(a,b,b,c,a,b)
site - c(1,1,4,4,1,4)
mat1 - cbind(cata, catb, doga, dogb, rata, ratb,
bata, batb)
Df - data.frame(site, id, mat1)
nn - levels(Df$id)
Df
nn
rate - c(2,3,4)
for (i in 1: length(nn)) {
dd- subset(Df, id==nn[i])
scat - sum(c(dd$cata,dd$catb), na.rm=T)
sdog - sum(c(dd$doga,dd$dogb), na.rm=T)
srat - sum(c(dd$rata, dd$ratb), na.rm=T)
sbat - sum(c(dd$bata,dd$batb), na.rm=T)
sss - c(scat,sdog, srat,sbat) * rate[i]
print(sss)
}
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] unable to restore saved data
Lately, when I try to open R I get the following error message: Error: object 'time' not found whilst loading namespace 'tseries' Fatal error: unable to restore saved data in .RData If I rename .RData to RData.RData and then try opening R again it works. Then I can load(RData.RData) without a problem. But if I try saving my workspace (as the default, .RData) and reload R it crashes all over again. I don't know how to get this 'time' object back. (I must have removed it by accident at some point.) Any ideas? Thanks in advance, -- Alessandro Gagliardi Integrative Neuroscience Program Rutgers University Mind Brain Analysis [EMAIL PROTECTED] The opposite of a correct statement is a false statement. But the opposite of a profound truth may well be another profound truth. -Niels Bohr __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Frequency Distribution
You could just do table(cut(...)) and cumsum(table(cut(...))). See the help pages for those functions. Example: R x - rnorm(1e4) R breaks - c(-Inf, -3:3, Inf) R table(cut(x, breaks)) (-Inf,-3] (-3,-2] (-2,-1](-1,0] (0,1] (1,2] (2,3] (3, Inf] 16 253 1389 3349 3419 1339 220 15 R cumsum(table(cut(x, breaks))) (-Inf,-3] (-3,-2] (-2,-1](-1,0] (0,1] (1,2] (2,3] (3, Inf] 16 269 1658 5007 8426 9765 9985 1 Andy From: Michael Zatorsky Thankyou William. I found the package and read through the documentation. I'm not a statistican, so it was largely over my head. I was looking for a command/function that described itself as performing a frequency distribution, and could not find anything obvious enough. What did you have in mind in the package that you thought may help? All I'm looking to do is ask it to give me frequencies and cumulative frequencies for the whole dataset, using intervale widths of 100 or 1000 (in much the same way the data would have to have been binned before producing a histogram. Regards Michael. --- William Asquith [EMAIL PROTECTED] wrote: You might be interested in the lmomco package that supports many nontraditional and traditional distributions. William A. On Aug 8, 2006, at 9:00 AM, Michael Zatorsky wrote: Hi, Could someone please suggest where I might find some instructions / tutorials / FAQs that describe how to create a frequency distribution and cumulative frequency distribution in R using different class withs. I have about a 2-million observations (distances between points ranging from sub-millimetre to about 400km, and I want to get a feel for how they are distributed). I'd like the output as a table / data rather than an graph. I've searched Google and R's help for obvious terms, and while I've found much information on graphing/plotting, I haven't hit on anything for this. (I only downloaded R about 2 hours ago, apologies if this is obviously documented somewhere I missed.) Regards Michael. Send instant messages to your online friends http:// au.messenger.yahoo.com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting- guide.html and provide commented, minimal, self-contained, reproducible code. Send instant messages to your online friends http://au.messenger.yahoo.com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] boxplot
When I tried: R x - c(5, rnorm(30)) R boxplot(x) R identify(x) [1] 1 (after I clicked on the obvious outlier, and R labeled it `1') it seems to work just fine. What do you mean by can't apply it boxplot? Andy From: Ana Patricia Martins Hello R-users and developers, Once again, I'm asking for your help. I've used identify to identify points in a scatter plot. However, I can't apple in the boxplot. I need to identify the outlier's id in the boxplot. Can anyone help me? Thanks in advance, Ana Patricia Martins [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Spline Extrapolation(NURBS) 2nd Post
Hi R-help
Does R supports Non Uniform Rational B Splines (NURBS) Extrapolation. If
yes then
please help me in knowing the way. Further, Is there any reference
regarding Spline Extrapolation (pls note I am not refering to Spline
Interpolation). If anybody has any idea as to how to do spline
extrapolation or any references
in regards to the same then please enlighten me either through the list or
at my personal id i.e.
[EMAIL PROTECTED]
thanks a lot in advance
-gaurav
DISCLAIMER AND CONFIDENTIALITY CAUTION:\ \ This message and ...{{dropped}}
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] debug print() commands not showing during file write loop
Hello, I have a for loop that takes about and hour to complete each loop. I added a print() command at the end of the looped code as a way to see its progress, but the output is suppressed until the entire for-loop is finished. why is that? can it be changed such that the print() output is echoed to the screen when it is processed? thanks in advance, -rich __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] How to draw the decision boundaries for LDA and Rpart object
Hello useR, Could you please tell me how to draw the decision boundaries in a scatterplot of the original data for a LDA or Rpart object. For example: library(rpart) fit.rpart - rpart(as.factor(group.id)~., data=data.frame(Data) ) How can I draw the cutting lines on the orignial Data? Or is there any built in functions that can read the rpart object 'fit.rpart' to do that? Thanks very much in advance! Leon [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] NLME: Problem with plotting ranef vs a factor
Your question is entirely too complex for me to try to answer in a reasonable amount of time, especially since your example in not self contained. If you would still like help on this, I suggest you try to generate a self contained example that is as simple as you can make it that illustrates your problem, as suggested in the posting guide www.R-project.org/posting-guide.html. With only a modest amount of luck, the things you try to simplify your example will lead to enlightenment. If that fails, please submit another question that is self contained, simple and clear. Doing so should substantially increase your chances of getting a quick, useful reply. I know this doesn't answer your question, but I hope it helps. Spencer Graves Greg Distiller wrote: Hi I am following the model building strategy that is outlined in the Pinheiro and Bates book wrt including covariates but am having a problem with the plot. Basically I am using 4 covariates (1 of them is continuous) and 3 of them are fine but the 4th one is being shown as a scatterplot despite the fact that it is a factor. I have explicitly declared this to be a factor (pcat-as.factor(pcat)) and have also checked by using the is.factor and the levels command that it is a factor. Yet despite this the plot command is not recognising it as a factor. Here is more information about my problem: I am reading in the data by: Data-read.csv(Data1_93_2.csv,header=T) attach(Data) Data1_93-transform(Data,log2game=log2(gamedens+1)) pcat-as.factor(pcat) Data1_93-groupedData(log2game ~ day | subjectno, data=Data1_93) detach(Data) Here is the code to check that the covariate called pcat is indeed a factor: levels(pcat) [1] 1 2 3 is.factor(pcat) [1] TRUE and then after the model is fitted I extract the random effects: D1C2.ran - ranef(mod11.103nlme,augFrame=T) and here is an extract from the object: C R day gamedens pcat site mutcat1 pdens0 log2game NA02_259 -1.016987007 0.0162825099 15.75000 23.51 Namaacha Mixed 15018 3.761099 NA02_073 -0.939355374 0.0132589702 10.5 23.750001 Namaacha Resistant6170 3.675543 M00_12 -0.775048474 0.0047124742 10.5 25.01 Mpumulanga Sensitive 17525 3.768326 M00_93 -0.555801118 0.0053872868 14.0 37.52 Mpumulanga Sensitive 332000 4.254319 NA02_053 -0.327990343 -0.0037659864 6.0 39.250001 Namaacha Resistant 65529 4.292481 Note that this output also seems to indicate that pcat is a factor as it is summarised correctly. I then generate plots for my random effects: plot(D1C2.ran,form= C ~site+mutcat2+pcat+pdens0) and the problem is that the panel for my random effects vs pcat is displayed as a scatterplot rather than as a boxplot. I am getting told to check warnings and these warnings look like: Warning messages: 1: at 0.99 2: radius 0.0001 3: all data on boundary of neighborhood. make span bigger 4: pseudoinverse used at 0.99 5: neighborhood radius 0.01 6: reciprocal condition number -1.#IND 7: zero-width neighborhood. make span bigger I do not get these warnings if I exclude the problematic variable pcat so must be something to do with this. Any ideas? Many thanks Greg [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
