Re: [R] Sum of Bernoullis with varying probabilities
Or perhaps its clearer (and saves a bit of space) to use apply...prod
here instead of exp...log:
fft(apply(mvfft(t(cbind(1-p,p,0,0,0))), 1, prod), inverse = TRUE)/5
On 10/7/06, Gabor Grothendieck [EMAIL PROTECTED] wrote:
One can get a one-line solution by taking the product of the FFTs.
For example, let p - 1:4/8 be the probabilities. Then the solution is:
fft(exp(rowSums(log(mvfft(t(cbind(1-p,p,0,0,0)), inverse = TRUE)/5
On 10/7/06, Ted Harding [EMAIL PROTECTED] wrote:
Hi again.
I had suspected that doing the calculation by a convolution
method might be both straightforward and efficient in R.
I've now located convolve() (in 'base'!!), and have a solution
using this function. Details below. Original statement of
problem, and method originally proposed, repeated below for
reference.
On 10/6/06, Ted Harding [EMAIL PROTECTED] wrote:
Hi Folks,
Given a series of n independent Bernoulli trials with
outcomes Yi (i=1...n) and Prob[Yi = 1] = Pi, I want
P = Prob[sum(Yi) = r] (r = 0,1,...,n)
I can certainly find a way to do it:
Let p be the vector c(P1,P2,...,Pn).
The cases r=0 and r=n are trivial (and also are exceptions
for the following routine).
For a given value of r in (1:(n-1)),
library(combinat)
Set - (1:n)
Subsets - combn(Set,r)
P - 0
for(i in (1:dim(Subsets)[2])){
ix - numeric(n)
ix[Subsets[,i]] - 1
P - P + prod((p^ix) * ((1-p)^(1-ix)))
}
The convolution method implemented in convolve computes, for
two series X[1],X[2],...,X[n] and Y[1],Y[2],...,Y[m] the
series (k = 1, 2, ... , n+m-1):
Z[k] = sum( X[k-m+i]*Y[i] )
over valid values of i, i.e. products of terms of X matched
with terms of a shifted Y, using the open method (see
?convolve).
This is not quite what is wanted for the type of convolution
needed to compute the distribution of the sum of two integer
random variables, since one needs to reverse one of the series
being convolved. This is the basis of the following:
Given a vector p of probabilities P1, P2, P3, for Y=1 in
successive trials, I need to convolve the successive Bernoulli
distributions (1-P1,P1), (1-P2,P2), ... , (1-Pn,Pn). Hence:
ps-cbind(1-p,p); u-ps[1,]; u1-u
for(i in (2:16)){
u-convolve(u,rev(ps[i,]),type=o)
}
In the case I was looking at, I had already obtained the vector
P of probabilities P(Sum(Y)=0), P(sum(Y)=1), ... by the method
previously described (repeated above); and here n=16.
Having obtained u by the above routine, I can now compare the
results, u with P:
cbind((0:16),u,P)
r uP
0 1.353353e-01 1.353353e-01
1 3.007903e-01 3.007903e-01
2 2.976007e-01 2.976007e-01
3 1.747074e-01 1.747074e-01
4 6.826971e-02 6.826971e-02
5 1.884969e-02 1.884969e-02
6 3.804371e-03 3.804371e-03
7 5.720398e-04 5.720398e-04
8 6.463945e-05 6.463945e-05
9 5.489454e-06 5.489454e-06
10 3.473997e-07 3.473997e-07
11 1.607822e-08 1.607822e-08
12 5.262533e-10 5.262532e-10
13 1.148626e-11 1.148618e-11
14 1.494650e-13 1.493761e-13
15 9.008700e-16 8.764298e-16
16 -1.896716e-17 1.313261e-19
so, apart from ignorable differences in the very small
probabilities for r11, they agree. I have to admit, though,
that I had to tread carefully (and experiment with Binomials)
to make sure of exactly how to introduce the reversal
(at rev(ps[i,])).
And the convolution method is *very* much faster than the
selection of all subsets method!
However, I wonder whether the subsets method may be the more
accurate of the two (not that it really matters).
Best wishes to all,
Ted.
E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 094 0861
Date: 07-Oct-06 Time: 22:03:27
-- XFMail --
__
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] Sum of Bernoullis with varying probabilities
On 08-Oct-06 Gabor Grothendieck wrote:
Or perhaps its clearer (and saves a bit of space) to use apply...prod
here instead of exp...log:
fft(apply(mvfft(t(cbind(1-p,p,0,0,0))), 1, prod), inverse = TRUE)/5
Yes, that's neat (in either version). With the example I have,
where length(p)=16, I applied your suggestion above in the form
v-fft(apply(mvfft(t(cbind(1-p,p,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))),
1,prod), inverse = TRUE)/17
(making the number of columns for cbind equal to 17 = 16+1). Now,
comparing this result with the P I got by brute force (subsets
selection method), and removing the imaginary parts from v:
cbind(r=(0:16),v=Re(v),P)
rvP
0 1.353353e-01 1.353353e-01
1 3.007903e-01 3.007903e-01
2 2.976007e-01 2.976007e-01
3 1.747074e-01 1.747074e-01
4 6.826971e-02 6.826971e-02
5 1.884969e-02 1.884969e-02
6 3.804371e-03 3.804371e-03
7 5.720398e-04 5.720398e-04
8 6.463945e-05 6.463945e-05
9 5.489454e-06 5.489454e-06
10 3.473997e-07 3.473997e-07
11 1.607822e-08 1.607822e-08
12 5.262532e-10 5.262532e-10
13 1.148620e-11 1.148618e-11
14 1.494031e-13 1.493761e-13
15 8.887779e-16 8.764298e-16
16 1.434973e-17 1.313261e-19
so this calculation gives a better result than convolve().
The only fly in the ointment (which comes down to how one
sets up the arguments to cbind(), so is quite easily handled
in general application) is the variable number of columns
required according to the length of p.
Your suggestion worked nicely!
Thanks,
ted.
On 10/7/06, Gabor Grothendieck [EMAIL PROTECTED] wrote:
One can get a one-line solution by taking the product of the FFTs.
For example, let p - 1:4/8 be the probabilities. Then the solution
is:
fft(exp(rowSums(log(mvfft(t(cbind(1-p,p,0,0,0)), inverse = TRUE)/5
On 10/7/06, Ted Harding [EMAIL PROTECTED] wrote:
Hi again.
I had suspected that doing the calculation by a convolution
method might be both straightforward and efficient in R.
I've now located convolve() (in 'base'!!), and have a solution
using this function. Details below. Original statement of
problem, and method originally proposed, repeated below for
reference.
On 10/6/06, Ted Harding [EMAIL PROTECTED] wrote:
Hi Folks,
Given a series of n independent Bernoulli trials with
outcomes Yi (i=1...n) and Prob[Yi = 1] = Pi, I want
P = Prob[sum(Yi) = r] (r = 0,1,...,n)
I can certainly find a way to do it:
Let p be the vector c(P1,P2,...,Pn).
The cases r=0 and r=n are trivial (and also are exceptions
for the following routine).
For a given value of r in (1:(n-1)),
library(combinat)
Set - (1:n)
Subsets - combn(Set,r)
P - 0
for(i in (1:dim(Subsets)[2])){
ix - numeric(n)
ix[Subsets[,i]] - 1
P - P + prod((p^ix) * ((1-p)^(1-ix)))
}
The convolution method implemented in convolve computes, for
two series X[1],X[2],...,X[n] and Y[1],Y[2],...,Y[m] the
series (k = 1, 2, ... , n+m-1):
Z[k] = sum( X[k-m+i]*Y[i] )
over valid values of i, i.e. products of terms of X matched
with terms of a shifted Y, using the open method (see
?convolve).
This is not quite what is wanted for the type of convolution
needed to compute the distribution of the sum of two integer
random variables, since one needs to reverse one of the series
being convolved. This is the basis of the following:
Given a vector p of probabilities P1, P2, P3, for Y=1 in
successive trials, I need to convolve the successive Bernoulli
distributions (1-P1,P1), (1-P2,P2), ... , (1-Pn,Pn). Hence:
ps-cbind(1-p,p); u-ps[1,]; u1-u
for(i in (2:16)){
u-convolve(u,rev(ps[i,]),type=o)
}
In the case I was looking at, I had already obtained the vector
P of probabilities P(Sum(Y)=0), P(sum(Y)=1), ... by the method
previously described (repeated above); and here n=16.
Having obtained u by the above routine, I can now compare the
results, u with P:
cbind((0:16),u,P)
r uP
0 1.353353e-01 1.353353e-01
1 3.007903e-01 3.007903e-01
2 2.976007e-01 2.976007e-01
3 1.747074e-01 1.747074e-01
4 6.826971e-02 6.826971e-02
5 1.884969e-02 1.884969e-02
6 3.804371e-03 3.804371e-03
7 5.720398e-04 5.720398e-04
8 6.463945e-05 6.463945e-05
9 5.489454e-06 5.489454e-06
10 3.473997e-07 3.473997e-07
11 1.607822e-08 1.607822e-08
12 5.262533e-10 5.262532e-10
13 1.148626e-11 1.148618e-11
14 1.494650e-13 1.493761e-13
15 9.008700e-16 8.764298e-16
16 -1.896716e-17 1.313261e-19
so, apart from ignorable differences in the very small
probabilities for r11, they agree. I have to admit, though,
that I had to tread carefully (and experiment with Binomials)
to make sure of exactly how to introduce the reversal
(at rev(ps[i,])).
And the convolution method is *very* much faster than the
selection of all subsets method!
Re: [R] is it possible to fill with a color or transparency gradient?
Eric Harley wrote: Hi all, Is there a way to fill a rectangle or polygon with a color and/or transparency gradient? This would be extremely useful for me in terms of adding some additional information to some plots I'm making, especially if I could define the gradient on my own by putting functions into rgb something like rgb( r=f(x,y), g=f(x,y), b=f(x,y), alpha=f(x,y) ). Not so important whether the coordinates are in terms of the plot axes or normalized to the polygon itself somehow. Ideally it would work not only for a fill color but also for shading lines. I haven't been using R very long, so it's possible that I'm just missing something, but I haven't found anything like this in the help files. I've tried to poke around in graphics, grid, and ggplot, without any luck so far. I really like some of the functionality in ggplot, and it does some nice things with continuous gradients for the color of scatter plot points, for example, but it each individual point (or grob) is always one solid color as far as I can tell. Hi Eric, You can fill a rectangle with a gradient (slices of colors) using gradient.rect in the plotrix package. Filling an arbitrary polygon would be a lot harder to program. You may also be interested in the color.scale function that assigns colors to numerical values. Jim __ 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] Installing Lindsey's packages
Hi Ken, My message was not well-focused. Here's the real problem: mk% R CMD INSTALL rmutil.tgz ERROR: cannot extract package from 'rmutil.tgz' On Oct 8, 2006, at 5:23 AM, Ken Knoblauch wrote: Hi Michael, What works for me is to retrieve the tgz from the garbage and to compile it from a terminal command line with R CMD INSTALL, this is assuming that you have the Development tools installed (of course). Hope that gets you there. best, Ken ** Dear r-helpers, I downloaded http://popgen.unimaas.nl/~jlindsey/rcode/rmutil.tar (it was originally .tgz, but got unzipped by my browser). Can anyone give me detailed instructions on installing this and Lindsey's other packages on R version 2.4.0 (2006-10-03)---(powerpc- apple-darwin8.7.0, locale: C)? -- Ken Knoblauch Inserm U371 Institut Cellule Souche et Cerveau Dept. of Cognitive Neuroscience 18 avenue du Doyen Lépine 69500 Bron France tel: +33 (0)4 72 91 34 77 fax: +33 (0)4 72 91 34 61 portable: +33 (0)6 84 10 64 10 http://www.lyon.inserm.fr/371/ _ Professor Michael Kubovy University of Virginia Department of Psychology USPS: P.O.Box 400400Charlottesville, VA 22904-4400 Parcels:Room 102Gilmer Hall McCormick RoadCharlottesville, VA 22903 Office:B011+1-434-982-4729 Lab:B019+1-434-982-4751 Fax:+1-434-982-4766 WWW:http://www.people.virginia.edu/~mk9y/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] [R-pkg] New packages pmg, gWidgets, gWidgetsRGtk2
hi, i got the same error messages. but after installing GTK libraries AND a restart of my windows machine, the package pmg and his GUI was starting correctly. greetings from austria helli Message: 19 Date: Fri, 6 Oct 2006 12:26:00 -0400 From: Gabor Grothendieck [EMAIL PROTECTED] Subject: Re: [R] [R-pkg] New packages pmg, gWidgets, gWidgetsRGtk2 To: j verzani [EMAIL PROTECTED] Cc: r-help@stat.math.ethz.ch Message-ID: [EMAIL PROTECTED] Content-Type: text/plain; charset=ISO-8859-1; format=flowed Hi, I need installation instructions. library(pmg) seems not to be enough. Thanks. library(pmg) Loading pmg() Loading required package: gWidgets Loading required package: gWidgetsRGtk2 Loading required package: RGtk2 Error in dyn.load(x, as.logical(local), as.logical(now)) : unable to load shared library 'C:/PROGRA~1/R/R-24~1.0/library/RGtk2/libs/RGtk2.dll': LoadLibrary failure: The specified procedure could not be found. Error: package 'RGtk2' could not be loaded Error: .onAttach failed in 'attachNamespace' Loading required package: pmggWidgetsRGtk PMG needs gWidgetsError in assign(x, value, envir = ns, inherits = FALSE) : could not find function gwindow In addition: Warning message: there is no package called 'pmggWidgetsRGtk' in: library(package, lib.loc = lib.loc, character.only = TRUE, logical = TRUE, Error: .onLoad failed in 'loadNamespace' for 'pmg' Error: package/namespace load failed for 'pmg' On 10/6/06, j verzani [EMAIL PROTECTED] wrote: I'd like to announce three new packages on CRAN: pmg, gWidgets, and __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] [R-pkg] New packages pmg, gWidgets, gWidgetsRGtk2
[EMAIL PROTECTED] writes:
Finally, for some reason I'm trying to track down, under Windows, pmg
is crashing when a cairoDevice graphics device is being closed. This
affects the plotnotebook and the Lattice Explorer. I have an update to
gWidgetsRGtk2 that may fix it, but I haven't had a chance to fully
test it. It can be installed with
This appears not to be Windows-specific. I see it on Fedora Core 5
too:
plot(x = ashina$vas.active, y = ashina$vas.plac, type = p,
+ pch = 1, col = ashina$grp)
*** caught segfault ***
address 0xfc, cause 'memory not mapped'
Traceback:
1: dev.off([EMAIL PROTECTED])
2: handler(...)
3: function (...) {handler(...)return(TRUE)}(list(obj = S4
object of class gGraphicsRGtk, action = NULL), pointer:
0xc327a60)
Possible actions:
1: abort (with core dump)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace
(2.4.0 Patched, upon closing the plotnotebook)
--
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] Simulate p-value in lme4
Dear r-helpers,
Spencer Graves and Manual Morales proposed the following methods to
simulate p-values in lme4:
preliminary
require(lme4)
require(MASS)
summary(glm(y ~ lbase*trt + lage + V4, family = poisson, data =
epil), cor = FALSE)
epil2 - epil[epil$period == 1, ]
epil2[period] - rep(0, 59); epil2[y] - epil2[base]
epil[time] - 1; epil2[time] - 4
epil2 - rbind(epil, epil2)
epil2$pred - unclass(epil2$trt) * (epil2$period 0)
epil2$subject - factor(epil2$subject)
epil3 - aggregate(epil2, list(epil2$subject, epil2$period 0),
function(x) if(is.numeric(x)) sum(x) else x[1])
simulation (SG)
o - with(epil3, order(subject, y))
epil3. - epil3[o,]
norep - with(epil3., subject[-1]!=subject[-dim(epil3)[1]])
subj1 - which(c(T, norep))
subj.pred - epil3.[subj1, c(subject, pred)]
subj. - as.character(subj.pred$subject)
pred. - subj.pred$pred
names(pred.) - subj.
iter - 10
chisq.sim - rep(NA, iter)
set.seed(1)
for(i in 1:iter){
s.i - sample(subj.)
# Randomize subject assignments to 'pred' groups
epil3.$pred - pred.[s.i][epil3.$subject]
fit1 - lmer(y ~ pred+(1 | subject),
family = poisson, data = epil3.)
fit0 - lmer(y ~ 1+(1 | subject),
family = poisson, data = epil3.)
chisq.sim[i] - anova(fit0, fit1)[2, Chisq]
}
simulation (MM)
iter - 10
chisq.sim - rep(NA, iter)
Zt - slot(model1,Zt) # see ?lmer-class
n.grps - dim(ranef(model1)[[1]])[1]
sd.ran.effs - sqrt(VarCorr(model1)[[1]][1])
X - slot(model1,X) # see ?lmer-class
fix.effs - matrix(rep(fixef(model1),dim(X)[1]), nrow=dim(X)[1],
byrow=T)
model.parms - X*fix.effs # This gives parameters for each case
# Generate predicted values
pred.vals - as.vector(apply(model.parms, 1, sum))
for(i in 1:iter) {
rand.new - as.vector(rnorm(grps,0, sd.ran.effs)) #grps should be
n.grps?
rand.vals - as.vector(rand.new%*%Zt) # Assign random effects
mu - pred.vals+rand.vals # Expected mean
resp - rpois(length(mu), exp(mu))
sim.data - data.frame(slot(model2,frame), resp) # Make data frame
sim.model1 - lmer(resp~1+(1|subject), data=sim.data,
family=poisson)
sim.model2 - lmer(resp~pred+(1|subject), data=sim.data,
family=poisson)
chisq.sim[i] - anova(sim.model1,sim.model2)$Chisq[[2]]
}
end
Unfortunately the latter fails (even if I replace grps with n.grps):
Error in slot(model2, frame) : object model2 not found
In any event, I would be eager to hear more discussion of the pros
and cons of these approaches.
_
Professor Michael Kubovy
University of Virginia
Department of Psychology
USPS: P.O.Box 400400Charlottesville, VA 22904-4400
Parcels:Room 102Gilmer Hall
McCormick RoadCharlottesville, VA 22903
Office:B011+1-434-982-4729
Lab:B019+1-434-982-4751
Fax:+1-434-982-4766
WWW:http://www.people.virginia.edu/~mk9y/
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Error Correcting Codes, Simplex
Hello, Is there a way in R to construct an (error correcting) binary code e.g. for an source alphabet containing integers from 1 to say 255 with the property that each pair of distinct codewords of length m is at Hamming distance exactly m/2 ? I was suggested to use so called simplex codes, which should be fairly standard, but I haven't found a direct way via R packages to do so, that's why I ask whether there might be in indirect way to solve this problem. Example: v1 =c(1,2,3,4) v2 =c(1,2,5,6) similarity(v1,v2)=0.5, (because 2 out of 4 elements are equal). Obviously, a binary representation of would yield a different similarity of: binary(v1) =001 010 011 100 binary(v1) =001 010 101 110 similarity(binary(v1),binary(v2))= 9/12 Remark: The focus here is not on error correction, but rather the binary encoding retaining similarity of the elements of vectors. Many thanks, Bjoern [[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] Sum of Bernoullis with varying probabilities
On 10/8/06, Ted Harding [EMAIL PROTECTED] wrote: On 08-Oct-06 Gabor Grothendieck wrote: Or perhaps its clearer (and saves a bit of space) to use apply...prod here instead of exp...log: fft(apply(mvfft(t(cbind(1-p,p,0,0,0))), 1, prod), inverse = TRUE)/5 Yes, that's neat (in either version). With the example I have, where length(p)=16, I applied your suggestion above in the form v-fft(apply(mvfft(t(cbind(1-p,p,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))), 1,prod), inverse = TRUE)/17 (making the number of columns for cbind equal to 17 = 16+1). Now, comparing this result with the P I got by brute force (subsets selection method), and removing the imaginary parts from v: cbind(r=(0:16),v=Re(v),P) rvP 0 1.353353e-01 1.353353e-01 1 3.007903e-01 3.007903e-01 2 2.976007e-01 2.976007e-01 3 1.747074e-01 1.747074e-01 4 6.826971e-02 6.826971e-02 5 1.884969e-02 1.884969e-02 6 3.804371e-03 3.804371e-03 7 5.720398e-04 5.720398e-04 8 6.463945e-05 6.463945e-05 9 5.489454e-06 5.489454e-06 10 3.473997e-07 3.473997e-07 11 1.607822e-08 1.607822e-08 12 5.262532e-10 5.262532e-10 13 1.148620e-11 1.148618e-11 14 1.494031e-13 1.493761e-13 15 8.887779e-16 8.764298e-16 16 1.434973e-17 1.313261e-19 so this calculation gives a better result than convolve(). The only fly in the ointment (which comes down to how one sets up the arguments to cbind(), so is quite easily handled in general application) is the variable number of columns required according to the length of p. Try this: p - seq(4)/8 # test data Re(fft(apply(mvfft(t(cbind(1-p, p, 0*p%o%p[-1]))), 1, prod), inv = TRUE)/(length(p)+1)) __ 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] Sum of Bernoullis with varying probabilities
On 10/8/06, Gabor Grothendieck [EMAIL PROTECTED] wrote: On 10/8/06, Ted Harding [EMAIL PROTECTED] wrote: On 08-Oct-06 Gabor Grothendieck wrote: Or perhaps its clearer (and saves a bit of space) to use apply...prod here instead of exp...log: fft(apply(mvfft(t(cbind(1-p,p,0,0,0))), 1, prod), inverse = TRUE)/5 Yes, that's neat (in either version). With the example I have, where length(p)=16, I applied your suggestion above in the form v-fft(apply(mvfft(t(cbind(1-p,p,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))), 1,prod), inverse = TRUE)/17 (making the number of columns for cbind equal to 17 = 16+1). Now, comparing this result with the P I got by brute force (subsets selection method), and removing the imaginary parts from v: cbind(r=(0:16),v=Re(v),P) rvP 0 1.353353e-01 1.353353e-01 1 3.007903e-01 3.007903e-01 2 2.976007e-01 2.976007e-01 3 1.747074e-01 1.747074e-01 4 6.826971e-02 6.826971e-02 5 1.884969e-02 1.884969e-02 6 3.804371e-03 3.804371e-03 7 5.720398e-04 5.720398e-04 8 6.463945e-05 6.463945e-05 9 5.489454e-06 5.489454e-06 10 3.473997e-07 3.473997e-07 11 1.607822e-08 1.607822e-08 12 5.262532e-10 5.262532e-10 13 1.148620e-11 1.148618e-11 14 1.494031e-13 1.493761e-13 15 8.887779e-16 8.764298e-16 16 1.434973e-17 1.313261e-19 so this calculation gives a better result than convolve(). The only fly in the ointment (which comes down to how one sets up the arguments to cbind(), so is quite easily handled in general application) is the variable number of columns required according to the length of p. Try this: p - seq(4)/8 # test data Re(fft(apply(mvfft(t(cbind(1-p, p, 0*p%o%p[-1]))), 1, prod), inv = TRUE)/(length(p)+1)) And here is one minor improvement (rbind() rather than t(cbind()): p - 1:4/8 Re(fft(apply(mvfft(rbind(1-p, p, 0 * p[-1] %o% p)), 1, prod), inverse = TRUE) / (length(p) + 1)) __ 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] Sum of Bernoullis with varying probabilities
On 08-Oct-06 Gabor Grothendieck wrote: And here is one minor improvement (rbind() rather than t(cbind()): p - 1:4/8 Re(fft(apply(mvfft(rbind(1-p, p, 0 * p[-1] %o% p)), 1, prod), inverse = TRUE) / (length(p) + 1)) Really nice! Much better than any solution I might have found for myself. It really is worth posting to the list! This all still leaves me with one nagging question in the back of my mind -- Surely this is not the first time someone has tackled the problem (in R) of the distribution of the sum of Bernoulli RVs with unequal probabilities? R Site Search on Bernoulli did not seem to throw anything up. (Or maybe they just did it, and didn't think it worth saying anything about--but it's not really trivial). Best wishes, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 08-Oct-06 Time: 14:36:43 -- XFMail -- __ 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] Windows/MAC difference (console)
Hello, a colleague of mine uses R on his Mac and he has quite a nice feature: When he starts writing a part of the function like plot( the programming is showing him at the bottom the kind of arguments you can input to the function. He told me that he has not installed any further stuff than the base R program. Is this a feature you only have in the MAC version? I really like it. I tried some editors, but there were all not the thing I would like to have. I am only searching for such a display and perhaps syntax highlighting in the console. But I do not want to have all this GUI stuff like you have in R commander. Do you have any suggestion? Thanks in advance, lisra [[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] Select range of dates
And here is a fourth: as.numeric(format(dd, %Y)) + (quarters(dd) Q1) On 10/8/06, Gabor Grothendieck [EMAIL PROTECTED] wrote: Here are three alternative ways to get the fiscal year as a numeric value assuming: dd - as.Date(x$Date,%d/%m/%Y) # add one to year if month is past March as.numeric(format(dd, %Y)) + (format(dd, %m) 03) # same but using POSIXlt # (Even though there are no time zones involved I have seen # situations where the time zone nevertheless crept in # where one would least expect it so even though I believe # this is correct I would nevertheless triple check # all your results if you use this one) with(as.POSIXlt(dd), 1900 + year + (mon 2)) # this makes use of the zoo yearqtr class # which represents dates as year + 0, .25, .5 or .75 for # the four quarters library(zoo) as.numeric(ceiling(as.yearqtr(dd))) On 10/7/06, Ian Broom [EMAIL PROTECTED] wrote: Hello, This is likely fairly silly question, and I apologize to whomever takes the time to respond. I am a relatively new user of R, on Windows XP, version 2.3.1. Say I have a data table that looks like the following: x Date Location Amount Blue Green 1 01/01/2001 Central1817 TRUE FALSE 2 01/02/2001 Central 20358 FALSE TRUE 3 05/08/2001 Central 16245 FALSE TRUE 4 02/02/2002 Western 112 TRUE FALSE 5 21/03/2002 Western 98756 TRUE FALSE 6 01/04/2002 Western 1598414 FALSE TRUE 7 07/01/2001 Western1255 FALSE TRUE 8 20/10/2003 Central 16289 TRUE FALSE 9 21/10/2003 Eastern 1 FALSE TRUE 10 22/10/2003 Eastern 98737 FALSE TRUE 11 23/10/2003 Eastern 198756 TRUE FALSE 12 24/10/2003 Eastern 98756 FALSE TRUE 13 25/10/2003 Eastern 65895 TRUE FALSE 14 26/10/2003 Eastern 2142266 FALSE TRUE 15 27/10/2003North 98756 TRUE FALSE 16 28/10/2003North 548236 FALSE TRUE and I want to do some summaries by Fiscal year (or FY quarter). Reading manuals and such, I cobbled this less than satisfactory bit together to start to build a factor representing a fiscal year split: y-as.Date(x$Date,%d/%m/%Y) y-as.matrix(y) y$FY0203-ifelse((y=(as.Date(2002-03-21)))(y=(as.Date (2003-04-01))),TRUE,FALSE) the values seem correct, but aside from ugly, the data is not in the proper format - with some more effort I might fix that... But my question is: is there a more simple function available to select a range of dates? Or, at least, a more elegant approach? [[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.
Re: [R] Windows/MAC difference (console)
The Tinn-R editor on sourceforge can do that. There is also a useful intro here on setting up Tinn-R: http://genetics.agrsci.dk/~sorenh/misc/Rlive/index.html On 10/8/06, Lina Jansen [EMAIL PROTECTED] wrote: Hello, a colleague of mine uses R on his Mac and he has quite a nice feature: When he starts writing a part of the function like plot( the programming is showing him at the bottom the kind of arguments you can input to the function. He told me that he has not installed any further stuff than the base R program. Is this a feature you only have in the MAC version? I really like it. I tried some editors, but there were all not the thing I would like to have. I am only searching for such a display and perhaps syntax highlighting in the console. But I do not want to have all this GUI stuff like you have in R commander. Do you have any suggestion? Thanks in advance, lisra [[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] [R-pkgs] new version IPSUR_0.1-1 now available on CRAN
Dear All, The latest version of IPSUR (last one for awhile) is now available on CRAN. What's New in 0.1-1: 1. Important bug fixes 2. Test of equality for several proportions... 3. Enter table for single proportion test... 4. Enter table for multiple proportions test... 5. Bar graph plot by groups (stacked and side-by side) 6. Enter table for bar graphs... Visit the IPSUR website for the latest bug fixes, updates, and new features: http://www.cc.ysu.edu/~gjkerns/IPSUR/packagehttp://www.cc.ysu.edu/%7Egjkerns/IPSUR/package Regards, Jay [[alternative HTML version deleted]] ___ R-packages mailing list R-packages@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-packages __ 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] latex and anova.lme problem
Dear R-helpers, When I try anova(txtE2.lme, txtE2.lme1) Model df AIC BIC logLik Test L.Ratio p-value txtE2.lme 1 10 8590 8638 -4285 txtE2.lme1 2 7 8591 8624 -4288 1 vs 26.79 0.0789 latex(anova(txtE2.lme, txtE2.lme1)) Error: object n.group not found I don't even see n.group as one of the arguments of latex() I checked to see class(anova(txtE2.lme, txtE2.lme1)) [1] anova.lme data.frame A bit more information (which I don't know is relevant): methods(anova.lme) no methods were found Warning message: function 'anova.lme' appears not to be generic in: methods(anova.lme) methods(latex) [1] latex.bystats latex.bystats2 latex.default [4] latex.describe latex.describe.single latex.function [7] latex.list latex.summary.formula.cross latex.summary.formula.response [10] latex.summary.formula.reverse *sessionInfo() R version 2.4.0 (2006-10-03) powerpc-apple-darwin8.7.0 locale: C attached base packages: [1] datasets methods stats graphics grDevices utils base other attached packages: Hmisc chron xtablegeepack glmmML nlme lme4 Matrixlattice 3.1-12.3-81.3-2 1.0-10 0.65-3 3.1-77 0.9975-1 0.9975-2 0.14-9 MASSJGR iplots JavaGD rJava 7.2-29 1.4-111.0-40.3-5 0.4-10 _ Professor Michael Kubovy University of Virginia Department of Psychology USPS: P.O.Box 400400Charlottesville, VA 22904-4400 Parcels:Room 102Gilmer Hall McCormick RoadCharlottesville, VA 22903 Office:B011+1-434-982-4729 Lab:B019+1-434-982-4751 Fax:+1-434-982-4766 WWW:http://www.people.virginia.edu/~mk9y/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] error in save.image
Greetings R-ians: Not infrequently of late I get this error message when trying to save what I'm working on. Error in save.image(C:/Documents and Settings/ ... /.RData) : image could not be renamed and is left in C:/Documents and Settings/ ... /.RDataTmp10 If I try again, it almost always is successful. I don't recall having had this situation earlier, but it does seem more frequent lately. Have others observed this? Could it be because my workspace has become large? (Large being relative: 660KB) I'm using Version 2.3.1 (2006-06-01) running Windows XP on a Dell with 2Meg. Thanks. Charles Annis, P.E. PS - Yes, I have upgraded to R version 2.4.0 (2006-10-03) but do not want to change in the middle of a project. [EMAIL PROTECTED] phone: 561-352-9699 eFax: 614-455-3265 http://www.StatisticalEngineering.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] Problem in getting 632plus error using randomForest by ipred!
Hello! I'm Taeho, a graduate student in South Korea. In order to get .632+ bootstrap error using random forest, I have tried to use 'ipred' package; more specifically the function 'errorest' has been used. Following the guidelines, I made a simple command line like below: error-errorest(class ~ ., data=data, model=randomForest, estimator = 632plus)$err however, I got an error message saying that: randomForest.default(m, y, ...) : Can't have empty classes in y. The data I used includes three classes and each of the classes has 5, 2, and 4 samples, respectively. Other than this fact about the data, nothing is different from the example 'iris' data. (By the way, I could get a simple OOB error for this data without any problem just using 'randomForest' function only.) What do you think my problem is? I think the problem is that the number of samples is too small to run 'errorest' for getting .632 plus bootstrap error. Is this correct? If so, what would you recommend to fix this problem? Thank you for reading this and for your kind answer in advance. 새로운 기부 문화의 씨앗, 해피빈 http://happybean.naver.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] latex and anova.lme problem
Michael Kubovy wrote: Dear R-helpers, When I try anova(txtE2.lme, txtE2.lme1) Model df AIC BIC logLik Test L.Ratio p-value txtE2.lme 1 10 8590 8638 -4285 txtE2.lme1 2 7 8591 8624 -4288 1 vs 26.79 0.0789 latex(anova(txtE2.lme, txtE2.lme1)) Error: object n.group not found Michael, As far as I know, no one has implemented a latex method for anova.lme objects. You are trying to use a default latex method for a complex object. n.group is an argument to latex.default in the Hmisc package Frank I don't even see n.group as one of the arguments of latex() I checked to see class(anova(txtE2.lme, txtE2.lme1)) [1] anova.lme data.frame A bit more information (which I don't know is relevant): methods(anova.lme) no methods were found Warning message: function 'anova.lme' appears not to be generic in: methods(anova.lme) methods(latex) [1] latex.bystats latex.bystats2 latex.default [4] latex.describe latex.describe.single latex.function [7] latex.list latex.summary.formula.cross latex.summary.formula.response [10] latex.summary.formula.reverse *sessionInfo() R version 2.4.0 (2006-10-03) powerpc-apple-darwin8.7.0 locale: C attached base packages: [1] datasets methods stats graphics grDevices utils base other attached packages: Hmisc chron xtablegeepack glmmML nlme lme4 Matrixlattice 3.1-12.3-81.3-2 1.0-10 0.65-3 3.1-77 0.9975-1 0.9975-2 0.14-9 MASSJGR iplots JavaGD rJava 7.2-29 1.4-111.0-40.3-5 0.4-10 _ Professor Michael Kubovy University of Virginia Department of Psychology USPS: P.O.Box 400400Charlottesville, VA 22904-4400 Parcels:Room 102Gilmer Hall McCormick RoadCharlottesville, VA 22903 Office:B011+1-434-982-4729 Lab:B019+1-434-982-4751 Fax:+1-434-982-4766 WWW:http://www.people.virginia.edu/~mk9y/ __ 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. -- Frank E Harrell Jr Professor and Chair School of Medicine Department of Biostatistics Vanderbilt 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] 'weaver' package problem
Hi Seth, The possibility of caching computations would be a great boon when one is iteratively refining a paper; so I'm most grateful for your work on this. Unfortunately I have a problem to report: **installing** source(http://bioconductor.org/biocLite.R;) biocLite(weaver) Running getBioC version 0.1.8 with R version 2.4.0 Running biocinstall version 1.9.8 with R version 2.4.0 Your version of R requires version 1.9 of Bioconductor. trying URL 'http://bioconductor.org/packages/1.9/bioc/bin/macosx/ universal/contrib/2.4/weaver_1.0.0.tgz' Content type 'application/x-gzip' length 101153 bytes opened URL == downloaded 98Kb The downloaded packages are in /tmp/RtmpO6uHHW/downloaded_packages **loading weaver** library(weaver) Loading required package: digest Loading required package: tools Loading required package: codetools **testing Sweave** testfile - system.file(Sweave, Sweave-test-1.Rnw, package = utils) ## enforce par(ask=FALSE) options(par.ask.default=FALSE) ## create a LaTeX file Sweave(testfile) Writing to file Sweave-test-1.tex Processing code chunks ... 1 : print term verbatim 2 : term hide 3 : echo print term verbatim 4 : term verbatim 5 : echo term verbatim 6 : echo term verbatim eps pdf 7 : echo term verbatim eps pdf You can now run LaTeX on 'Sweave-test-1.tex' **testing weaver** weaver(testfile) Error in weaver(testfile) : unused argument(s) (/Library/Frameworks/ R.framework/Resources/library/utils/Sweave/Sweave-test-1.Rnw) **session info** R version 2.4.0 (2006-10-03) powerpc-apple-darwin8.7.0 locale: C attached base packages: [1] tools grid splines datasets methods stats graphics grDevices [9] utils base other attached packages: weaver codetools digesteffects Hmisc chron xtablegeepack glmmML 1.0.00.0-30.2.21.0-93.1-12.3-8 1.3-2 1.0-10 0.65-3 nlme lme4 Matrixlattice MASS JGR iplots JavaGD rJava 3.1-77 0.9975-1 0.9975-2 0.14-9 7.2-29 1.4-11 1.0-40.3-5 0.4-10 _ Professor Michael Kubovy University of Virginia Department of Psychology USPS: P.O.Box 400400Charlottesville, VA 22904-4400 Parcels:Room 102Gilmer Hall McCormick RoadCharlottesville, VA 22903 Office:B011+1-434-982-4729 Lab:B019+1-434-982-4751 Fax:+1-434-982-4766 WWW:http://www.people.virginia.edu/~mk9y/ __ 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] 'weaver' package problem
Hi Michael, Michael Kubovy [EMAIL PROTECTED] writes: Hi Seth, The possibility of caching computations would be a great boon when one is iteratively refining a paper; so I'm most grateful for your work on this. Unfortunately I have a problem to report: **testing Sweave** testfile - system.file(Sweave, Sweave-test-1.Rnw, package = utils) **testing weaver** weaver(testfile) Error in weaver(testfile) : unused argument(s) (/Library/Frameworks/ R.framework/Resources/library/utils/Sweave/Sweave-test-1.Rnw) Actually, you want: Sweave(testfile, driver=weaver()) I will add something to the man page for weaver to make this more clear. You might want to have a look at the vignette included with the package: library(Biobase) library(weaver) openVignette(weaver) + seth __ 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] latex and anova.lme problem
It works for me. I am also using R-2.4.0 library(Hmisc) library(nlme) ?lme fm1 - lme(distance ~ age, data = Orthodont) # random is ~ age fm2 - lme(distance ~ age + Sex, data = Orthodont, random = ~ 1) anova(fm1,fm2) Model df AIC BIClogLik Test L.Ratio p-value fm1 1 6 454.6367 470.6173 -221.3183 fm2 2 5 447.5125 460.7823 -218.7562 1 vs 2 5.124178 0.0236 Warning message: Fitted objects with different fixed effects. REML comparisons are not meaningful. in: anova.lme (fm1, fm2) fm1f - lme(distance ~ age+Sex, data = Orthodont) ## to avoid above warning message anova(fm1f,fm2) Model df AIC BIClogLik Test L.Ratio p-value fm1f 1 10 453.6604 480.2000 -216.8302 fm2 2 5 447.5125 460.7823 -218.7562 1 vs 2 3.852152 0.5709 tmp.lme - latex(anova(fm1f,fm2)) getwd() ## to find out where the anova.tex file was placed. Therefore it is necessary for you to provide commented, minimal, self-contained, reproducible code. 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.
Re: [R] 'weaver' package problem
Michael Kubovy [EMAIL PROTECTED] writes: Hi Seth, The possibility of caching computations would be a great boon when one is iteratively refining a paper; so I'm most grateful for your work on this. Unfortunately I have a problem to report: **installing** source(http://bioconductor.org/biocLite.R;) biocLite(weaver) Running getBioC version 0.1.8 with R version 2.4.0 Running biocinstall version 1.9.8 with R version 2.4.0 Your version of R requires version 1.9 of Bioconductor. trying URL 'http://bioconductor.org/packages/1.9/bioc/bin/macosx/ universal/contrib/2.4/weaver_1.0.0.tgz' Content type 'application/x-gzip' length 101153 bytes opened URL == downloaded 98Kb The downloaded packages are in /tmp/RtmpO6uHHW/downloaded_packages **loading weaver** library(weaver) Loading required package: digest Loading required package: tools Loading required package: codetools **testing Sweave** testfile - system.file(Sweave, Sweave-test-1.Rnw, package = utils) ## enforce par(ask=FALSE) options(par.ask.default=FALSE) ## create a LaTeX file Sweave(testfile) Writing to file Sweave-test-1.tex Processing code chunks ... 1 : print term verbatim 2 : term hide 3 : echo print term verbatim 4 : term verbatim 5 : echo term verbatim 6 : echo term verbatim eps pdf 7 : echo term verbatim eps pdf You can now run LaTeX on 'Sweave-test-1.tex' **testing weaver** weaver(testfile) Error in weaver(testfile) : unused argument(s) (/Library/Frameworks/ R.framework/Resources/library/utils/Sweave/Sweave-test-1.Rnw) Er, that would not appear to be the intended usage, as per http://www.bioconductor.org/packages/1.9/bioc/vignettes/weaver/inst/doc/weaver_howTo.pdf#search=%22weaver%20bioconductor%22 -- 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] Simulate p-value in lme4
On Sun, 2006-10-08 at 07:34 -0400, Michael Kubovy wrote:
Dear r-helpers,
Spencer Graves and Manual Morales proposed the following methods to
simulate p-values in lme4:
preliminary
require(lme4)
require(MASS)
summary(glm(y ~ lbase*trt + lage + V4, family = poisson, data =
epil), cor = FALSE)
epil2 - epil[epil$period == 1, ]
epil2[period] - rep(0, 59); epil2[y] - epil2[base]
epil[time] - 1; epil2[time] - 4
epil2 - rbind(epil, epil2)
epil2$pred - unclass(epil2$trt) * (epil2$period 0)
epil2$subject - factor(epil2$subject)
epil3 - aggregate(epil2, list(epil2$subject, epil2$period 0),
function(x) if(is.numeric(x)) sum(x) else x[1])
simulation (SG)
o - with(epil3, order(subject, y))
epil3. - epil3[o,]
norep - with(epil3., subject[-1]!=subject[-dim(epil3)[1]])
subj1 - which(c(T, norep))
subj.pred - epil3.[subj1, c(subject, pred)]
subj. - as.character(subj.pred$subject)
pred. - subj.pred$pred
names(pred.) - subj.
iter - 10
chisq.sim - rep(NA, iter)
set.seed(1)
for(i in 1:iter){
s.i - sample(subj.)
# Randomize subject assignments to 'pred' groups
epil3.$pred - pred.[s.i][epil3.$subject]
fit1 - lmer(y ~ pred+(1 | subject),
family = poisson, data = epil3.)
fit0 - lmer(y ~ 1+(1 | subject),
family = poisson, data = epil3.)
chisq.sim[i] - anova(fit0, fit1)[2, Chisq]
}
simulation (MM)
iter - 10
chisq.sim - rep(NA, iter)
Zt - slot(model1,Zt) # see ?lmer-class
n.grps - dim(ranef(model1)[[1]])[1]
sd.ran.effs - sqrt(VarCorr(model1)[[1]][1])
X - slot(model1,X) # see ?lmer-class
fix.effs - matrix(rep(fixef(model1),dim(X)[1]), nrow=dim(X)[1],
byrow=T)
model.parms - X*fix.effs # This gives parameters for each case
# Generate predicted values
pred.vals - as.vector(apply(model.parms, 1, sum))
for(i in 1:iter) {
rand.new - as.vector(rnorm(grps,0, sd.ran.effs)) #grps should be
#n.grps?
Yes. Change grps to n.grps.
rand.vals - as.vector(rand.new%*%Zt) # Assign random effects
mu - pred.vals+rand.vals # Expected mean
resp - rpois(length(mu), exp(mu))
sim.data - data.frame(slot(model2,frame), resp) # Make data frame
sim.model1 - lmer(resp~1+(1|subject), data=sim.data,
family=poisson)
sim.model2 - lmer(resp~pred+(1|subject), data=sim.data,
family=poisson)
chisq.sim[i] - anova(sim.model1,sim.model2)$Chisq[[2]]
}
end
Unfortunately the latter fails (even if I replace grps with n.grps):
Error in slot(model2, frame) : object model2 not found
You need to specify the models fit0 and fit1 from SG's example as model1
and model2. E.g., model1 - fit0, etc.
In any event, I would be eager to hear more discussion of the pros
and cons of these approaches.
One practical problem with my approach (MM) is that the fitting
algorithms for lmer would often choke - in particular for those
randomly generated data sets that by chance included variance components
close to zero.
In any case, the MCMC approach advocated by Douglas Bates is *much*
faster. That's the approach I've been using since he suggested a
function for estimating P-values from MCMC samples for factors (or model
comparisons) with greater than 2 levels.
See:http://tolstoy.newcastle.edu.au/R/e2/help/06/09/1003.html
and the very long thread that accompanies it.
HTH,
Manuel
--
Manuel A. Morales
http://mutualism.williams.edu
signature.asc
Description: This is a digitally signed message part
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] latex and anova.lme problem
Frank E Harrell Jr f.harrell at vanderbilt.edu writes: As far as I know, no one has implemented a latex method for anova.lme objects. Well, I have, based on your latex/Hmisc, but it's a bit to my private taste, so I won't upload it to CRAN. Also has latex.glm, .lme, .summary.lme. http://www.menne-biomed.de/download/dmisc.zip http://www.menne-biomed.de/download/dmisc_0.3.tar.gz Dieter __ 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] Generating bivariate or multivariate data with known parameter values
Greetings, I'm interested in generating data from various bivariate or mulitivariate distributions (e.g. gamma, t, etc), where I can specify the parameter values, including the correlations among the variables. I haven't been able to dig anything up on the faq, but I probably missed something. A nudge in the right direction would be appreciated. David -- David Kaplan, Ph.D. Professor Department of Educational Psychology University of Wisconsin - Madison Educational Sciences, Room 1061 1025 W. Johnson Street Madison, WI 53706 email: [EMAIL PROTECTED] Web: http://www.education.wisc.edu/edpsych/facstaff/kaplan/kaplan.htm Phone: 608-262-0836 Fax: 608-262-0843 __ 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] Generating bivariate or multivariate data with known parameter values
you can do that using the copula package, e.g., library(copula) bivGamma - mvdc(claytonCopula(2), c(gamma, gamma), list(list(shape = 1, rate = 2), list(shape = 2, rate = 1))) dat - rmvdc(bivGamma, 1000) regarding correlation, the Clayton copula with alpha = 2 gives Kendall's tau = 0.5, i.e., cor(dat, method = kendall) I hope it helps. 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 David Kaplan [EMAIL PROTECTED]: Greetings, I'm interested in generating data from various bivariate or mulitivariate distributions (e.g. gamma, t, etc), where I can specify the parameter values, including the correlations among the variables. I haven't been able to dig anything up on the faq, but I probably missed something. A nudge in the right direction would be appreciated. David -- David Kaplan, Ph.D. Professor Department of Educational Psychology University of Wisconsin - Madison Educational Sciences, Room 1061 1025 W. Johnson Street Madison, WI 53706 email: [EMAIL PROTECTED] Web: http://www.education.wisc.edu/edpsych/facstaff/kaplan/kaplan.htm Phone: 608-262-0836 Fax: 608-262-0843 __ 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. 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.
[R] Probability of exceedance function question
I'm trying to calculate a cumulative area distribution (graph) of drainage areas. This is defined as P(A A*). Simple in principle. I can do this in excel, with COUNTIF, which will count the number of cells in the row area that have area A, then determine, for each cell in the row area, how many cells exceede that area, then dividing that number by the total number of cells, which gives me the probability that drainage area A exceeds drainage area A*. E.g, drainage area of 6 sq meters (One DEM grid cell) has a high probability of exceedance(.99), while a drainage area of 100,000 square meters has a low probability of exceedance (.001). I wish to plot this relationship, and we all know that excel is not the tool of choice when working with hundreds of thousands of records. I'd like to port the CAD into a few R functions that I've already developed for other tests as well. So my challenge, in R, is to (1)count the number of rows in column Area that have AREA(*), (2) determine, by row, how many rows have an area greater than the area given in that one row (3) divide step 2 by number of rows (how can I do a row count and port that to a variable, as I have to do this on 10 datasets?) Thanks for any advice you can offer to this endevour __ 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] Size problem with two dotcharts side by side
Dear all,
I'm trying to produce two dotcharts side-by-side within a Sweave
document. When I'm compiling this example:
\documentclass{article}
\begin{document}
fig=T,width=8,height=4=
par(mfrow = c(1, 2), cex = 0.7)
for(i in 1:2) dotchart(1:10)
@
fig=T,width=8,height=4=
par(mfrow = c(1, 2), cex = 0.7)
for(i in 1:2) hist(1:10)
@
fig=T,width=8,height=4=
par(mfrow = c(1, 2), cex = 0.7)
for(i in 1:2) plot(1:10)
@
\end{document}
I obtain the following:
http://pbil.univ-lyon1.fr/members/lobry/mini.pdf
that is with the second dotchart (on the right) smaller than the
one on the left. I do not have the same problem with
hist() or plot().
Any idea of what should I do for my two dotcharts to be of
the same size?
Best,
sessionInfo()
R version 2.4.0 (2006-10-03)
powerpc-apple-darwin8.7.0
locale:
fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
attached base packages:
[1] methods stats graphics grDevices utils datasets
[7] base
--
Jean R. Lobry([EMAIL PROTECTED])
Laboratoire BBE-CNRS-UMR-5558, Univ. C. Bernard - LYON I,
43 Bd 11/11/1918, F-69622 VILLEURBANNE CEDEX, FRANCE
allo : +33 472 43 27 56 fax: +33 472 43 13 88
http://pbil.univ-lyon1.fr/members/lobry/
__
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] merge and polylist
On Sat, 7 Oct 2006, Mihai Nica wrote: Greetings: I would like to kindly ask for a little help. The rough code is: Maybe R-sig-geo would be a more relevant list. Comments inline. # dat=data.frame(read.delim(file=all.txt, header = TRUE, sep = \t, quote=\, dec=.,na.strings = NA)) We do not know if dat is in the same order as the shapefile. If dat$cod is malformed wrt. nc$att.data$AREA (very strange choice, in ESRI generated files often the polygon area), you are asking for trouble. nc=read.shape(astae.shp, dbf.data=TRUE, verbose=TRUE) mappolys=Map2poly(nc) submap - subset(mappolys, nc$att.data$NAME!=Honolulu, HI) nc$att.data=subset(nc$att.data, nc$att.data$NAME!=Honolulu, HI) In situations like this, overwriting the input is not advisable, bercause you destroy your ability to check that the output corresponds to your intentions. nc$att.data[,1]=as.numeric(paste(nc$att.data$MSACMSA)) #attributes(nc$att.data) nc$att.data=merge(nc$att.data, dat, by.x=AREA, by.y=cod, all.x=TRUE, sort=FALSE) Ditto. #attributes(nc$att.data) tmp=file(tmp) write.polylistShape(submap, nc$att.data, tmp) Any good reason for not using the sp class framework? The objects you are using here are very low-level and messy. library(maptools) nc_1 - readShapePoly(system.file(shapes/sids.shp, package=maptools)[1], ID=FIPS) # put shapefile in SpatialPolygonsDataFrame nc_2 - nc_1[coordinates(nc_1)[,1] -80,] # subset with [ method row.names(as(nc_2, data.frame)) as.character(nc_2$FIPS) tmpfl - paste(tempfile(), dbf, sep=.) download.file(http://spatial.nhh.no/misc/nc_xtra.dbf;, tmpfl, mode=wb) nc.df - read.dbf(tmpfl) # extra data keyed on CNTY_ID nc_2$CNTY_ID nc.df$CNTY_ID nc_df2 - merge(as(nc_2, data.frame), nc.df, by=CNTY_ID, sort=FALSE) all.equal(nc_df2$CNTY_ID, nc_2$CNTY_ID) row.names(nc_df2) - nc_df2$FIPS # re-instate IDs nc_3 - SpatialPolygonsDataFrame(as(nc_2, SpatialPolygons), data=nc_df2) # (if data row names do not match polygon IDs, will reorder, or fail if # any differ or absent writePolyShape(nc_3, nc_3) This still isn't as tidy as it could be, but gives much more control than the original old-style classes. Roger #_ All works fine, but merge() changes the rownames and the link between the polygons and the corresponding rows is lost. I tried numerous other solutions (such as to paste back the old rownames), to no avail. After a few days, here I am. Please, if you have a moment, send a tip. Thanks, mihai -- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 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.
Re: [R] Probability of exceedance function question
On Sun, 8 Oct 2006, Thomas P. Colson wrote: I'm trying to calculate a cumulative area distribution (graph) of drainage areas. This is defined as P(A A*). Simple in principle. I can do this in excel, with COUNTIF, which will count the number of cells in the row area that have area A, then determine, for each cell in the row area, how many cells exceede that area, then dividing that number by the total number of cells, which gives me the probability that drainage area A exceeds drainage area A*. Is this ecdf() of the vector or its suitable subset? If so, it runs very fast even for large data sets. For plotting, bear in mind that you are generating a lot of output, though: t0 - runif(10) system.time(t1 - ecdf(t0)) [1] 0.222 0.022 0.248 0.000 0.000 system.time(plot(t1, pch=.)) [1] 1.089 0.079 1.186 0.000 0.000 isn't at all bad! E.g, drainage area of 6 sq meters (One DEM grid cell) has a high probability of exceedance(.99), while a drainage area of 100,000 square meters has a low probability of exceedance (.001). I wish to plot this relationship, and we all know that excel is not the tool of choice when working with hundreds of thousands of records. I'd like to port the CAD into a few R functions that I've already developed for other tests as well. So my challenge, in R, is to (1)count the number of rows in column Area that have AREA(*), (2) determine, by row, how many rows have an area greater than the area given in that one row (3) divide step 2 by number of rows (how can I do a row count and port that to a variable, as I have to do this on 10 datasets?) Thanks for any advice you can offer to this endevour __ 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. -- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 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] drawing a curly bracket in a plot
I need to re-draw a graph that consists of a few lines. Next to
each line is a curly bracket - like { - with its ends near the
ends of the line. (At the point of the bracket is some sort of
label, but that is done with text() or something like that.)
Some brackets are horizontal, some are vertical, and some are
diagonal.
I can't find any code for this. Has anyone done it?
Jon
__
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] Probability of exceedance function question
I am able to convert the flow accumulation grid into an area (for each pixel) grid, then import this into R as an ASCII file. The plot(ecdf) function in R seems to plot the opposite: curve starts at probability 0, for drainage area 0, should be the other way? About 150,000 data points in these sets, ecdf curve plots in about 15 seconds. Could the problem be, how I'm importing the data from ascii grid? Cellsize is 20 ft and z is the drainage area, for each cell (flow weighted) area - read.table(file = c:/temp/area.asc, sep = , na.strings = -, skip = 6) area - area[,-ncol(area)] xLLcorner - 1985649.0700408898 yLLcorner - 841301.04004059616 cellsize -20 xURcorner - xLLcorner + (cellsize * (ncol(area) - 1)) xLRcorner - xURcorner xULcorner - xLLcorner yULcorner - yLLcorner + (cellsize * (nrow(area) - 1)) yURcorner - yULcorner yLRcorner - yLLcorner coordsa - expand.grid(y = seq(yULcorner, yLRcorner, by = -20),x = seq(xULcorner, xLRcorner, by = 20)) area- data.frame(coordsa, tmin = as.vector(c(area,recursive = T))) names(area)-c(x,y,z) Plot(ecdf(area$z)) -Original Message- From: Roger Bivand [mailto:[EMAIL PROTECTED] Sent: Sunday, October 08, 2006 4:37 PM To: Thomas P. Colson Cc: r-help@stat.math.ethz.ch Subject: Re: [R] Probability of exceedance function question On Sun, 8 Oct 2006, Thomas P. Colson wrote: I'm trying to calculate a cumulative area distribution (graph) of drainage areas. This is defined as P(A A*). Simple in principle. I can do this in excel, with COUNTIF, which will count the number of cells in the row area that have area A, then determine, for each cell in the row area, how many cells exceede that area, then dividing that number by the total number of cells, which gives me the probability that drainage area A exceeds drainage area A*. Is this ecdf() of the vector or its suitable subset? If so, it runs very fast even for large data sets. For plotting, bear in mind that you are generating a lot of output, though: t0 - runif(10) system.time(t1 - ecdf(t0)) [1] 0.222 0.022 0.248 0.000 0.000 system.time(plot(t1, pch=.)) [1] 1.089 0.079 1.186 0.000 0.000 isn't at all bad! E.g, drainage area of 6 sq meters (One DEM grid cell) has a high probability of exceedance(.99), while a drainage area of 100,000 square meters has a low probability of exceedance (.001). I wish to plot this relationship, and we all know that excel is not the tool of choice when working with hundreds of thousands of records. I'd like to port the CAD into a few R functions that I've already developed for other tests as well. So my challenge, in R, is to (1)count the number of rows in column Area that have AREA(*), (2) determine, by row, how many rows have an area greater than the area given in that one row (3) divide step 2 by number of rows (how can I do a row count and port that to a variable, as I have to do this on 10 datasets?) Thanks for any advice you can offer to this endevour __ 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. -- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 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.
Re: [R] Block comments in R?
I wrote:
R documentation comments really belong
in [.Rd] files where help() can find them.
Barry Rowlingson [EMAIL PROTECTED] replied:
R documentation comments belong in .Rd files at the moment, but how
joyous would it be if they could be included in the .R files?
How joyous? About as joyous as a root canal job without anaesthetic.
I used to be a fan of including thorough documentation in code.
Then JavaDoc hit the world, and many of the languages I use are following
suit with their own EDoc or PlDoc or whateverDoc mess.
I have read far more Java than I ever wanted to, and the more Java I read
the more I *hate* JavaDoc, and the more I am convinced that if you want
to mix documentation and code you need *really* sophisticated tools (like
Web in its various incarnations) or really simple tools (like Haskell's
Bird Tracks, a notation I have adapted to my own use for several other
languages).
.Rd files are semisophisticated; if JavaDoc is a reliable guide, then
shoving that stuff into .R files would be horrible.
Do I need to point out the single biggest difference between JavaDoc
and .Rd? Maybe I do. .Rd files are *USEFUL*. (Because of references,
examples, consistency checking, c, and because they can describe closely
related GROUPS of functions rather than being nailed to specific methods.
The right documentation about the right topics at the right level of
detail.) JavaDoc-style documentation seems to systematically encourage
bulky documentation of low utility.
Okay, this is all part of my incessant whining to make R more like
Python, but I've found managing separate .Rd and .R files a pain. If .R
files could have embedded documentation you'd have one source for code,
documentation, tests etc. I did play about with this in the Splus days,
attaching documentation strings to functions with attributes, but it
was
just kludgy without a proper mechanism.
Let me point out that right now there is NOTHING stopping anyone mixing
.Rd and .R and test cases as they wish. How? Here's how:
$poomat.1
.TH POOMAT 1 Oct 2006 Version 1.0 USER COMMANDS
.SH NAME
poomat - poor man's Tangle
.SH SYNOPSIS
.B poomat
file ...
.SH DESCRIPTION
.B poomat
extracts several interleaved files (such as documentation,
source code, and test cases) from a single file. You can
pack several files together using shell archives, tar files,
or ZIP files, but those are distribution formats, not meant
to be edited as single files.
.B poomat
lets you scatter a file in pieces, interleaved with other
pieces, so that a function, the documentation for the function,
and the test cases for the function can all be in one place.
.SH OPTIONS
None.
.B poomat
concatenates its input files just like
.IR cat (1)
or
.IR awk (1)
and writes to files named in the input.
.SH INPUT LANGUAGE
The input to
.B poomat
is a sequence of chunks. Each chunk is introduced by a line
consisting of a dollar sign in column 1, immediately followed
by a file name. The first chunk for any file name creates a
new file; remaining chunks for the same file are appended to it.
.SH BUGS
Unlike
.IR tangle (1),
.IR ctangle (1),
and other Literate Programming tools, there is no facility for
re-ordering chunks. Nor is there any macro facility.
.PP
It is up to you to make sure that the file names are portable.
Stick to 8.3 file names without any directory affixes and you
should be right.
$INSTALL
Edit the first line of the poomat file to refer to the right
version of awk (nawk, gawk, mawk) for your system, and then
move the poomat file to some directory in your $PATH.
$poomat
#!/usr/ucb/nawk -f
BEGIN { output = /dev/stdout }
/^[$]/ { output = substr($0, 2); next }
{ print output }
Yes, I do mean the whole thing is a three-line AWK script.
Yes, I do mean it is language-independent, and doesn't NEED to be built
into R or anything else.
Yes, I do mean that the source code, documentation, and test cases get
separated as part of the build process. So?
For people interested in doing stuff like that with C or C++, it really
doesn't take a lot of work to make poomat notice a
/[.][chCH][pxPX+]?[pxPX+]?$/ suffix and emit a #line directive.
OK, here it is:
/^[$]/ {
output = substr($0, 2)
if (/[.][chCH][pxPX+]?[pxPX+]?$/)
print #line, FNR+1, \ FILENAME \ output
next
}
So your C or C++ debugger can refer back to the original file.
Maybe that should be in the official version, but R doesn't need it.
Re: [R] drawing a curly bracket in a plot
Jonathan, I would approach this task by drawing a curly bracket in a grob and then rotate it. I haven't made the attempt. 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.
Re: [R] Block comments in R?
I must say I agree with Richard O'Keefe who wrote: I wrote: R documentation comments really belong in [.Rd] files where help() can find them. Barry Rowlingson [EMAIL PROTECTED] replied: R documentation comments belong in .Rd files at the moment, but how joyous would it be if they could be included in the .R files? How joyous? About as joyous as a root canal job without anaesthetic. etc. snip (a) I don't speak Python and I don't believe I ever will; I'm not sure what Python does; I don't want to know. R does what I want to do. I don't want to have to learn Python syntax to do something with which I am very happy doing the R way currently. (b) Modularity is a ***Good Thing***. Code should be code, documentation should be documentation. They should click together but be kept separate. (c) The R documentation system is neat, efficient, well thought out, and well constructed. It is an intrinsic part of the nature of R. (It's one of the --- many --- reasons I like R better than Splus, which I no longer use.) Changing the R documentation system to something (Monty?) Pythonesque would change the nature of R and make it uglier. cheers, Rolf __ 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] Block comments in R?
Rolf Turner wrote: I must say I agree with Richard O'Keefe who wrote: I wrote: R documentation comments really belong in [.Rd] files where help() can find them. Barry Rowlingson [EMAIL PROTECTED] replied: R documentation comments belong in .Rd files at the moment, but how joyous would it be if they could be included in the .R files? How joyous? About as joyous as a root canal job without anaesthetic. etc. snip (a) I don't speak Python and I don't believe I ever will; I'm not sure what Python does; I don't want to know. R does what I want to do. I don't want to have to learn Python syntax to do something with which I am very happy doing the R way currently. To be fair, I don't think anyone proposed that. As I read it, the suggestion was to take Python's idea of inline documentation into R, not necessarily an identical implementation. (b) Modularity is a ***Good Thing***. Code should be code, documentation should be documentation. They should click together but be kept separate. I tend to like comments in code. Rd-style documentation has a different purpose, but not all that different. I tend not to put so many comments into my R code because I spend the time writing .Rd files, but then it's inconvenient to see the docs when I'm working on the code. (c) The R documentation system is neat, efficient, well thought out, and well constructed. It is an intrinsic part of the nature of R. (It's one of the --- many --- reasons I like R better than Splus, which I no longer use.) Changing the R documentation system to something (Monty?) Pythonesque would change the nature of R and make it uglier. I'd rather wait to see an actual proposal before I decided if it was uglier or more beautiful. Current .Rd documentation has some obvious problems: - the parser strips comments out of examples when it runs them - there's no way to put images into the documentation - the keywords aren't much use - there's isn't a definition anywhere of what the format really is, so it's hard to know whether something will work other than by trying it -- and it may break with the next release. - there's no way to link from a help man page to a vignette or other form of documentation. None of these would be addressed by inline help, but other people may have other complaints. I think it was DSC 2001 where R Core decided to replace the .Rd system with something better; I would say it's a reasonable prediction that that won't happen soon, but incremental improvements to .Rd have happened, and we should think about other ones. Duncan Murdoch __ 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] Block comments in R?
Apparently misunderstanding the argument,
BBands [EMAIL PROTECTED] wrote:
Should R be editor specific? Or should it be editor neutral?
R is, and should remain, editor-neutral.
Amongst other things, it should NOT acquire misfeatures
in order to support editors that happen not to support comment-region.
In my view blocks comments are a desirable, editor-neutral approach.
Block comments are indeed editor-neutral,
but they do not solve any problem that R currently has,
and they ARE in practice highly error-prone.
Note that most of the more recent languages have some form of block
comment capability.
http://en.wikipedia.org/wiki/Comment#Summary
That list is not a list of more recent languages; it is a list containing
both old languages and new ones. Some of the entries (such as the one for
Algol 60), are clearly wrong. (In fact the Algol one is doubly wrong,
because the language always was Algol, not ALGOL.) So is the Scheme entry
wrong: there are no block comments in any official version of Scheme and
never have been. The list is missing Burroughs Algol (% end of line
comments), Prolog (%...\n and /*...*/), Pop-2 (!...! or !...\n), and many
others. If you want to be pedantic, the list is technically wrong about
C and C++. We also find COBOL, with some interesting variants, missing.
Texinfo is another one that's wrong. Quite a lot of entries are wrong.
If more recent languages just means languages that are still in use,
it's a pity APL isn't there; the APL lamp character (comments are supposed
to be illuminating, no?) U+235D is cute.
We find
- languages which have copied PL/I (/* ... */)
- languages which have copied BCPL (//)
- languages which have copied Pascal (* ... *)
- languages which have copied sh (#)
- languages which have copied Burroughs Algol (%)
- Fortran (only ever end-of-line comments, first C and now !)
- BASIC (only end-of-line ' and REM comments, and you would certainly
have to call VB.net recent)
- eclectic languages
- languages designed for high reliability, notably Ada and Eiffel (--)
Basically, most language designers have pretty much blindly followed
designs from the 1960s, with the notable exception of Ada and Eiffel.
The Haskell experience is instructive. The Haskell designers thought they
could put in nesting block comments {-...-} and make them work. But it
took several iterations of fiddling with the details, and they *still*
don't work in every case.
__
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] where should R be installed?
I am a new user of Linux (long time user of Windows) whose only training is from books, articles, and just playing with Linux. I have built and installed R from source, and also installed from RPMs (SuSE).I am wondering where in the file system experienced users of R and Linux typically install R. Are there any conventions, or any reasons to choose one location over another? I know that this is, strictly speaking, not an R question but I thought I might ask here for guidance that might help me organize things for ease of use and maintenance, upgrading of packages etc. Any suggestions gratefully accepted. Dan Daniel Nordlund Bothell, WA __ 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] read.zoo question
I have comma delimited asci data with each row being in the format : 2006-01-24 02:41:24.00011,1.2293,5,1.2295,7 . . . . and i'm trying to use read.zoo ( which is similar to read.table ) to read in the data. the data goes all the way out to milliseconds and i can't figure out what to put for the format field. if i put nothing, then read.zoo gets rid of the the whole time and just keeps the date but i want to keep the time. i'm sure gabor knows but he may not be on so maybe someone else knows ? I've done ?format and ?read.zoo and i've read through the zoo article in the journal of statistical software but i can't figure it out. even if someone could point me in the right direction for looking ( i.e a ?whatever ) that would be a great help. thanks. This is not an offer (or solicitation of an offer) to buy/sell the securities/instruments mentioned or an official confirmation. Morgan Stanley may deal as principal in or own or act as market maker for securities/instruments mentioned or may advise the issuers. This is not research and is not from MS Research but it may refer to a research analyst/research report. Unless indicated, these views are the author's and may differ from those of Morgan Stanley research or others in the Firm. We do not represent this is accurate or complete and we may not update this. Past performance is not indicative of future returns. For additional information, research reports and important disclosures, contact me or see https://secure.ms.com/servlet/cls. You should not use e-mail to request, authorize or effect the purchase or sale of any security or instrument, to send transfer instructions, or to effect any other transactions. We cannot guarantee that any such requests received via ! e-mail will be processed in a timely manner. This communication is solely for the addressee(s) and may contain confidential information. We do not waive confidentiality by mistransmission. Contact me if you do not wish to receive these communications. In the UK, this communication is directed in the UK to those persons who are market counterparties or intermediate customers (as defined in the UK Financial Services Authority's rules). [[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] read.zoo question
Try
read.zoo(myfile, sep = ,, FUN = as.POSIXct)
or
to.chron - function(x) {
s - do.call(rbind, strsplit(format(x), ))
chron(dates(s[,1], format = Y-M-D), times(s[,2]))
}
read.zoo(myfile, sep = ,, FUN = to.chron)
depending on which class you want.
On 10/8/06, Leeds, Mark (IED) [EMAIL PROTECTED] wrote:
I have comma delimited asci data with each row being in the format :
2006-01-24 02:41:24.00011,1.2293,5,1.2295,7
.
.
.
.
and i'm trying to use read.zoo ( which is similar to read.table ) to
read in the data. the data goes all the way out to milliseconds and i
can't figure out what to put for the format field. if i put nothing,
then read.zoo gets rid of the the whole time and just keeps the date but
i want to keep the time. i'm sure gabor knows but he may not be on so
maybe someone else knows ? I've done ?format and ?read.zoo and i've read
through the zoo article in the journal of statistical software but i
can't figure it out. even if someone could point me in the right
direction for looking ( i.e a ?whatever ) that would be a great help.
thanks.
This is not an offer (or solicitation of an offer) to buy/sell the
securities/instruments mentioned or an official confirmation. Morgan Stanley
may deal as principal in or own or act as market maker for
securities/instruments mentioned or may advise the issuers. This is not
research and is not from MS Research but it may refer to a research
analyst/research report. Unless indicated, these views are the author's and
may differ from those of Morgan Stanley research or others in the Firm. We
do not represent this is accurate or complete and we may not update this.
Past performance is not indicative of future returns. For additional
information, research reports and important disclosures, contact me or see
https://secure.ms.com/servlet/cls. You should not use e-mail to request,
authorize or effect the purchase or sale of any security or instrument, to
send transfer instructions, or to effect any other transactions. We cannot
guarantee that any such requests received vi!
a !
e-mail will be processed in a timely manner. This communication is solely
for the addressee(s) and may contain confidential information. We do not
waive confidentiality by mistransmission. Contact me if you do not wish to
receive these communications. In the UK, this communication is directed in
the UK to those persons who are market counterparties or intermediate
customers (as defined in the UK Financial Services Authority's rules).
[[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.
Re: [R] where should R be installed?
I usually install R somewhere off my home directory because it saves me the trouble of having to become root and install somewhere in the system. I'm the only person who uses the system anyway. -roger On 10/8/06, Daniel Nordlund [EMAIL PROTECTED] wrote: I am a new user of Linux (long time user of Windows) whose only training is from books, articles, and just playing with Linux. I have built and installed R from source, and also installed from RPMs (SuSE).I am wondering where in the file system experienced users of R and Linux typically install R. Are there any conventions, or any reasons to choose one location over another? I know that this is, strictly speaking, not an R question but I thought I might ask here for guidance that might help me organize things for ease of use and maintenance, upgrading of packages etc. Any suggestions gratefully accepted. Dan Daniel Nordlund Bothell, WA __ 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. -- Roger D. Peng | http://www.biostat.jhsph.edu/~rpeng/ __ 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] lattice: adding text to plots
Hi Everyone,
I know there was a thread recently that dealt with a similar issue,
but this one is a little different.
I have the following data frame
DF - data.frame(x = 1:12, y1 = rnorm(12), y2 = rnorm(12), g =
gl(2,6), h = rep(c(1, 2), 6), f = c(rep(c(1-1,1-2),3),
rep(c(2-1,2-2),3)))
I essentially want this plot
xyplot(y1+y2~x|g*h, data=DF, type=l)
However, now I want to add a different text to each panel, the text
being from column f of the data frame. In other words, I want text
1-1 in the panel where g=1 and h=1 and so on. I tried to pass groups
and subscript to the panel function but could not get what I was
looking for. The following attempt does not work.
xyplot(y1+y2~x|g*h, data=DF, type=l, auto.key=TRUE,
panel=function(x,y,..., groups, subscripts){panel.xyplot(x,y,...);
panel.text(x=4,y=0, labels=DF$f[subscripts])})
My R version is 2.2.1 and lattice version is 0.12-11 (sorry they are
not the latest ones, these are on the server).
--
Ritwik Sinha
Graduate Student
Epidemiology and Biostatistics
Case Western Reserve University
[EMAIL PROTECTED] | +12163682366 | http://darwin.cwru.edu/~rsinha
__
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] lattice: adding text to plots
Could you explain what does not work means. It seems to produce a
graph with x-y numbers on it in R 2.4.0 on Windows.
At any rate, I would have done it like this although I think you can
leave off the [1] on subscripts and it will still work.
library(lattice)
library(grid)
xyplot(y1 + y2 ~ x | g * h, data = DF, type = l,
panel = function(x, y, subscripts, groups, ...) {
panel.xyplot(x, y, ...)
grid.text(DF$f[subscripts[1]], .1, .9)
})
On 10/9/06, Ritwik Sinha [EMAIL PROTECTED] wrote:
Hi Everyone,
I know there was a thread recently that dealt with a similar issue,
but this one is a little different.
I have the following data frame
DF - data.frame(x = 1:12, y1 = rnorm(12), y2 = rnorm(12), g =
gl(2,6), h = rep(c(1, 2), 6), f = c(rep(c(1-1,1-2),3),
rep(c(2-1,2-2),3)))
I essentially want this plot
xyplot(y1+y2~x|g*h, data=DF, type=l)
However, now I want to add a different text to each panel, the text
being from column f of the data frame. In other words, I want text
1-1 in the panel where g=1 and h=1 and so on. I tried to pass groups
and subscript to the panel function but could not get what I was
looking for. The following attempt does not work.
xyplot(y1+y2~x|g*h, data=DF, type=l, auto.key=TRUE,
panel=function(x,y,..., groups, subscripts){panel.xyplot(x,y,...);
panel.text(x=4,y=0, labels=DF$f[subscripts])})
My R version is 2.2.1 and lattice version is 0.12-11 (sorry they are
not the latest ones, these are on the server).
--
Ritwik Sinha
Graduate Student
Epidemiology and Biostatistics
Case Western Reserve University
[EMAIL PROTECTED] | +12163682366 | http://darwin.cwru.edu/~rsinha
__
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.
