[R] Propensity score and three treatments
Dear All, I would like to find something ( references, code,..) to implement a comparison of three treatments in an observational study using the 'Propensity Score'. Any help is much appreciated. Thanks! Giovanni -- dr. Giovanni Parrinello Department of Biotecnologies Medical Statistics Unit University of Brescia Viale Europa, 11 25123 Brescia email: [EMAIL PROTECTED] Phone: +390303717528 Fax: +390303717488 __ 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] Different result from nls in R-2.2.1 and R-2.3.1
The way nls looks for variables in its formula is very unusual. For
non-parameters it has
for(var in varNames[!varIndex])
mf[[var]] - eval(as.name(var), data)
inside nls, and that means that the scoping rules are first 'data' then
the body of nls(). That's surely not what one would expect here:
fitting functions using model.frame look in 'data' and then the
environment of the formula (if non-null) or the parent frame.
So to avoid surprises, put fixed values in 'data'. (Your example did not
work in 2.2.1 when run from a function.)
I'll try to find a bug fix for 2.4.0.
On Thu, 21 Sep 2006, Frede Aakmann Tøgersen wrote:
Short story: January 2006 I did some analysis in R-2.2.1 using nls. Repeating
the exercise in R-2.3.1 yesterday produced somewhat different results.
After some debugging I found that either nls is the problem or that mine
understanding of environments or scoping rules is lacking something.
This is a short reproducing example.
x - seq(0,5,len=20)
n - 1
y - 2*x^2 + n + rnorm(x)
xy - data.frame(x=x,y=y)
myf - function(x,a,b,n){
res - a*x^b + n
## a print for debugging purpose
print(n)
res
}
## This works as I expect it to do in R-2.2.1 but doesn't work in R-2.3.1.
## n is somehow sat to nrow(xy) inside nls()
## Note that x and y is defined in the dataframe xy, whereas n is found in the
global environment.
fit - nls(y ~ myf(x,a,b,n), data=xy, start=c(a=1,b=1), trace=TRUE)
## this works in both versions
## x,y,n found in the .GlobalEnv
fit - nls(y ~ myf(x,a,b,n), start=c(a=1,b=1), trace=TRUE)
## this works in both versions.
## x, y, n found in dataframe xyn
xyn - data.frame(xy,n=n)
fit - nls(y ~ myf(x,a,b,n), data=xyn, start=c(a=1,b=1), trace=TRUE)
## this works in both versions
## Now using the variable .n instead of n
## .n is found in .GlobaEnv
.n - 1
fit - nls(y ~ myf(x,a,b,.n), data=xy, start=c(a=1,b=1), trace=TRUE)
In my real case and the example above, I do have three or more parameters of
which fitting is done only on few of theme. Is this a problem? Or should I ask,
why is this a problem in R-2.3.1 but not in R-2.2.1?
Is my problem related to this difference between lines of code from nls:
R-2.2.1: mf - as.list(eval(mf, parent.frame()))
R-2.3.1: mf - eval.parent(mf)
n - nrow(mf)
mf - as.list(mf)
where n is being defined in the scope of nls in the latest version?
Best regards
Frede Aakmann Tøgersen
Danish Institute of Agricultural Sciences
Research Centre Foulum
Dept. of Genetics and Biotechnology
Blichers Allé 20, P.O. BOX 50
DK-8830 Tjele
Phone: +45 8999 1900
Direct: +45 8999 1878
E-mail: [EMAIL PROTECTED]
Web: http://www.agrsci.org
This email may contain information that is confidential.
Any use or publication of this email without written permission from DIAS is
not allowed.
If you are not the intended recipient, please notify DIAS immediately and
delete this email.
__
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.
--
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] Update to Dillo browser question
Hi I asked about if there is any way of opening URLs from the help browser in the same window of the same dillo browser - here is the answer. Just to reiterate: dillo is for me the perfect browser for the help of R when you use ?... Rainer Original Message Subject: Re: [Dillo-dev] Opening new URL in same instance and -s option Date: Thu, 21 Sep 2006 21:31:57 -0400 From: Jorge Arellano Cid [EMAIL PROTECTED] To: Rainer M Krug [EMAIL PROTECTED] References: [EMAIL PROTECTED] On Thu, Sep 21, 2006 at 02:10:28PM +0200, Rainer M Krug wrote: Hi Hi. I have two questions first: I would like to be able to have only one instance of dillo open and when I try to open a new url (from outside dillo), that this url is opened in the old instance. Is this possible, and if yes, how? Not now. The idea is to implement this with dpip (dillo plugin protocol) but it has not being done yet, becuase of lack of manpower. second: I have seen emails concerning running dillo in server mode (dillo -s server). I tried it, but dillo didn't know -s. Is this feature imp[lemented, and if yes, how can I access it? Not in the official dillo. Thanks a lot for a brilliant lightning fast browser, :-) -- Cheers Jorge.- -- Rainer M. Krug, Dipl. Phys. (Germany), MSc Conservation Biology (UCT) Department of Conservation Ecology and Entomology University of Stellenbosch Matieland 7602 South Africa Tel:+27 - (0)72 808 2975 (w) Fax:+27 - (0)21 808 3304 Cell: +27 - (0)83 9479 042 email: [EMAIL PROTECTED] [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] Matrix from a vector?
Hi
other approach is to use embed
embed(c(rhsel,rhsel),length(rhsel))[-1,]
which is a little bit quicker
rhsel-rnorm(1000)
n - length(rhsel)
system.time({
i - sapply(1:n, function(i) (1:n - i)%%n + 1)
x2-matrix(rhsel[i], n, n)
})
system.time(x1-embed(c(rhsel,rhsel),length(rhsel))[-1,])
all.equal(x1,x2)
HTH
Petr
On 21 Sep 2006 at 14:53, Sundar Dorai-Raj wrote:
Date sent: Thu, 21 Sep 2006 14:53:52 -0500
From: Sundar Dorai-Raj [EMAIL PROTECTED]
Organization: PDF Solutions, Inc.
To: kone [EMAIL PROTECTED]
Copies to: r-help@stat.math.ethz.ch
Subject:Re: [R] Matrix from a vector?
kone said the following on 9/21/2006 2:30 PM:
Hi,
Is there some function, which generates this kind of n x n -matrix
from a vector?
rhset
[1] 1792 256 13312 512 1024 2048 8192 4096
m=matrix(nrow=length(rhset),ncol=length(rhset))
for(i in 1:length(rhset))
+ {
+ m[,i]=rhset
+ rhset=c(rhset[length(rhset)], rhset[2:length(rhset)-1])
+ }
m
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1792 4096 8192 2048 1024 512 13312 256
[2,] 256 1792 4096 8192 2048 1024 512 13312
[3,] 13312 256 1792 4096 8192 2048 1024 512
[4,] 512 13312 256 1792 4096 8192 2048 1024
[5,] 1024 512 13312 256 1792 4096 8192 2048
[6,] 2048 1024 512 13312 256 1792 4096 8192
[7,] 8192 2048 1024 512 13312 256 1792 4096
[8,] 4096 8192 2048 1024 512 13312 256 1792
Atte Tenkanen
University of Turku, Finland
__
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.
Does this work for you?
rhsel - c(1792, 256, 13312, 512, 1024, 2048, 8192, 4096)
n - length(rhsel)
i - sapply(1:n, function(i) (1:n - i)%%n + 1)
matrix(rhsel[i], n, n)
HTH,
--sundar
__
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] Creating Movies with R
Dear All,
I'd like to know if it is possible to create animations with R.
To be specific, I attach a code I am using for my research to plot
some analytical results in 3D using the lattice package. It is not
necessary to go through the code.
Simply, it plots some 3D density profiles at two different times
selected by the user.
I wonder if it is possible to use the data generated for different
times to create something like an .avi file.
Here is the script:
rm(list=ls())
library(lattice)
# I start defining the analytical functions needed to get the density
as a function of time
expect_position - function(t,lam1,lam2,pos_ini,vel_ini)
{1/(lam1-lam2)*(lam1*exp(lam2*t)-lam2*exp(lam1*t))*pos_ini+
1/(lam1-lam2)*(exp(lam1*t)-exp(lam2*t))*vel_ini
}
sigma_pos-function(t,q,lam1,lam2)
{
q/(lam1-lam2)^2*(
(exp(2*lam1*t)-1)/(2*lam1)-2/(lam1+lam2)*(exp(lam1*t+lam2*t)-1) +
(exp(2*lam2*t)-1)/(2*lam2) )
}
rho_x-function(x,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(x-expect_position)^2/sigma_pos)
}
Now the physical parameters
tau-0.1
beta-1/tau
St-tau ### since I am in dimensionless units and tau is already in
units of 1/|alpha|
D=2e-2
q-2*beta^2*D
### Now the grid in space and time
time-5 # time extent
tsteps-501 # time steps
newtime-seq(0,time,len=tsteps)
Now the things specific for the dynamics along x
lam1- -beta/2*(1+sqrt(1+4*St))
lam2- -beta/2*(1-sqrt(1+4*St))
xmin- -0.5
xmax-0.5
x0-0.1
vx0-x0
nx-101 ## grid intervals along x
newx-seq(xmin,xmax,len=nx) # grid along x
# M1 - do.call(g, c(list(x = newx), mypar))
mypar-c(q,lam1,lam2)
sig_xx-do.call(sigma_pos,c(list(t=newtime),mypar))
mypar-c(lam1,lam2,x0,vx0)
exp_x-do.call(expect_position,c(list(t=newtime),mypar))
#rho_x-function(x,expect_position,sigma_pos)
#NB: at t=0, the density blows up, since I have a delta as the initial state!
# At any t0, instead, the result is finite.
#for this reason I now redefine time by getting rid of the istant t=0
to work out
# the density
rho_x_t-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_x[i],sig_xx[i])
myrho_x-do.call(rho_x,c(list(x=newx),mypar))
rho_x_t[ i-1, ]-myrho_x
}
### Now I also define a scaled density
rho_x_t_scaled-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_x[i],sig_xx[i])
myrho_x-do.call(rho_x,c(list(x=newx),mypar))
rho_x_t_scaled[ i-1, ]-myrho_x/max(myrho_x)
}
###Now I deal with the dynamics along y
lam1- -beta/2*(1+sqrt(1-4*St))
lam2- -beta/2*(1-sqrt(1-4*St))
ymin- 0
ymax- 1
y0-ymax
vy0- -y0
mypar-c(q,lam1,lam2)
sig_yy-do.call(sigma_pos,c(list(t=newtime),mypar))
mypar-c(lam1,lam2,y0,vy0)
exp_y-do.call(expect_position,c(list(t=newtime),mypar))
# now I introduce the function giving the density along y: this has to
include the BC of zero
# density at wall
rho_y-function(y,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y-expect_position)^2/sigma_pos)-
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y+expect_position)^2/sigma_pos)
}
newy-seq(ymin,ymax,len=nx) # grid along y with the same # of points
as the one along x
rho_y_t-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_y[i],sig_yy[i])
myrho_y-do.call(rho_y,c(list(y=newy),mypar))
rho_y_t[ i-1, ]-myrho_y
}
rho_y_t_scaled-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_y[i],sig_yy[i])
myrho_y-do.call(rho_y,c(list(y=newy),mypar))
rho_y_t_scaled[ i-1, ]-myrho_y/max(myrho_y)
}
# The following 2 plots are an example of the plots I'd like to use to
make an animation
g - expand.grid(x = newx, y = newy)
instant-100
mydens-rho_x_t[ instant, ]%o%rho_y_t[ instant, ]/(max(rho_x_t[
instant, ]%o%rho_y_t[ instant, ]))
lentot-nx^2
dim(mydens)-c(lentot,1)
g$z-mydens
jpeg(dens-t-3.jpeg)
print(wireframe(z ~ x * y, g, drape = TRUE,shade=TRUE,
scales = list(arrows = FALSE),pretty=FALSE, aspect = c(1,1), colorkey = TRUE
,zoom=0.8, main=expression(Density at t=2), zlab =
list(expression(density),rot = 90),distance=0.0,
perspective=TRUE,#screen = list(z = 150, x = -55,y= 0)
,zlim=range(c(0,1
dev.off()
instant-300
mydens-rho_x_t[ instant, ]%o%rho_y_t[ instant, ]/(max(rho_x_t[
instant, ]%o%rho_y_t[ instant, ]))
lentot-nx^2
dim(mydens)-c(lentot,1)
g$z-mydens
jpeg(dens-t-3.jpeg)
print(wireframe(z ~ x * y, g, drape = TRUE,shade=TRUE,
scales = list(arrows = FALSE),pretty=FALSE, aspect = c(1,1), colorkey = TRUE
,zoom=0.8, main=expression(Density at t=3), zlab =
list(expression(density),rot = 90),distance=0.0,
perspective=TRUE,#screen = list(z = 150, x = -55,y= 0)
,zlim=range(c(0,1
dev.off()
Kind Regards
Lorenzo
__
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] Exponentiate a matrix
Out of curiosity and possibly for later use: is there an R-function that
does matrix logarithms?
Best, Ingmar
On 9/21/06 8:58 PM, Dimitrios Rizopoulos
[EMAIL PROTECTED] wrote:
Quoting Douglas Bates [EMAIL PROTECTED]:
On 9/21/06, Dimitrios Rizopoulos [EMAIL PROTECTED] wrote:
Quoting Duncan Murdoch [EMAIL PROTECTED]:
On 9/21/2006 10:40 AM, Doran, Harold wrote:
Suppose I have a square matrix P
P - matrix(c(.3,.7, .7, .3), ncol=2)
I know that
P * P
Returns the element by element product, whereas
P%*%P
Returns the matrix product.
Now, P^2 also returns the element by element product. But, is there a
slick way to write
P %*% P %*% P
Obviously, P^3 does not return the result I expect.
I don't think there's anything built in, but it's easy to write your own:
I think there was function mtx.exp() in the Malmig package, but it
seems that this package has been withdrawn from CRAN. An old version
appears to exist in:
http://r.meteo.uni.wroc.pl/src/contrib/Descriptions/Malmig.html
Best,
Dimitris
Is that function for matrix powers or for the exponential of a matrix
(which is what I initally thought that Harold wanted)? There is a
function expm in the Matrix package, patterned on the octave function
of the same name, the calculates the matrix exponential for a square
matrix.
this function calculates the n-th power of a matrix, and this is what
I thought Harold wanted, i.e.,
P %*% P %*% P %*% P
should be equal to
mtx.exp(P, 4)
%^% - function(mat, pow) {
stopifnot(length(pow) == 1, all.equal(pow, round(pow)), nrow(mat) ==
ncol(mat))
pow - round(pow)
if (pow 0) {
mat - solve(mat)
pow - abs(pow)
}
result - diag(nrow(mat))
while (pow 0) {
result - result %*% mat
pow - pow - 1
}
result
}
Now P %^% 3 will give you the matrix cube.
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.
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-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.
--
Ingmar Visser
Department of Psychology, University of Amsterdam
Roetersstraat 15, 1018 WB Amsterdam
The Netherlands
http://users.fmg.uva.nl/ivisser/
tel: +31-20-5256735
__
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] Extracting and using Lattice x or y axis
Dear list, I am using Lattice graphics together with Grid. In order to get completely good results, I would like to copy the X and Y axes from the Lattice graphic and draw these in another Grid Viewport. Is there a way to achieve this? I would be gratefull for your help, Albart Coster __ 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] Creating Movies with R
You can make animated gifs using the write.gif function in the caTools package.
On 22/09/06, Lorenzo Isella [EMAIL PROTECTED] wrote:
Dear All,
I'd like to know if it is possible to create animations with R.
To be specific, I attach a code I am using for my research to plot
some analytical results in 3D using the lattice package. It is not
necessary to go through the code.
Simply, it plots some 3D density profiles at two different times
selected by the user.
I wonder if it is possible to use the data generated for different
times to create something like an .avi file.
Here is the script:
rm(list=ls())
library(lattice)
# I start defining the analytical functions needed to get the density
as a function of time
expect_position - function(t,lam1,lam2,pos_ini,vel_ini)
{1/(lam1-lam2)*(lam1*exp(lam2*t)-lam2*exp(lam1*t))*pos_ini+
1/(lam1-lam2)*(exp(lam1*t)-exp(lam2*t))*vel_ini
}
sigma_pos-function(t,q,lam1,lam2)
{
q/(lam1-lam2)^2*(
(exp(2*lam1*t)-1)/(2*lam1)-2/(lam1+lam2)*(exp(lam1*t+lam2*t)-1) +
(exp(2*lam2*t)-1)/(2*lam2) )
}
rho_x-function(x,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(x-expect_position)^2/sigma_pos)
}
Now the physical parameters
tau-0.1
beta-1/tau
St-tau ### since I am in dimensionless units and tau is already in
units of 1/|alpha|
D=2e-2
q-2*beta^2*D
### Now the grid in space and time
time-5 # time extent
tsteps-501 # time steps
newtime-seq(0,time,len=tsteps)
Now the things specific for the dynamics along x
lam1- -beta/2*(1+sqrt(1+4*St))
lam2- -beta/2*(1-sqrt(1+4*St))
xmin- -0.5
xmax-0.5
x0-0.1
vx0-x0
nx-101 ## grid intervals along x
newx-seq(xmin,xmax,len=nx) # grid along x
# M1 - do.call(g, c(list(x = newx), mypar))
mypar-c(q,lam1,lam2)
sig_xx-do.call(sigma_pos,c(list(t=newtime),mypar))
mypar-c(lam1,lam2,x0,vx0)
exp_x-do.call(expect_position,c(list(t=newtime),mypar))
#rho_x-function(x,expect_position,sigma_pos)
#NB: at t=0, the density blows up, since I have a delta as the initial state!
# At any t0, instead, the result is finite.
#for this reason I now redefine time by getting rid of the istant t=0
to work out
# the density
rho_x_t-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_x[i],sig_xx[i])
myrho_x-do.call(rho_x,c(list(x=newx),mypar))
rho_x_t[ i-1, ]-myrho_x
}
### Now I also define a scaled density
rho_x_t_scaled-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_x[i],sig_xx[i])
myrho_x-do.call(rho_x,c(list(x=newx),mypar))
rho_x_t_scaled[ i-1, ]-myrho_x/max(myrho_x)
}
###Now I deal with the dynamics along y
lam1- -beta/2*(1+sqrt(1-4*St))
lam2- -beta/2*(1-sqrt(1-4*St))
ymin- 0
ymax- 1
y0-ymax
vy0- -y0
mypar-c(q,lam1,lam2)
sig_yy-do.call(sigma_pos,c(list(t=newtime),mypar))
mypar-c(lam1,lam2,y0,vy0)
exp_y-do.call(expect_position,c(list(t=newtime),mypar))
# now I introduce the function giving the density along y: this has to
include the BC of zero
# density at wall
rho_y-function(y,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y-expect_position)^2/sigma_pos)-
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y+expect_position)^2/sigma_pos)
}
newy-seq(ymin,ymax,len=nx) # grid along y with the same # of points
as the one along x
rho_y_t-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_y[i],sig_yy[i])
myrho_y-do.call(rho_y,c(list(y=newy),mypar))
rho_y_t[ i-1, ]-myrho_y
}
rho_y_t_scaled-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_y[i],sig_yy[i])
myrho_y-do.call(rho_y,c(list(y=newy),mypar))
rho_y_t_scaled[ i-1, ]-myrho_y/max(myrho_y)
}
# The following 2 plots are an example of the plots I'd like to use to
make an animation
g - expand.grid(x = newx, y = newy)
instant-100
mydens-rho_x_t[ instant, ]%o%rho_y_t[ instant, ]/(max(rho_x_t[
instant, ]%o%rho_y_t[ instant, ]))
lentot-nx^2
dim(mydens)-c(lentot,1)
g$z-mydens
jpeg(dens-t-3.jpeg)
print(wireframe(z ~ x * y, g, drape = TRUE,shade=TRUE,
scales = list(arrows = FALSE),pretty=FALSE, aspect = c(1,1), colorkey = TRUE
,zoom=0.8, main=expression(Density at t=2), zlab =
list(expression(density),rot = 90),distance=0.0,
perspective=TRUE,#screen = list(z = 150, x = -55,y= 0)
,zlim=range(c(0,1
dev.off()
instant-300
mydens-rho_x_t[ instant, ]%o%rho_y_t[ instant, ]/(max(rho_x_t[
instant, ]%o%rho_y_t[ instant, ]))
lentot-nx^2
dim(mydens)-c(lentot,1)
g$z-mydens
jpeg(dens-t-3.jpeg)
print(wireframe(z ~ x * y, g, drape = TRUE,shade=TRUE,
scales = list(arrows = FALSE),pretty=FALSE, aspect = c(1,1), colorkey = TRUE
,zoom=0.8, main=expression(Density at t=3), zlab =
list(expression(density),rot = 90),distance=0.0,
perspective=TRUE,#screen = list(z = 150, x = -55,y= 0)
,zlim=range(c(0,1
dev.off()
Kind Regards
Lorenzo
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting
[R] Compiling a contingency table of counts by case
I have asked a similar question before but this time the problem is somewhat more involved. I have the following data: case;name;x 1;Joe;1 1;Mike;1 1;Zoe;1 2;Joe;1 2;Mike;0 2;Zoe;1 2;John;1 3;Mike;1 3;Zoe;0 3;Karl;0 I would like to count the number of case in which any two name a. both have x=1, b. the first has x=0 - the second has x=1, c. the first has x=1 - the second has x=0, d. both have x=0, The difficulty is that the number of names and their identity changes from case to case. Thanks a lot for you help, Serguei Kaniovski -- ___ Austrian Institute of Economic Research (WIFO) Name: Serguei Kaniovski P.O.Box 91 Tel.: +43-1-7982601-231 Arsenal Objekt 20 Fax: +43-1-7989386 1103 Vienna, Austria Mail: [EMAIL PROTECTED] http://www.wifo.ac.at/Serguei.Kaniovski __ 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] 'help' information not modified when I modify man files
John Tillinghast wrote: I am updating the Bioconductor package, LMGene. Thus I am modifying someone else's package, editing or writing new documentation, etc. When I modify the man files in LMGene and install the library, it doesn't change the 'help' that R gives you. The 'help' you actually get in R is the same as before. I haven't been able to find an explanation in Writing R Extensions. The package I'm working with also contains a 'help' directory, which is not mentioned in WRE. When installing a package, how do I make sure that the 'man' files get used to generate the 'help' files? The files in ./man are used automatically. Are you sure you are using the modified version of package LMGene rather than the old one (perhaps installed into a different library?)? Uwe Ligges [[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] Merge problem
Hello all, I have read as many merge issues as I possibly could tonight and although I presume this is a small error, I have not found the solution to my problem. I'm trying to merge two data sets: dat0 and TransTable. As you can see below, dat0 has 8000 rows, whereas TransTable has 47296 rows. I would expect when I merge the two data sets, with all.x=F, and all.y=F, that the intersection would yield 8000 rows, considering dat0 is a subset of TransTable. However, I get a neat little surprise when I check the dimensions of the resultant data frame - dat0merge, the merged data frame has 8007 rows! How can this be? Where did these extra 7 rows come from? This appears to defy logic! Thank you in advance for your help. I've put my code below for reference. Tova Fuller dim(dat0) [1] 8000 60 dim(TransTable) [1] 47296 9 dat0merge=merge(TransTable,dat0, by.x=Target,by.y=TargetID,all.x=F,all.y=F) dim(dat0merge) [1] 8007 68 __ 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] Compiling a contingency table of counts by case
dat - read.delim(clipboard, sep=;) dat - dat[order(dat$case, dat$name), ] res - apply(combinations(nlevels(dat$name), 2), 1, function(x) with(dat[dat$name %in% levels(dat$name)[x],], table(unlist(sapply(split(x, case), function(y) ifelse(length(y) == 2, paste(y, collapse=), NA)) names(res) - apply(combinations(nlevels(dat$name), 2), 1, function(x) paste(levels(dat$name)[x], collapse=.)) res $Joe.John 11 1 $Joe.Karl character(0) $Joe.Mike 10 11 1 1 $Joe.Zoe 11 2 $John.Karl character(0) $John.Mike 10 1 $John.Zoe 11 1 $Karl.Mike 01 1 $Karl.Zoe 00 1 $Mike.Zoe 01 10 11 1 1 1 --- Jacques VESLOT CNRS UMR 8090 I.B.L (2ème étage) 1 rue du Professeur Calmette B.P. 245 59019 Lille Cedex Tel : 33 (0)3.20.87.10.44 Fax : 33 (0)3.20.87.10.31 http://www-good.ibl.fr --- Serguei Kaniovski a écrit : I have asked a similar question before but this time the problem is somewhat more involved. I have the following data: case;name;x 1;Joe;1 1;Mike;1 1;Zoe;1 2;Joe;1 2;Mike;0 2;Zoe;1 2;John;1 3;Mike;1 3;Zoe;0 3;Karl;0 I would like to count the number of case in which any two name a. both have x=1, b. the first has x=0 - the second has x=1, c. the first has x=1 - the second has x=0, d. both have x=0, The difficulty is that the number of names and their identity changes from case to case. Thanks a lot for you help, Serguei Kaniovski __ 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 problem
Dear Tova There is no reason why the merged data.frame should have exactely 8000 or less rows. The all=FALSE options only says that now new rows are generated for cases that appear only in one of the two data.frames. Have a look at this sample dat1 - data.frame(a = c(1,2,3,4), b = letters[1:4]) dat2 - data.frame(a = c(1,2,3,4,5,6,7,8,1), b = LETTERS[1:9]) merge(dat1,dat2, by = a, all = FALSE) 1 1a A 2 1a I 3 2b B 4 3c C 5 4d D Since 1 appears twice in the large data.frame it is repeated as the help page ?merge says: If there is more than one match, all possible matches contribute one row each. To compare have a look what the option all = TRUE changes merge(dat1,dat2, by = a, all = TRUE) Probably in your large data frame some rows have identical target ids and get repeated. It should be easy to check it with unique() Hope this helps Christoph -- Credit and Surety PML study: visit our web page www.cs-pml.org -- Christoph Buser [EMAIL PROTECTED] Seminar fuer Statistik, LEO C13 ETH Zurich 8092 Zurich SWITZERLAND phone: x-41-44-632-4673 fax: 632-1228 http://stat.ethz.ch/~buser/ -- Tova Fuller writes: Hello all, I have read as many merge issues as I possibly could tonight and although I presume this is a small error, I have not found the solution to my problem. I'm trying to merge two data sets: dat0 and TransTable. As you can see below, dat0 has 8000 rows, whereas TransTable has 47296 rows. I would expect when I merge the two data sets, with all.x=F, and all.y=F, that the intersection would yield 8000 rows, considering dat0 is a subset of TransTable. However, I get a neat little surprise when I check the dimensions of the resultant data frame - dat0merge, the merged data frame has 8007 rows! How can this be? Where did these extra 7 rows come from? This appears to defy logic! Thank you in advance for your help. I've put my code below for reference. Tova Fuller dim(dat0) [1] 8000 60 dim(TransTable) [1] 47296 9 dat0merge=merge(TransTable,dat0, by.x=Target,by.y=TargetID,all.x=F,all.y=F) dim(dat0merge) [1] 8007 68 __ 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] Merge problem
On Fri, 22 Sep 2006, Tova Fuller wrote: Hello all, I have read as many merge issues as I possibly could tonight and although I presume this is a small error, I have not found the solution to my problem. I'm trying to merge two data sets: dat0 and TransTable. As you can see below, dat0 has 8000 rows, whereas TransTable has 47296 rows. I would expect when I merge the two data sets, with all.x=F, and all.y=F, that the intersection would yield 8000 rows, considering dat0 is a subset of TransTable. However, I get a neat little surprise when I check the dimensions of the resultant data frame - dat0merge, the merged data frame has 8007 rows! How can this be? Where did these extra 7 rows come from? This appears to defy logic! Not the help page, though. joined together. If there is more than one match, all possible matches contribute one row each. I presume you think you are doing a 1:1 match, but probably you have multiple matches for up to 7 ids. This is one of the commonest misconceptions about merge that you will find frequently in the list archives. Thank you in advance for your help. I've put my code below for reference. Tova Fuller dim(dat0) [1] 8000 60 dim(TransTable) [1] 47296 9 dat0merge=merge(TransTable,dat0, by.x=Target,by.y=TargetID,all.x=F,all.y=F) dim(dat0merge) [1] 8007 68 __ 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. -- 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] how to store recursive results
One suggestion: as a recursive list. For example have a look at the 'party' package and the BinaryTree class. I hope this helps. Bálint On 22/09/06, X.H Chen [EMAIL PROTECTED] wrote: Hi all, How to store recursive resutls from a function for each step without using global operators -? Thanks ahead. Xiaohui Chen Dept. of Statistics UBC, Canada _ Don't waste time standing in line—try shopping online. Visit Sympatico / MSN __ 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] how to store recursive results
Note that - is not necessarily global:
if (exists(x)) rm(x)
f - function() {
x - 2
g - function() x - 3
g()
x
}
f() # 3
exists(x) # FALSE
On 9/22/06, X.H Chen [EMAIL PROTECTED] wrote:
Hi all,
How to store recursive resutls from a function for each step without using
global operators -? Thanks ahead.
Xiaohui Chen
Dept. of Statistics
UBC, Canada
_
Don't waste time standing in line—try shopping online. Visit Sympatico / MSN
__
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] Creating Movies with R
On 9/22/2006 3:23 AM, Lorenzo Isella wrote:
Dear All,
I'd like to know if it is possible to create animations with R.
To be specific, I attach a code I am using for my research to plot
some analytical results in 3D using the lattice package. It is not
necessary to go through the code.
Simply, it plots some 3D density profiles at two different times
selected by the user.
I wonder if it is possible to use the data generated for different
times to create something like an .avi file.
There's an example of an animation in the rgl package (see
?rgl.snapshot); it needs ImageMagick to put the frames together. It
supports GIF; I think it also supports MPEG and maybe some other
animated formats.
Duncan Murdoch
Here is the script:
rm(list=ls())
library(lattice)
# I start defining the analytical functions needed to get the density
as a function of time
expect_position - function(t,lam1,lam2,pos_ini,vel_ini)
{1/(lam1-lam2)*(lam1*exp(lam2*t)-lam2*exp(lam1*t))*pos_ini+
1/(lam1-lam2)*(exp(lam1*t)-exp(lam2*t))*vel_ini
}
sigma_pos-function(t,q,lam1,lam2)
{
q/(lam1-lam2)^2*(
(exp(2*lam1*t)-1)/(2*lam1)-2/(lam1+lam2)*(exp(lam1*t+lam2*t)-1) +
(exp(2*lam2*t)-1)/(2*lam2) )
}
rho_x-function(x,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(x-expect_position)^2/sigma_pos)
}
Now the physical parameters
tau-0.1
beta-1/tau
St-tau ### since I am in dimensionless units and tau is already in
units of 1/|alpha|
D=2e-2
q-2*beta^2*D
### Now the grid in space and time
time-5 # time extent
tsteps-501 # time steps
newtime-seq(0,time,len=tsteps)
Now the things specific for the dynamics along x
lam1- -beta/2*(1+sqrt(1+4*St))
lam2- -beta/2*(1-sqrt(1+4*St))
xmin- -0.5
xmax-0.5
x0-0.1
vx0-x0
nx-101 ## grid intervals along x
newx-seq(xmin,xmax,len=nx) # grid along x
# M1 - do.call(g, c(list(x = newx), mypar))
mypar-c(q,lam1,lam2)
sig_xx-do.call(sigma_pos,c(list(t=newtime),mypar))
mypar-c(lam1,lam2,x0,vx0)
exp_x-do.call(expect_position,c(list(t=newtime),mypar))
#rho_x-function(x,expect_position,sigma_pos)
#NB: at t=0, the density blows up, since I have a delta as the initial state!
# At any t0, instead, the result is finite.
#for this reason I now redefine time by getting rid of the istant t=0
to work out
# the density
rho_x_t-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_x[i],sig_xx[i])
myrho_x-do.call(rho_x,c(list(x=newx),mypar))
rho_x_t[ i-1, ]-myrho_x
}
### Now I also define a scaled density
rho_x_t_scaled-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_x[i],sig_xx[i])
myrho_x-do.call(rho_x,c(list(x=newx),mypar))
rho_x_t_scaled[ i-1, ]-myrho_x/max(myrho_x)
}
###Now I deal with the dynamics along y
lam1- -beta/2*(1+sqrt(1-4*St))
lam2- -beta/2*(1-sqrt(1-4*St))
ymin- 0
ymax- 1
y0-ymax
vy0- -y0
mypar-c(q,lam1,lam2)
sig_yy-do.call(sigma_pos,c(list(t=newtime),mypar))
mypar-c(lam1,lam2,y0,vy0)
exp_y-do.call(expect_position,c(list(t=newtime),mypar))
# now I introduce the function giving the density along y: this has to
include the BC of zero
# density at wall
rho_y-function(y,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y-expect_position)^2/sigma_pos)-
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y+expect_position)^2/sigma_pos)
}
newy-seq(ymin,ymax,len=nx) # grid along y with the same # of points
as the one along x
rho_y_t-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_y[i],sig_yy[i])
myrho_y-do.call(rho_y,c(list(y=newy),mypar))
rho_y_t[ i-1, ]-myrho_y
}
rho_y_t_scaled-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_y[i],sig_yy[i])
myrho_y-do.call(rho_y,c(list(y=newy),mypar))
rho_y_t_scaled[ i-1, ]-myrho_y/max(myrho_y)
}
# The following 2 plots are an example of the plots I'd like to use to
make an animation
g - expand.grid(x = newx, y = newy)
instant-100
mydens-rho_x_t[ instant, ]%o%rho_y_t[ instant, ]/(max(rho_x_t[
instant, ]%o%rho_y_t[ instant, ]))
lentot-nx^2
dim(mydens)-c(lentot,1)
g$z-mydens
jpeg(dens-t-3.jpeg)
print(wireframe(z ~ x * y, g, drape = TRUE,shade=TRUE,
scales = list(arrows = FALSE),pretty=FALSE, aspect = c(1,1), colorkey = TRUE
,zoom=0.8, main=expression(Density at t=2), zlab =
list(expression(density),rot = 90),distance=0.0,
perspective=TRUE,#screen = list(z = 150, x = -55,y= 0)
,zlim=range(c(0,1
dev.off()
instant-300
mydens-rho_x_t[ instant, ]%o%rho_y_t[ instant, ]/(max(rho_x_t[
instant, ]%o%rho_y_t[ instant, ]))
lentot-nx^2
dim(mydens)-c(lentot,1)
g$z-mydens
jpeg(dens-t-3.jpeg)
print(wireframe(z ~ x * y, g, drape = TRUE,shade=TRUE,
scales = list(arrows = FALSE),pretty=FALSE, aspect = c(1,1), colorkey = TRUE
,zoom=0.8, main=expression(Density at t=3), zlab =
list(expression(density),rot = 90),distance=0.0,
perspective=TRUE,#screen = list(z = 150, x = -55,y= 0)
,zlim=range(c(0,1
dev.off()
Re: [R] Creating Movies with R
There is also a write.gif function in package caTools.
See the example there.
Best,
Matthias
Dear All,
I'd like to know if it is possible to create animations with R.
To be specific, I attach a code I am using for my research to
plot some analytical results in 3D using the lattice package.
It is not necessary to go through the code.
Simply, it plots some 3D density profiles at two different
times selected by the user.
I wonder if it is possible to use the data generated for
different times to create something like an .avi file.
Here is the script:
rm(list=ls())
library(lattice)
# I start defining the analytical functions needed to get the
density as a function of time
expect_position - function(t,lam1,lam2,pos_ini,vel_ini)
{1/(lam1-lam2)*(lam1*exp(lam2*t)-lam2*exp(lam1*t))*pos_ini+
1/(lam1-lam2)*(exp(lam1*t)-exp(lam2*t))*vel_ini
}
sigma_pos-function(t,q,lam1,lam2)
{
q/(lam1-lam2)^2*(
(exp(2*lam1*t)-1)/(2*lam1)-2/(lam1+lam2)*(exp(lam1*t+lam2*t)-1) +
(exp(2*lam2*t)-1)/(2*lam2) )
}
rho_x-function(x,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(x-expect_position)^2/sigma_pos)
}
Now the physical parameters
tau-0.1
beta-1/tau
St-tau ### since I am in dimensionless units and tau is
already in units of 1/|alpha|
D=2e-2
q-2*beta^2*D
### Now the grid in space and time
time-5 # time extent
tsteps-501 # time steps
newtime-seq(0,time,len=tsteps)
Now the things specific for the dynamics along x
lam1- -beta/2*(1+sqrt(1+4*St))
lam2- -beta/2*(1-sqrt(1+4*St))
xmin- -0.5
xmax-0.5
x0-0.1
vx0-x0
nx-101 ## grid intervals along x
newx-seq(xmin,xmax,len=nx) # grid along x
# M1 - do.call(g, c(list(x = newx), mypar))
mypar-c(q,lam1,lam2)
sig_xx-do.call(sigma_pos,c(list(t=newtime),mypar))
mypar-c(lam1,lam2,x0,vx0)
exp_x-do.call(expect_position,c(list(t=newtime),mypar))
#rho_x-function(x,expect_position,sigma_pos)
#NB: at t=0, the density blows up, since I have a delta as
the initial state!
# At any t0, instead, the result is finite.
#for this reason I now redefine time by getting rid of the
istant t=0 to work out # the density
rho_x_t-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_x[i],sig_xx[i])
myrho_x-do.call(rho_x,c(list(x=newx),mypar))
rho_x_t[ i-1, ]-myrho_x
}
### Now I also define a scaled density
rho_x_t_scaled-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_x[i],sig_xx[i])
myrho_x-do.call(rho_x,c(list(x=newx),mypar))
rho_x_t_scaled[ i-1, ]-myrho_x/max(myrho_x) }
###Now I deal with the dynamics along y
lam1- -beta/2*(1+sqrt(1-4*St))
lam2- -beta/2*(1-sqrt(1-4*St))
ymin- 0
ymax- 1
y0-ymax
vy0- -y0
mypar-c(q,lam1,lam2)
sig_yy-do.call(sigma_pos,c(list(t=newtime),mypar))
mypar-c(lam1,lam2,y0,vy0)
exp_y-do.call(expect_position,c(list(t=newtime),mypar))
# now I introduce the function giving the density along y:
this has to include the BC of zero # density at wall
rho_y-function(y,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y-expect_position)^2/sigma_pos)-
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y+expect_position)^2/sigma_pos)
}
newy-seq(ymin,ymax,len=nx) # grid along y with the same # of
points as the one along x
rho_y_t-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_y[i],sig_yy[i])
myrho_y-do.call(rho_y,c(list(y=newy),mypar))
rho_y_t[ i-1, ]-myrho_y
}
rho_y_t_scaled-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_y[i],sig_yy[i])
myrho_y-do.call(rho_y,c(list(y=newy),mypar))
rho_y_t_scaled[ i-1, ]-myrho_y/max(myrho_y) }
# The following 2 plots are an example of the plots I'd like
to use to make an animation
g - expand.grid(x = newx, y = newy)
instant-100
mydens-rho_x_t[ instant, ]%o%rho_y_t[ instant,
]/(max(rho_x_t[ instant, ]%o%rho_y_t[ instant, ]))
lentot-nx^2
dim(mydens)-c(lentot,1)
g$z-mydens
jpeg(dens-t-3.jpeg)
print(wireframe(z ~ x * y, g, drape = TRUE,shade=TRUE, scales
= list(arrows = FALSE),pretty=FALSE, aspect = c(1,1),
colorkey = TRUE ,zoom=0.8, main=expression(Density at t=2),
zlab = list(expression(density),rot = 90),distance=0.0,
perspective=TRUE,#screen = list(z = 150, x = -55,y= 0)
,zlim=range(c(0,1
dev.off()
instant-300
mydens-rho_x_t[ instant, ]%o%rho_y_t[ instant,
]/(max(rho_x_t[ instant, ]%o%rho_y_t[ instant, ]))
lentot-nx^2
dim(mydens)-c(lentot,1)
g$z-mydens
jpeg(dens-t-3.jpeg)
print(wireframe(z ~ x * y, g, drape = TRUE,shade=TRUE, scales
= list(arrows = FALSE),pretty=FALSE, aspect = c(1,1),
colorkey = TRUE ,zoom=0.8, main=expression(Density at t=3),
zlab = list(expression(density),rot = 90),distance=0.0,
perspective=TRUE,#screen = list(z = 150, x = -55,y= 0)
,zlim=range(c(0,1
dev.off()
Kind Regards
Lorenzo
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
[R] A simple resampling problem
Dear UseRs I would like to show my students how to use resampling to solve the following simple problem: If a family has two children of which one is a boy, what is the probability that the other child is also a boy. The answer is (obviously) 1/3, and can be show easily using the usual methods. But I would like to get the students to think of resampling, by doing the following: Flip two coins repeatedly, denoted 0 and 1 (1 for boy, say). Discard those pairs that both contain 0. From those left over, count how many pairs are (1,1). Divide this number by the number available to choose from (i.e. all pairs, except (0,0)). This will then give 1/3 (more or less of course). Can somebody help me to code this efficiently, or elegantly (and smartly) in R, without loops etc. It is intended for first year students that are only starting to learn about statistics (or probability theory), and R of course. Thank you for your time. Regards Jacob Jacob L van Wyk Department of Statistics University of Johannesburg, APK P O Box 524 Auckland Park 2006 South Africa Tel: +27 11 489 3080 Fax: +27 11 489 2832 [[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] delete a entire vector of a dataframe
This works too: t.d$V712 - NULL On Thursday 21 September 2006 22:28, Gavin Simpson wrote: On Thu, 2006-09-21 at 20:01 +0200, Thomas Preuth wrote: delete a entire vector of a dataframe Hello, i want to delete a vector and tried rm (t.d$V712). This did not work, message was, could not find variable. I thought the $ defines the vectro in a dataframe, when I just type t.d$V712 the content of this vector is displayed. Greetings, Thomas You can't do that, and that is not what the error message said exactly - which should have told you something was wrong with your thinking as it also said 1: remove: variable $ was not found. Instead, copy over the object, minus the column you want to delete: dat - as.data.frame(matrix(rnorm(100), nrow = 10)) names(dat) - paste(Var, 1:10, sep = _) dat # now we don't want column Var_6 dat - dat[, -6] # or if we don't know which column is Var_6 you could do not.want - which(names(dat) %in% Var_7) # now don't want Var_7 dat - dat[, -not.want] dat This can be extended to many variables: not.want - which(names(dat) %in% c(Var_10, Var_2, Var_8)) dat - dat[, -not.want] dat # only 1, 3, 4, 5, 9 left HTH G -- Adrian DUSA Arhiva Romana de Date Sociale Bd. Schitu Magureanu nr.1 050025 Bucuresti sectorul 5 Romania Tel./Fax: +40 21 3126618 \ +40 21 3120210 / int.101 __ 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] A simple resampling problem
try the following: B - 1 index - rowSums(matrix(rbinom(B*2, 1, 0.5), ncol = 2)) sum(index == 2) / sum(index 0) 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 - Original Message - From: Jacob van Wyk [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Sent: Friday, September 22, 2006 12:57 PM Subject: [R] A simple resampling problem Dear UseRs I would like to show my students how to use resampling to solve the following simple problem: If a family has two children of which one is a boy, what is the probability that the other child is also a boy. The answer is (obviously) 1/3, and can be show easily using the usual methods. But I would like to get the students to think of resampling, by doing the following: Flip two coins repeatedly, denoted 0 and 1 (1 for boy, say). Discard those pairs that both contain 0. From those left over, count how many pairs are (1,1). Divide this number by the number available to choose from (i.e. all pairs, except (0,0)). This will then give 1/3 (more or less of course). Can somebody help me to code this efficiently, or elegantly (and smartly) in R, without loops etc. It is intended for first year students that are only starting to learn about statistics (or probability theory), and R of course. Thank you for your time. Regards Jacob Jacob L van Wyk Department of Statistics University of Johannesburg, APK P O Box 524 Auckland Park 2006 South Africa Tel: +27 11 489 3080 Fax: +27 11 489 2832 [[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. 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] A simple resampling problem
Try this: set.seed(1) n - 100 first - sample(0:1, n, replace = TRUE) second - sample(0:1, n, replace = TRUE) sum(first == 1 second == 1) / sum(first == 1 | second == 1) # The last could be written more compactly like this # as 1 and 0 will be treated as TRUE and FALSE by the # logical operators and then TRUE and FALSE will be # regarded as 1 and 0 by sum sum(first second) / sum(first | second) On 9/22/06, Jacob van Wyk [EMAIL PROTECTED] wrote: Dear UseRs I would like to show my students how to use resampling to solve the following simple problem: If a family has two children of which one is a boy, what is the probability that the other child is also a boy. The answer is (obviously) 1/3, and can be show easily using the usual methods. But I would like to get the students to think of resampling, by doing the following: Flip two coins repeatedly, denoted 0 and 1 (1 for boy, say). Discard those pairs that both contain 0. From those left over, count how many pairs are (1,1). Divide this number by the number available to choose from (i.e. all pairs, except (0,0)). This will then give 1/3 (more or less of course). Can somebody help me to code this efficiently, or elegantly (and smartly) in R, without loops etc. It is intended for first year students that are only starting to learn about statistics (or probability theory), and R of course. Thank you for your time. Regards Jacob Jacob L van Wyk Department of Statistics University of Johannesburg, APK P O Box 524 Auckland Park 2006 South Africa Tel: +27 11 489 3080 Fax: +27 11 489 2832 [[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] Propensity score and three treatments
Giovanni Parrinello wrote:
Dear All,
I would like to find something ( references, code,..) to implement a
comparison of three
treatments in an observational study using the 'Propensity Score'.
Any help is much appreciated. Thanks!
Giovanni
@Article{tch05use,
author = {Tchernis, Rusty and {Horvitz-Lennon}, Marcela and
Normand, {Sharon-Lise} T.},
title ={On the use of discrete choice models for
causal infere
nce},
journal = Stat in Med,
year = 2005,
volume = 24,
pages ={2197-2212},
annote = {causal inference;discrete choice models;matching
estimator;multi-valued treatment propensity model using discrete
choice model}
}
In one application I created 2 propensity scores and adjusted for those
using spline functions of the logits. See
@ARTICLE{mar94con,
author = {Mark, D. B. and Nelson, C. L. and Califf, R. M. and
Harrell, F. E.
and Lee, K. L. and Jones, R. H. and Fortin, D. F. and Stack,
R. S.
and Glower, D. D. and Smith, L. R. and {DeLong}, E. R. and
Smith,
P. K. and Reves, J. G. and Jollis, J. G. and Tcheng, J. E. and
Muhlbaier, L. H. and Lowe, J. E. and Phillips, H. R. and
Pryor, D.
B.},
year = 1994,
title = {The continuing evolution of therapy for coronary artery disease:
{Initial} results from the era of coronary angioplasty},
journal = {Circulation},
volume = 89,
pages = {2015-2025},
annote = {CABG; PTCA; propensity; observational study; Cox model
applications; adjusted survival curves;two propensity
scores for three treatments}
}
I have many annotated rreferences about propensity scores in my
bibliographic database - see the bottom of biostat.mc.vanderbilt.edu/rms
--
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.
Re: [R] Compiling a contingency table of counts by case
Thanks Jacques, this works! Serguei -- ___ Austrian Institute of Economic Research (WIFO) Name: Serguei Kaniovski P.O.Box 91 Tel.: +43-1-7982601-231 Arsenal Objekt 20 Fax: +43-1-7989386 1103 Vienna, Austria Mail: [EMAIL PROTECTED] http://www.wifo.ac.at/Serguei.Kaniovski __ 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] Variable as color in a barplot
Dear wise ones, I have a problem assigning different colors to bars in a barplot. The data I'm using is the following dataframe (truncated) : L0 r n p t [...] 18 19 1 1 RFM 19 20 1 1 RFM 20 21 2 1 RFM 21 23 6 1 RIH 22 24 2 1 ROC 23 25 1 1 ROC 24 26 1 1 ROC 25 27 2 1 ROC 26 28 2 1 RFT 27 29 1 1 RFT 28 30 2 1 RFT 29 31 1 1 ROH [...] My barplot should display ascending bars according to L0$n. Their width depends on L0$p, which can be greater than 1. Last, I want a different color for each of the five values of L0$t. I use the following array to assign a colour for each type given in L0$t : coul - array() coul['ROC'] - 'khaki1' coul['RFM'] - 'palegreen' coul['RFT'] - 'lightorange' coul['RIH'] - 'blue' coul['ROH'] - 'lightblue' And here's the barplot command I'm using : barplot(sort(L0$n), ylim=c(0,10), width=L0$p[order(L0$n)], col=coul[L0$t[order(L0$n)]], main=Nombre de genres, sub=Livrets) The barplot uses the sort() and order() functions to display ascending bars with the right values. The width parameter works fine, but the color doesn't. I get : Erreur dans rect(y1, x1, y2, x2, ...) : nom de couleur incorrect (error in rect(y1, x1, y2, x2, ...) : incorrect color name) I tried to fiddle with coul[] or its indice but it didn't help. Does anybody know what's wrong ? Many thanks, Octave __ 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] Variable as color in a barplot
I think it is because your command: coul - array() means that the first element of coul is NA, even after you've added the colour names. On 22/09/06, Octave Julien [EMAIL PROTECTED] wrote: Dear wise ones, I have a problem assigning different colors to bars in a barplot. The data I'm using is the following dataframe (truncated) : L0 r n p t [...] 18 19 1 1 RFM 19 20 1 1 RFM 20 21 2 1 RFM 21 23 6 1 RIH 22 24 2 1 ROC 23 25 1 1 ROC 24 26 1 1 ROC 25 27 2 1 ROC 26 28 2 1 RFT 27 29 1 1 RFT 28 30 2 1 RFT 29 31 1 1 ROH [...] My barplot should display ascending bars according to L0$n. Their width depends on L0$p, which can be greater than 1. Last, I want a different color for each of the five values of L0$t. I use the following array to assign a colour for each type given in L0$t : coul - array() coul['ROC'] - 'khaki1' coul['RFM'] - 'palegreen' coul['RFT'] - 'lightorange' coul['RIH'] - 'blue' coul['ROH'] - 'lightblue' And here's the barplot command I'm using : barplot(sort(L0$n), ylim=c(0,10), width=L0$p[order(L0$n)], col=coul[L0$t[order(L0$n)]], main=Nombre de genres, sub=Livrets) The barplot uses the sort() and order() functions to display ascending bars with the right values. The width parameter works fine, but the color doesn't. I get : Erreur dans rect(y1, x1, y2, x2, ...) : nom de couleur incorrect (error in rect(y1, x1, y2, x2, ...) : incorrect color name) I tried to fiddle with coul[] or its indice but it didn't help. Does anybody know what's wrong ? Many thanks, Octave __ 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. -- = David Barron Said Business School University of Oxford Park End Street Oxford OX1 1HP [[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] 'help' information not modified when I modify man files
The problem was that I was using the unzipped distributed package, rather than the source. Once I got the source I was fine. See http://bioconductor.org/packages/1.8/bioc/html/LMGene.html for links to source and package. Installing from source works differently from installing from a package, and it was messing me up. Thanks John On 9/22/06, Uwe Ligges [EMAIL PROTECTED] wrote: John Tillinghast wrote: I am updating the Bioconductor package, LMGene. Thus I am modifying someone else's package, editing or writing new documentation, etc. When I modify the man files in LMGene and install the library, it doesn't change the 'help' that R gives you. The 'help' you actually get in R is the same as before. I haven't been able to find an explanation in Writing R Extensions. The package I'm working with also contains a 'help' directory, which is not mentioned in WRE. When installing a package, how do I make sure that the 'man' files get used to generate the 'help' files? The files in ./man are used automatically. Are you sure you are using the modified version of package LMGene rather than the old one (perhaps installed into a different library?)? Uwe Ligges [[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. [[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] extract data from lm object and then use again?
Hi list, I want to write a general function so that it would take an lm object, extract its data element, then use the data at another R function (eg, glm). I searched R-help list, and found this would do the trick of the first part: a.lm$call$data this would return a name object but could not be recognized as a data.frameby glm. I also tried call(as.character(a.lm$call$data)) or eval(call(as.character(a.lm$call$data))) neither works. By eval(call(...)), it acts as evaluating of a function, but what I want is just a data frame object which could be inserted into glm function. Anyone could help? Thanks, Mike [[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] extract data from lm object and then use again?
Hi, the data of the model fit is stored in lm$model and should work Good luck, Thilo On Friday 22 September 2006 16:45, Mike Wolfgang wrote: Hi list, I want to write a general function so that it would take an lm object, extract its data element, then use the data at another R function (eg, glm). I searched R-help list, and found this would do the trick of the first part: a.lm$call$data this would return a name object but could not be recognized as a data.frameby glm. I also tried call(as.character(a.lm$call$data)) or eval(call(as.character(a.lm$call$data))) neither works. By eval(call(...)), it acts as evaluating of a function, but what I want is just a data frame object which could be inserted into glm function. Anyone could help? Thanks, Mike [[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. -- Thilo Kellermann Department of Psychiatry und Psychotherapy RWTH Aachen University Pauwelstr. 30 52074 Aachen Tel.: +49 (0)241 / 8089977 Fax.: +49 (0)241 / 8082401 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] how to store recursive results
Hi Patrick,
Thanks for your suggestion. I find your method works for the functions with
integer paramters. For example,
If we have function f:
f-function(i)
{
if(i1)
i*f(i-1)
else
1
}
and then using:
ans - vector(list, n)
for(i in 1:5) {
ans[[i]] - f(i)
}
the ans should be:
[[1]]
[1] 1
[[2]]
[1] 2
[[3]]
[1] 6
[[4]]
[1] 24
[[5]]
[1] 120
But actually, we there is no such a i can be referrenced in f(), no
parametric function for example, this will be a problem. Anyway, thanks a
lot for your suggestions.
Cheers,
Xiaohui Chen
Dept. of Statistics
UBC, Canada
From: Patrick Burns [EMAIL PROTECTED]
To: X.H Chen [EMAIL PROTECTED]
Subject: Re: [R] how to store recursive results
Date: Fri, 22 Sep 2006 10:26:35 +0100
It isn't clear to me exactly what you are asking, but
I think that a list might be what you are after. Something
like:
ans - vector(list, n)
for(i in 1:n) {
ans[[i]] -
}
X.H Chen wrote:
Hi all,
How to store recursive resutls from a function for each step without using
global operators -? Thanks ahead.
Xiaohui Chen
Dept. of Statistics
UBC, Canada
_
Dont waste time standing in linetry shopping online. Visit Sympatico /
MSN
__
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] extract data from lm object and then use again?
On Fri, 2006-09-22 at 10:45 -0400, Mike Wolfgang wrote: Hi list, I want to write a general function so that it would take an lm object, extract its data element, then use the data at another R function (eg, glm). I searched R-help list, and found this would do the trick of the first part: a.lm$call$data this would return a name object but could not be recognized as a data.frameby glm. I also tried call(as.character(a.lm$call$data)) or eval(call(as.character(a.lm$call$data))) neither works. By eval(call(...)), it acts as evaluating of a function, but what I want is just a data frame object which could be inserted into glm function. Anyone could help? Thanks, Mike If the 'data' argument in lm() is used, then this approach could work: Iris2 - eval(lm(Sepal.Length ~ Species, data = iris)$call$data) str(Iris2) `data.frame': 150 obs. of 5 variables: $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ... $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ... $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ... $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ... $ Species : Factor w/ 3 levels setosa,versicolor,..: 1 1 1 1 1 1 1 1 1 1 ... However, note that you do not get the actual data used within the lm() function (the model frame) but the entire source data frame. What you likely want instead is the model frame containing the columns actually used in the model formula: Iris3 - lm(Sepal.Length ~ Species, data = iris)$model str(Iris3) `data.frame': 150 obs. of 2 variables: $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ... $ Species : Factor w/ 3 levels setosa,versicolor,..: 1 1 1 1 1 1 1 1 1 1 ... - attr(*, terms)=Classes 'terms', 'formula' length 3 Sepal.Length ~ Species .. ..- attr(*, variables)= language list(Sepal.Length, Species) .. ..- attr(*, factors)= int [1:2, 1] 0 1 .. .. ..- attr(*, dimnames)=List of 2 .. .. .. ..$ : chr [1:2] Sepal.Length Species .. .. .. ..$ : chr Species .. ..- attr(*, term.labels)= chr Species .. ..- attr(*, order)= int 1 .. ..- attr(*, intercept)= int 1 .. ..- attr(*, response)= int 1 .. ..- attr(*, .Environment)=length 15 environment .. ..- attr(*, predvars)= language list(Sepal.Length, Species) .. ..- attr(*, dataClasses)= Named chr [1:2] numeric factor .. .. ..- attr(*, names)= chr [1:2] Sepal.Length Species 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.
Re: [R] how to store recursive results
Hi Gabor,
Thanks for pointing out this for me. However, what I try to get is how to
construct such form a function f that
ret-f(...),
where ret contains the each recursive result from f, and meantime f consists
of no - operator. Do you have any idea how to implemet this. Thanks a lot
for your suggestions.
Cheer
Xiaohui Chen
Dept. of Statistics
UBC, Canada
From: Gabor Grothendieck [EMAIL PROTECTED]
To: X.H Chen [EMAIL PROTECTED]
CC: r-help@stat.math.ethz.ch
Subject: Re: [R] how to store recursive results
Date: Fri, 22 Sep 2006 06:49:22 -0400
Note that - is not necessarily global:
if (exists(x)) rm(x)
f - function() {
x - 2
g - function() x - 3
g()
x
}
f() # 3
exists(x) # FALSE
On 9/22/06, X.H Chen [EMAIL PROTECTED] wrote:
Hi all,
How to store recursive resutls from a function for each step without using
global operators -? Thanks ahead.
Xiaohui Chen
Dept. of Statistics
UBC, Canada
_
Don't waste time standing in linetry shopping online. Visit Sympatico /
MSN
__
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.
_
Buy what you want when you want it on Sympatico / MSN Shopping
__
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 store recursive results
Hi Bálint, Thanks very much for your suggestions. The party package is a little bit complicated to use. Do you have a much simpler example for this problem? Anyway, I am gonna try this package. Cheers, Xiaohui Chen Dept. of Statistics UBC, Canada From: Bálint Czúcz [EMAIL PROTECTED] To: X.H Chen [EMAIL PROTECTED] CC: r-help@stat.math.ethz.ch Subject: Re: [R] how to store recursive results Date: Fri, 22 Sep 2006 12:21:02 +0200 One suggestion: as a recursive list. For example have a look at the 'party' package and the BinaryTree class. I hope this helps. Bálint On 22/09/06, X.H Chen [EMAIL PROTECTED] wrote: Hi all, How to store recursive resutls from a function for each step without using global operators -? Thanks ahead. Xiaohui Chen Dept. of Statistics UBC, Canada _ Don't waste time standing in linetry shopping online. Visit Sympatico / MSN __ 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] How to retrieve results of most recent command?
In R, is there an automatic variable that stores the results of the most
recent command or commands? (I am thinking of a behavior like
Mathematica's % result-history substitution syntax.)
(I am using R 2.3.1 on Linux and R 2.3.1 on Windows XP.)
This is a pretty basic question, so I tried to do an extensive version
of the recommended pre-posting homework.
help.search('substitution')
failed to turn up anything relevant. This is more an environmental
question, not an R-language programming question.
help.search('history')
was interesting, but focused on command history. 'savehistory' appears
to be an atomic function, so the temporary-file manipulation done by
'history' seems required. Perhaps I could try store the desired line in
a string and then evaluate it, but this method feels terribly clunky.
On the other hand, that may be the price to pay for avoidance of
creeping featurism.
The top 10 RSiteSearch results also focus on saving the command
histories to a file or script. The 'Emacs Speaks Statistics' (ESS)
package seems to offer shortcuts for command history expansion, but not
results history. I'm a bit reluctant to start using Emacs again, having
quit it for VI over a decade ago.
212-933-3305 / [EMAIL PROTECTED]
NOTICE TO RECIPIENTS: Any information contained in or attached
to this message is intended solely for the use of the intended
recipient(s). If you are not the intended recipient of this
transmittal, you are hereby notified that you received this
transmittal in error and we request that you please delete and
destroy all copies and attachments in your possession, notify
the sender that you have received this communication in error,
and note that any review or dissemination of, or the taking of
any action in reliance on, this communication is expressly
prohibited.
Banc of America Securities LLC (BAS) does not accept time
sensitive, action-oriented messages or transaction orders,
including orders to purchase or sell securities, via e-mail.
Regular internet e-mail transmission cannot be guaranteed to be
secure or error-free. Therefore, we do not represent that this
information is complete or accurate and it should not be relied
upon as such. If you prefer to communicate with BAS using secure
(i.e., encrypted) e-mail transmission, please notify the sender.
Otherwise, you will be deemed to have consented to communicate
with BAS via regular internet e-mail transmission. Please note
that BAS reserves the right to intercept, monitor, and retain all
e-mail messages (including secure e-mail messages) sent to or
from its systems as permitted by applicable law.
---
IRS Circular 230 Disclosure:
Bank of America Corporation and its affiliates, including BAS,
(Bank of America) do not provide tax advice. Accordingly, any
statements contained herein as to tax matters were neither written
nor intended by the sender or Bank of America to be used and cannot
be used by any taxpayer for the purpose of avoiding tax penalties
that may be imposed on such taxpayer. If any person uses or refers
to any such tax statement in promoting, marketing or recommending a
partnership or other entity, investment plan or arrangement to any
taxpayer, then the statement expressed above is being delivered to
support the promotion or marketing of the transaction or matter
addressed and the recipient should seek advice based on its
particular circumstances from an independent tax advisor.
__
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 retrieve results of most recent command?
Yeh, Richard C wrote: In R, is there an automatic variable that stores the results of the most recent command or commands? (I am thinking of a behavior like Mathematica's % result-history substitution syntax.) Something like .Last.value: x=sqrt(2) .Last.value [1] 1.414214 which might get more usage if it was less fiddly to type. Barry __ 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 store recursive results
1. The point is that you can use - and still not pollute the
global environment with any variables so its not clear
why there should be any requirement not to use it.
Other possiblities are
2. to pass the info around through the return value or through an environment:
fact - function(n) {
if (n == 1) list(n, c(1 = 1))
else {
f - fact(n-1)
out - list(n * f[[1]], unlist(f))
names(out[[2]])[1] - n
out
}
}
fact(4)
3. or through an environment:
fact - function(n, e) {
if (n == 1) e[[1]] - 1
else n * (e[[format(n)]] - fact(n-1, e))
}
e - new.env()
fact(4, e)
as.list(e)
On 9/22/06, X.H Chen [EMAIL PROTECTED] wrote:
Hi Gabor,
Thanks for pointing out this for me. However, what I try to get is how to
construct such form a function f that
ret-f(...),
where ret contains the each recursive result from f, and meantime f consists
of no - operator. Do you have any idea how to implemet this. Thanks a lot
for your suggestions.
Cheer
Xiaohui Chen
Dept. of Statistics
UBC, Canada
From: Gabor Grothendieck [EMAIL PROTECTED]
To: X.H Chen [EMAIL PROTECTED]
CC: r-help@stat.math.ethz.ch
Subject: Re: [R] how to store recursive results
Date: Fri, 22 Sep 2006 06:49:22 -0400
Note that - is not necessarily global:
if (exists(x)) rm(x)
f - function() {
x - 2
g - function() x - 3
g()
x
}
f() # 3
exists(x) # FALSE
On 9/22/06, X.H Chen [EMAIL PROTECTED] wrote:
Hi all,
How to store recursive resutls from a function for each step without using
global operators -? Thanks ahead.
Xiaohui Chen
Dept. of Statistics
UBC, Canada
_
Don't waste time standing in line—try shopping online. Visit Sympatico /
MSN
__
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.
_
Buy what you want when you want it on Sympatico / MSN Shopping
http://shopping.sympatico.msn.ca/content/shp/?ctId=2,ptnrid=176,ptnrdata=081805
__
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] Compiling a contingency table of counts by case
I have asked a similar question before but this time the problem is somewhat more involved. I have the following data: case;name;x 1;Joe;1 1;Mike;1 1;Zoe;1 2;Joe;1 2;Mike;0 2;Zoe;1 2;John;1 3;Mike;1 3;Zoe;0 3;Karl;0 I would like to count the number of case in which any two name a. both have x=1, b. the first has x=0 - the second has x=1, c. the first has x=1 - the second has x=0, d. both have x=0, One way is to use the reshape package: dat - read.delim(clipboard, sep=;) dat - rename(dat, c(x=value)) cast(dat, name ~ case) (which doesn't give you exactly what you want, but I think might be more illuminating) Hadley __ 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] extract data from lm object and then use again?
On Fri, 22 Sep 2006, Thilo Kellermann wrote: Hi, the data of the model fit is stored in lm$model and should work Not reliably. In the first place, you should use the accessor function model.frame(model) rather than model$model, which works even if the model was fitted with model=FALSE. But even then, glm(formula(model), data=model.frame(model)) will not work reliably. Consider model - lm(log(Volume)~log(Height)+log(Girth),data=trees) The model frame has variables called eg log(Volume) rather than Volume. When you need the source data frame you need to do something like eval(model$call$data, environment(formula(model))) and even this might not work, eg if the model had no data argument. However, if the model had no data argument then the variables must be available in environment(formula(model)), in which case any data frame of the right size will do. If there are no missing observations or the model was fitted with na.action=na.exclude then a fairly reliable approach is to use eval(model$call$data, environment(formula(model))) if it is not NULL and to fall back to model.frame(model). This is what termplot() does. -thomas __ 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 retrieve results of most recent command?
That's perfect! Thanks!
I now see that .Last.value is mentioned in the 'Introduction to R',
footnote 5. Sorry, I hadn't read it that carefully.
(Before posting, I even did read
help('.Last')
but thought that that function wasn't what I wanted.
I guess I wasn't using the right keywords.
help.search('last')
returns .Last.value as the top hit.
(I also apologize for the long disclaimer at the end. It's
automatically generated.)
212-933-3305 / [EMAIL PROTECTED]
-Original Message-
From: Barry Rowlingson [mailto:[EMAIL PROTECTED]
Sent: Friday, September 22, 2006 11:59 AM
To: Yeh, Richard C
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] How to retrieve results of most recent command?
Yeh, Richard C wrote:
In R, is there an automatic variable that stores the results of the
most
recent command or commands? (I am thinking of a behavior like
Mathematica's % result-history substitution syntax.)
Something like .Last.value:
x=sqrt(2)
.Last.value
[1] 1.414214
which might get more usage if it was less fiddly to type.
Barry
NOTICE TO RECIPIENTS: Any information contained in or attached
to this message is intended solely for the use of the intended
recipient(s). If you are not the intended recipient of this
transmittal, you are hereby notified that you received this
transmittal in error and we request that you please delete and
destroy all copies and attachments in your possession, notify
the sender that you have received this communication in error,
and note that any review or dissemination of, or the taking of
any action in reliance on, this communication is expressly
prohibited.
Banc of America Securities LLC (BAS) does not accept time
sensitive, action-oriented messages or transaction orders,
including orders to purchase or sell securities, via e-mail.
Regular internet e-mail transmission cannot be guaranteed to be
secure or error-free. Therefore, we do not represent that this
information is complete or accurate and it should not be relied
upon as such. If you prefer to communicate with BAS using secure
(i.e., encrypted) e-mail transmission, please notify the sender.
Otherwise, you will be deemed to have consented to communicate
with BAS via regular internet e-mail transmission. Please note
that BAS reserves the right to intercept, monitor, and retain all
e-mail messages (including secure e-mail messages) sent to or
from its systems as permitted by applicable law.
---
IRS Circular 230 Disclosure:
Bank of America Corporation and its affiliates, including BAS,
(Bank of America) do not provide tax advice. Accordingly, any
statements contained herein as to tax matters were neither written
nor intended by the sender or Bank of America to be used and cannot
be used by any taxpayer for the purpose of avoiding tax penalties
that may be imposed on such taxpayer. If any person uses or refers
to any such tax statement in promoting, marketing or recommending a
partnership or other entity, investment plan or arrangement to any
taxpayer, then the statement expressed above is being delivered to
support the promotion or marketing of the transaction or matter
addressed and the recipient should seek advice based on its
particular circumstances from an independent tax advisor.
__
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] Compiling a contingency table of counts by case
what's different from: with(dat, tapply(x, list(name,case), sum)) 1 2 3 Joe 1 1 NA John NA 1 NA Karl NA NA 0 Mike 1 0 1 and how to deal with this table ? --- Jacques VESLOT CNRS UMR 8090 I.B.L (2ème étage) 1 rue du Professeur Calmette B.P. 245 59019 Lille Cedex Tel : 33 (0)3.20.87.10.44 Fax : 33 (0)3.20.87.10.31 http://www-good.ibl.fr --- hadley wickham a écrit : I have asked a similar question before but this time the problem is somewhat more involved. I have the following data: case;name;x 1;Joe;1 1;Mike;1 1;Zoe;1 2;Joe;1 2;Mike;0 2;Zoe;1 2;John;1 3;Mike;1 3;Zoe;0 3;Karl;0 I would like to count the number of case in which any two name a. both have x=1, b. the first has x=0 - the second has x=1, c. the first has x=1 - the second has x=0, d. both have x=0, One way is to use the reshape package: dat - read.delim(clipboard, sep=;) dat - rename(dat, c(x=value)) cast(dat, name ~ case) (which doesn't give you exactly what you want, but I think might be more illuminating) Hadley __ 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] Exponentiate a matrix
I am getting a bit rusty on some of these things, but I seem to recall
that there is a numerical advantage (speed and/or accuracy?) to
diagonalizing:
expM - function(X,e) { v - La.svd(X); v$u %*% diag(v$d^e) %*% v$vt }
P - matrix(c(.3,.7, .7, .3), ncol=2)
P %*% P %*% P
[,1] [,2]
[1,] 0.468 0.532
[2,] 0.532 0.468
expM(P,3)
[,1] [,2]
[1,] 0.468 0.532
[2,] 0.532 0.468
I think this also works for non-integer, negative, large, and complex
exponents:
expM(P, 1.5)
[,1] [,2]
[1,] 0.3735089 0.6264911
[2,] 0.6264911 0.3735089
expM(P, 1000)
[,1] [,2]
[1,] 0.5 0.5
[2,] 0.5 0.5
expM(P, -3)
[,1][,2]
[1,] -7.3125 8.3125
[2,] 8.3125 -7.3125
expM(P, 3+.5i)
[,1] [,2]
[1,] 0.4713+0.0141531i 0.5287-0.0141531i
[2,] 0.5287-0.0141531i 0.4713+0.0141531i
Paul Gilbert
Doran, Harold wrote:
Suppose I have a square matrix P
P - matrix(c(.3,.7, .7, .3), ncol=2)
I know that
P * P
Returns the element by element product, whereas
P%*%P
Returns the matrix product.
Now, P^2 also returns the element by element product. But, is there a
slick way to write
P %*% P %*% P
Obviously, P^3 does not return the result I expect.
Thanks,
Harold
[[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.
La version française suit le texte anglais.
This email may contain privileged and/or confidential inform...{{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 store recursive results
On Fri, 22 Sep 2006, X.H Chen wrote:
Hi Gabor,
Thanks for pointing out this for me. However, what I try to get is how to
construct such form a function f that
ret-f(...),
where ret contains the each recursive result from f, and meantime f consists
of no - operator. Do you have any idea how to implemet this. Thanks a lot
for your suggestions.
It depends on the situation. You can always pass the results back in a
list or vector, eg
cumfactorial-function(n){
if (n==0)
1
else c(1, n*cumfactorial(n-1))
}
If you want to get the results out then you have to either accumulate and
return them like this or use -, since return() and - are the only ways
to get results out of a function. As long as you don't use - to assign
to variables outside the function it is a perfectly reasonable thing to do
If you were doing something like a Fibonacci sequence then assigning would
be preferable
fib-function(n){
memo-new.env(hash=TRUE)
fibrec-function(m){
if (m=2) return(1)
vm-paste(v,m,sep=)
if(exists(vm,envir=memo,inherits=FALSE))
return(get(vm,envir=memo,inherits=FALSE))
rval-fibrec(m-1)+fibrec(m-2)
assign(vm,rval,envir=memo)
rval
}
fibrec(n)
sapply(ls(envir=memo),get, envir=memo)
}
since the memoization converts the algorithm from exponential time to
linear time.
__
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] Compiling a contingency table of counts by case
what's different from: with(dat, tapply(x, list(name,case), sum)) 1 2 3 Joe 1 1 NA John NA 1 NA Karl NA NA 0 Mike 1 0 1 and how to deal with this table ? Well, the syntax is easier (once you have the data in the correct, molten, form), and more flexible for other tasks. It is surely better to learn a general purpose tool rather than a tool for a specific task. To use that table to answer the original question, you just need to look column by column for the desired patterns of 0's and 1's, eg. in case 1, Joe, Mike and Zoe all had ones. Perhaps I misunderstood the original question. Hadley __ 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] extract data from lm object and then use again?
On Fri, 22 Sep 2006, Thomas Lumley wrote: On Fri, 22 Sep 2006, Thilo Kellermann wrote: Hi, the data of the model fit is stored in lm$model and should work Not reliably. In the first place, you should use the accessor function model.frame(model) rather than model$model, which works even if the model was fitted with model=FALSE. But even then, glm(formula(model), data=model.frame(model)) will not work reliably. Consider model - lm(log(Volume)~log(Height)+log(Girth),data=trees) The model frame has variables called eg log(Volume) rather than Volume. When you need the source data frame you need to do something like eval(model$call$data, environment(formula(model))) and even this might not work, eg if the model had no data argument. However, if the model had no data argument then the variables must be available in environment(formula(model)), in which case any data frame of the right size will do. to be picky ... must have been available. You or some other command could very easily have changed them. That's actually why we store the model frame by default: there is no other 100% reliable way to get at the data used in the fitting (as distinct from the data originally supplied). If there are no missing observations or the model was fitted with na.action=na.exclude then a fairly reliable approach is to use eval(model$call$data, environment(formula(model))) if it is not NULL and to fall back to model.frame(model). This is what termplot() does. -- 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] $theta of frailty in coxph
Dear all, Does the frailty.object$history[[1]]$theta returns the Variance of random effect? Why is the value different? Here is an example with kidney data: library(survival) data(kidney) frailty.object-coxph(Surv(time, status)~ age + sex + disease + frailty(id), kidney) frailty.object Call: coxph(formula = Surv(time, status) ~ age + sex + disease + frailty(id), data = kidney) coef se(coef) se2Chisq DF p age 0.00318 0.0111 0.0111 0.08 1 7.8e-01 sex -1.48314 0.3582 0.3582 17.14 1 3.5e-05 diseaseGN0.08796 0.4064 0.4064 0.05 1 8.3e-01 diseaseAN0.35079 0.3997 0.3997 0.77 1 3.8e-01 diseasePKD -1.43111 0.6311 0.6311 5.14 1 2.3e-02 frailty(id) 0.00 0 9.3e-01 Iterations: 6 outer, 28 Newton-Raphson Variance of random effect= 5e-07 I-likelihood = -179.1 Degrees of freedom for terms= 1 1 3 0 Likelihood ratio test=17.6 on 5 df, p=0.00342 n= 76 frailty.object$history[[1]]$theta [1] 5e-09 Thanks in advance. Mohammad Ehsanul Karim Statistical Research and Training, University of Dhaka [[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] Re : Statitics Textbook - any recommendation?
Excuses me Dear Martin, I have make a mistake when posting my message ! Sincerly ! - Message d'origine De : Martin Maechler [EMAIL PROTECTED] À : justin bem [EMAIL PROTECTED] Cc : Martin Maechler [EMAIL PROTECTED] Envoyé le : Vendredi, 22 Septembre 2006, 7h39mn 25s Objet : Re: Re : [R] Statitics Textbook - any recommendation? justin == justin bem [EMAIL PROTECTED] on Thu, 21 Sep 2006 19:00:24 + (GMT) writes: justin Why need to learn Statistic in R in the same book ? justin I think it is not very efficient. you can use a book justin to learn how to use R in statistic, and you get justin statistic skills by reading many others books. Good point. Why are you sending this only to me rather than back to R-help? Bonnes salutations de Zurich! Martin justin - Message d'origine De : Martin Maechler justin [EMAIL PROTECTED] À : Iuri Gavronski justin [EMAIL PROTECTED]; [EMAIL PROTECTED] Cc : justin [EMAIL PROTECTED]; justin r-help@stat.math.ethz.ch Envoyé le : Jeudi, 21 justin Septembre 2006, 9h16mn 21s Objet : Re: [R] Statitics justin Textbook - any recommendation? GS == Gavin Simpson [EMAIL PROTECTED] on Thu, 21 Sep 2006 00:08:17 +0100 writes: GS On Wed, 2006-09-20 at 18:56 -0400, Charles Annis, GS P.E. wrote: Recommending a good book on statistics is like recommending a good book on sports: Which sports? A good book for learning statistical concepts (and learning R at the same time), one that assumes you understand algebra but are new to statistics, is Peter Dalgaard's _Introductory Statistics with R_ (Springer 2002). The writing is relaxed and succinct, not condescending as some texts might appear to a newcomer. It's just a good book. GS I couldn't agree more. A number of my colleagues have GS bought Peter Dalgaard's book to a) learn some R and b) GS learn some statistics. They have found it very useful GS indeed. justin Yes! I'm pretty sure it has been the first book of justin its kind (Intro Stats + R), and in my view is justin still the best. justin Martin Maechler, ETH Zurich Charles Annis, P.E. [EMAIL PROTECTED] phone: 561-352-9699 eFax: 614-455-3265 http://www.StatisticalEngineering.com -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Mike Nielsen Sent: Wednesday, September 20, 2006 6:36 PM To: Berton Gunter Cc: r-help@stat.math.ethz.ch Subject: Re: [R] Statitics Textbook - any recommendation? Excellent characterization. MASS is a very good book, but I'm not sure I would describe it as a statistics textbook, much less one of the basic variety. While I certainly wouldn't presume to speak for Prof. Ripley and Dr. Venables, it seems unlikely their intent in writing MASS was to teach statistics, but rather, as the name of the book might suggest, to explain how S+ (and R) can be applied to modern statistical techniques. My experience with this book is that it assumes considerable background knowledge. By all means, buy MASS, but if you need guidance on the how and why of statistical techniques, you may wish to shop Amazon to find a supplement. Regards, Mike On 9/20/06, Berton Gunter [EMAIL PROTECTED] wrote: Not withstanding Prof. Heiberger's admirable enthusiasm, I think the canonical answer is probably MASS (Modern Applied Statistics with S) by Venables and Ripley. It is very comprehensive, but depending on your background, you may find it too telegraphic. -- Bert Gunter Genentech Non-Clinical Statistics South San Francisco, CA The business of the statistician is to catalyze the scientific learning process. - George E. P. Box -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Iuri Gavronski Sent: Wednesday, September 20, 2006 1:22 PM To: r-help@stat.math.ethz.ch Subject: [R] Statitics Textbook - any recommendation? I would like to buy a basic statistics book (experimental design, sampling, ANOVA, regression, etc.) with examples in R. Or download it in PDF or html format.I went to the CRAN contributed documentation, but there were only R textbooks, that is, textbooks where R is the focus, not the statistics. And I would like to find the opposite. Other text I am trying to find is multivariate data analysis (EFA, cluster, mult regression, MANOVA, etc.) with examples with R.Any recommendation? Thank you in advance, Iuri. __ R-help@stat.math.ethz.ch mailing list
[R] inequality with NA
Dear everybody!
take a-c(5,3,NA,6).
if(a[1]!=NA){b-7}
if(a[3]!=5){b-7}
if(a[3]!=NA){b-7}
if(a[3]==NA){b-7}
will alltogeather return
Fehler in if (a[1] != NA) { : Fehlender Wert, wo TRUE/FALSE nötig ist
(or simularly). Somehow this is logical. But how else should I get out,
whether a certain vector-component has an existing value?
Thank you in advance!
Yours,
Mag. Ferri Leberl
__
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] inequality with NA
Mag. Ferri Leberl wrote:
Dear everybody!
take a-c(5,3,NA,6).
if(a[1]!=NA){b-7}
if(a[3]!=5){b-7}
if(a[3]!=NA){b-7}
if(a[3]==NA){b-7}
will alltogeather return
Fehler in if (a[1] != NA) { : Fehlender Wert, wo TRUE/FALSE nötig ist
(or simularly). Somehow this is logical. But how else should I get out,
whether a certain vector-component has an existing value?
see ?is.na
Hence:
if(!is.na(a[1])) b - 7
Uwe Ligges
Thank you in advance!
Yours,
Mag. Ferri Leberl
__
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] inequality with NA
On Fri, 2006-09-22 at 20:16 +0200, Mag. Ferri Leberl wrote:
Dear everybody!
take a-c(5,3,NA,6).
if(a[1]!=NA){b-7}
if(a[3]!=5){b-7}
if(a[3]!=NA){b-7}
if(a[3]==NA){b-7}
will alltogeather return
Fehler in if (a[1] != NA) { : Fehlender Wert, wo TRUE/FALSE nötig ist
(or simularly). Somehow this is logical. But how else should I get out,
whether a certain vector-component has an existing value?
Thank you in advance!
Yours,
Mag. Ferri Leberl
NA is not defined, so you cannot predictably perform equality/inequality
tests with it. There are specific functions in place for dealing with
this.
See ?is.na and ?na.omit
a
[1] 5 3 NA 6
a[is.na(a)]
[1] NA
a[!is.na(a)]
[1] 5 3 6
You can also use which() to find the indices:
which(is.na(a))
[1] 3
which(!is.na(a))
[1] 1 2 4
Finally, use na.omit() to remove all NA's:
na.omit(a)
[1] 5 3 6
attr(,na.action)
[1] 3
attr(,class)
[1] omit
Note that the object attribute 'na.action' shows that a[3] was removed:
a.omit - na.omit(a)
as.vector(attr(a.omit, na.action))
[1] 3
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] behavior of [-.foo
Can someone help me understand the following behavior of [- ?
If I define a simple class based on a matrix, the [- operation only
inserts into the first column:
x - matrix(rnorm(10),nrow=5,ncol=2)
class(x) - foo
[-.foo - function(x, i, j, value) {
+ if(missing(i)) i - 1:nrow(x)
+ if(missing(j)) j - 1:ncol(x)
+
+ x - unclass(x)
+ x - NextMethod(.Generic)
+ class(x) - foo
+ x
+ }
x[] - 100.0
x
[,1] [,2]
[1,] 100 -0.1465296
[2,] 100 -0.2615796
[3,] 100 -0.8882629
[4,] 100 -0.2886357
[5,] 100 -0.9565273
attr(,class)
[1] foo
Based on the behavior of [- for a matrix, I would have thought that the
data for the whole object would be replaced.
for instance:
y - matrix(rnorm(10),nrow=5,ncol=2)
y
[,1] [,2]
[1,] -0.55297049 -1.1896488
[2,] 0.06157438 -0.6628254
[3,] -0.28184208 -2.5260177
[4,] 0.61204398 -0.3492488
[5,] 0.43971216 1.8990789
y[] - 100
y
[,1] [,2]
[1,] 100 100
[2,] 100 100
[3,] 100 100
[4,] 100 100
[5,] 100 100
Thanks,
Whit
code for above:
x - matrix(rnorm(10),nrow=5,ncol=2)
x
class(x) - foo
[-.foo - function(x, i, j, value) {
if(missing(i)) i - 1:nrow(x)
if(missing(j)) j - 1:ncol(x)
x - unclass(x)
x - NextMethod(.Generic)
class(x) - foo
x
}
x[] - 100.0
x
R.Version()
$platform
[1] i686-pc-linux-gnu
$arch
[1] i686
$os
[1] linux-gnu
$system
[1] i686, linux-gnu
$status
[1]
$major
[1] 2
$minor
[1] 3.1
$year
[1] 2006
$month
[1] 06
$day
[1] 01
$`svn rev`
[1] 38247
$language
[1] R
$version.string
[1] Version 2.3.1 (2006-06-01)
This e-mail message is intended only for the named recipient(s) above. It may
contain confidential information. If you are not the intended recipient you are
hereby notified that any dissemination, distribution or copying of this e-mail
and any attachment(s) is strictly prohibited. If you have received this e-mail
in error, please immediately notify the sender by replying to this e-mail and
delete the message and any attachment(s) from your system. Thank you.
__
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] Creating Movies with R
If you run R on Linux, then you can run the ImageMagick command called
convert. I place this in an R function to use a sequence of PNG plots as
movie frames:
make.mov.plotcol3d - function(){
unlink(plotcol3d.mpg)
system(convert -delay 10 plotcol3d*.png plotcol3d.mpg)
}
Examples can be seen here:
http://biostat.mc.vanderbilt.edu/JrhRgbColorSpace
Look for the 'Download Movie' links.
Cheers,
Jeff
Lorenzo Isella wrote:
Dear All,
I'd like to know if it is possible to create animations with R.
To be specific, I attach a code I am using for my research to plot
some analytical results in 3D using the lattice package. It is not
necessary to go through the code.
Simply, it plots some 3D density profiles at two different times
selected by the user.
I wonder if it is possible to use the data generated for different
times to create something like an .avi file.
Here is the script:
rm(list=ls())
library(lattice)
# I start defining the analytical functions needed to get the density
as a function of time
expect_position - function(t,lam1,lam2,pos_ini,vel_ini)
{1/(lam1-lam2)*(lam1*exp(lam2*t)-lam2*exp(lam1*t))*pos_ini+
1/(lam1-lam2)*(exp(lam1*t)-exp(lam2*t))*vel_ini
}
sigma_pos-function(t,q,lam1,lam2)
{
q/(lam1-lam2)^2*(
(exp(2*lam1*t)-1)/(2*lam1)-2/(lam1+lam2)*(exp(lam1*t+lam2*t)-1) +
(exp(2*lam2*t)-1)/(2*lam2) )
}
rho_x-function(x,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(x-expect_position)^2/sigma_pos)
}
Now the physical parameters
tau-0.1
beta-1/tau
St-tau ### since I am in dimensionless units and tau is already in
units of 1/|alpha|
D=2e-2
q-2*beta^2*D
### Now the grid in space and time
time-5 # time extent
tsteps-501 # time steps
newtime-seq(0,time,len=tsteps)
Now the things specific for the dynamics along x
lam1- -beta/2*(1+sqrt(1+4*St))
lam2- -beta/2*(1-sqrt(1+4*St))
xmin- -0.5
xmax-0.5
x0-0.1
vx0-x0
nx-101 ## grid intervals along x
newx-seq(xmin,xmax,len=nx) # grid along x
# M1 - do.call(g, c(list(x = newx), mypar))
mypar-c(q,lam1,lam2)
sig_xx-do.call(sigma_pos,c(list(t=newtime),mypar))
mypar-c(lam1,lam2,x0,vx0)
exp_x-do.call(expect_position,c(list(t=newtime),mypar))
#rho_x-function(x,expect_position,sigma_pos)
#NB: at t=0, the density blows up, since I have a delta as the initial state!
# At any t0, instead, the result is finite.
#for this reason I now redefine time by getting rid of the istant t=0
to work out
# the density
rho_x_t-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_x[i],sig_xx[i])
myrho_x-do.call(rho_x,c(list(x=newx),mypar))
rho_x_t[ i-1, ]-myrho_x
}
### Now I also define a scaled density
rho_x_t_scaled-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_x[i],sig_xx[i])
myrho_x-do.call(rho_x,c(list(x=newx),mypar))
rho_x_t_scaled[ i-1, ]-myrho_x/max(myrho_x)
}
###Now I deal with the dynamics along y
lam1- -beta/2*(1+sqrt(1-4*St))
lam2- -beta/2*(1-sqrt(1-4*St))
ymin- 0
ymax- 1
y0-ymax
vy0- -y0
mypar-c(q,lam1,lam2)
sig_yy-do.call(sigma_pos,c(list(t=newtime),mypar))
mypar-c(lam1,lam2,y0,vy0)
exp_y-do.call(expect_position,c(list(t=newtime),mypar))
# now I introduce the function giving the density along y: this has to
include the BC of zero
# density at wall
rho_y-function(y,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y-expect_position)^2/sigma_pos)-
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y+expect_position)^2/sigma_pos)
}
newy-seq(ymin,ymax,len=nx) # grid along y with the same # of points
as the one along x
rho_y_t-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_y[i],sig_yy[i])
myrho_y-do.call(rho_y,c(list(y=newy),mypar))
rho_y_t[ i-1, ]-myrho_y
}
rho_y_t_scaled-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_y[i],sig_yy[i])
myrho_y-do.call(rho_y,c(list(y=newy),mypar))
rho_y_t_scaled[ i-1, ]-myrho_y/max(myrho_y)
}
# The following 2 plots are an example of the plots I'd like to use to
make an animation
g - expand.grid(x = newx, y = newy)
instant-100
mydens-rho_x_t[ instant, ]%o%rho_y_t[ instant, ]/(max(rho_x_t[
instant, ]%o%rho_y_t[ instant, ]))
lentot-nx^2
dim(mydens)-c(lentot,1)
g$z-mydens
jpeg(dens-t-3.jpeg)
print(wireframe(z ~ x * y, g, drape = TRUE,shade=TRUE,
scales = list(arrows = FALSE),pretty=FALSE, aspect = c(1,1), colorkey = TRUE
,zoom=0.8, main=expression(Density at t=2), zlab =
list(expression(density),rot = 90),distance=0.0,
perspective=TRUE,#screen = list(z = 150, x = -55,y= 0)
,zlim=range(c(0,1
dev.off()
instant-300
mydens-rho_x_t[ instant, ]%o%rho_y_t[ instant, ]/(max(rho_x_t[
instant, ]%o%rho_y_t[ instant, ]))
lentot-nx^2
dim(mydens)-c(lentot,1)
g$z-mydens
jpeg(dens-t-3.jpeg)
print(wireframe(z ~ x * y, g, drape = TRUE,shade=TRUE,
scales = list(arrows = FALSE),pretty=FALSE, aspect = c(1,1), colorkey = TRUE
,zoom=0.8, main=expression(Density
Re: [R] behavior of [-.foo
Try this:
x - matrix(rnorm(10),nrow=5,ncol=2)
class(x) - foo
[-.foo - function(x, i = TRUE, j = TRUE, ..., value) {
x - unclass(x)
x - NextMethod()
class(x) - foo
x
}
x[] - 100.0
On 9/22/06, Armstrong, Whit [EMAIL PROTECTED] wrote:
Can someone help me understand the following behavior of [- ?
If I define a simple class based on a matrix, the [- operation only
inserts into the first column:
x - matrix(rnorm(10),nrow=5,ncol=2)
class(x) - foo
[-.foo - function(x, i, j, value) {
+ if(missing(i)) i - 1:nrow(x)
+ if(missing(j)) j - 1:ncol(x)
+
+ x - unclass(x)
+ x - NextMethod(.Generic)
+ class(x) - foo
+ x
+ }
x[] - 100.0
x
[,1] [,2]
[1,] 100 -0.1465296
[2,] 100 -0.2615796
[3,] 100 -0.8882629
[4,] 100 -0.2886357
[5,] 100 -0.9565273
attr(,class)
[1] foo
Based on the behavior of [- for a matrix, I would have thought that the
data for the whole object would be replaced.
for instance:
y - matrix(rnorm(10),nrow=5,ncol=2)
y
[,1] [,2]
[1,] -0.55297049 -1.1896488
[2,] 0.06157438 -0.6628254
[3,] -0.28184208 -2.5260177
[4,] 0.61204398 -0.3492488
[5,] 0.43971216 1.8990789
y[] - 100
y
[,1] [,2]
[1,] 100 100
[2,] 100 100
[3,] 100 100
[4,] 100 100
[5,] 100 100
Thanks,
Whit
code for above:
x - matrix(rnorm(10),nrow=5,ncol=2)
x
class(x) - foo
[-.foo - function(x, i, j, value) {
if(missing(i)) i - 1:nrow(x)
if(missing(j)) j - 1:ncol(x)
x - unclass(x)
x - NextMethod(.Generic)
class(x) - foo
x
}
x[] - 100.0
x
R.Version()
$platform
[1] i686-pc-linux-gnu
$arch
[1] i686
$os
[1] linux-gnu
$system
[1] i686, linux-gnu
$status
[1]
$major
[1] 2
$minor
[1] 3.1
$year
[1] 2006
$month
[1] 06
$day
[1] 01
$`svn rev`
[1] 38247
$language
[1] R
$version.string
[1] Version 2.3.1 (2006-06-01)
This e-mail message is intended only for the named recipient(s) above. It may
contain confidential information. If you are not the intended recipient you
are hereby notified that any dissemination, distribution or copying of this
e-mail and any attachment(s) is strictly prohibited. If you have received
this e-mail in error, please immediately notify the sender by replying to
this e-mail and delete the message and any attachment(s) from your system.
Thank you.
__
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] Creating Movies with R
An alternative that I've used a few times is the jpg() function to
create the sequence of images, and then converting these to an mpeg
movie using mencoder distributed with mplayer. This works on both
windows and linux. I have a pretty self-contained example file
written up that I can send to anyone who is interested. Oddly, the
most challenging part was creating a sequence of file names that would
be correctly ordered - for this I use:
lex - function(N){
## produce vector of N lexicograpically ordered strings
ndig - nchar(N)
substr(formatC((1:N)/10^ndig,digits=ndig,format=f),3,1000)
}
On Fri, 22 Sep 2006, Jeffrey Horner wrote:
Date: Fri, 22 Sep 2006 13:46:52 -0500
From: Jeffrey Horner [EMAIL PROTECTED]
To: Lorenzo Isella [EMAIL PROTECTED], r-help@stat.math.ethz.ch
Subject: Re: [R] Creating Movies with R
If you run R on Linux, then you can run the ImageMagick command called
convert. I place this in an R function to use a sequence of PNG plots as
movie frames:
make.mov.plotcol3d - function(){
unlink(plotcol3d.mpg)
system(convert -delay 10 plotcol3d*.png plotcol3d.mpg)
}
Examples can be seen here:
http://biostat.mc.vanderbilt.edu/JrhRgbColorSpace
Look for the 'Download Movie' links.
Cheers,
Jeff
Lorenzo Isella wrote:
Dear All,
I'd like to know if it is possible to create animations with R.
To be specific, I attach a code I am using for my research to plot
some analytical results in 3D using the lattice package. It is not
necessary to go through the code.
Simply, it plots some 3D density profiles at two different times
selected by the user.
I wonder if it is possible to use the data generated for different
times to create something like an .avi file.
Here is the script:
rm(list=ls())
library(lattice)
# I start defining the analytical functions needed to get the density
as a function of time
expect_position - function(t,lam1,lam2,pos_ini,vel_ini)
{1/(lam1-lam2)*(lam1*exp(lam2*t)-lam2*exp(lam1*t))*pos_ini+
1/(lam1-lam2)*(exp(lam1*t)-exp(lam2*t))*vel_ini
}
sigma_pos-function(t,q,lam1,lam2)
{
q/(lam1-lam2)^2*(
(exp(2*lam1*t)-1)/(2*lam1)-2/(lam1+lam2)*(exp(lam1*t+lam2*t)-1) +
(exp(2*lam2*t)-1)/(2*lam2) )
}
rho_x-function(x,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(x-expect_position)^2/sigma_pos)
}
Now the physical parameters
tau-0.1
beta-1/tau
St-tau ### since I am in dimensionless units and tau is already in
units of 1/|alpha|
D=2e-2
q-2*beta^2*D
### Now the grid in space and time
time-5 # time extent
tsteps-501 # time steps
newtime-seq(0,time,len=tsteps)
Now the things specific for the dynamics along x
lam1- -beta/2*(1+sqrt(1+4*St))
lam2- -beta/2*(1-sqrt(1+4*St))
xmin- -0.5
xmax-0.5
x0-0.1
vx0-x0
nx-101 ## grid intervals along x
newx-seq(xmin,xmax,len=nx) # grid along x
# M1 - do.call(g, c(list(x = newx), mypar))
mypar-c(q,lam1,lam2)
sig_xx-do.call(sigma_pos,c(list(t=newtime),mypar))
mypar-c(lam1,lam2,x0,vx0)
exp_x-do.call(expect_position,c(list(t=newtime),mypar))
#rho_x-function(x,expect_position,sigma_pos)
#NB: at t=0, the density blows up, since I have a delta as the initial
state!
# At any t0, instead, the result is finite.
#for this reason I now redefine time by getting rid of the istant t=0
to work out
# the density
rho_x_t-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_x[i],sig_xx[i])
myrho_x-do.call(rho_x,c(list(x=newx),mypar))
rho_x_t[ i-1, ]-myrho_x
}
### Now I also define a scaled density
rho_x_t_scaled-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_x[i],sig_xx[i])
myrho_x-do.call(rho_x,c(list(x=newx),mypar))
rho_x_t_scaled[ i-1, ]-myrho_x/max(myrho_x)
}
###Now I deal with the dynamics along y
lam1- -beta/2*(1+sqrt(1-4*St))
lam2- -beta/2*(1-sqrt(1-4*St))
ymin- 0
ymax- 1
y0-ymax
vy0- -y0
mypar-c(q,lam1,lam2)
sig_yy-do.call(sigma_pos,c(list(t=newtime),mypar))
mypar-c(lam1,lam2,y0,vy0)
exp_y-do.call(expect_position,c(list(t=newtime),mypar))
# now I introduce the function giving the density along y: this has to
include the BC of zero
# density at wall
rho_y-function(y,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y-expect_position)^2/sigma_pos)-
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y+expect_position)^2/sigma_pos)
}
newy-seq(ymin,ymax,len=nx) # grid along y with the same # of points
as the one along x
rho_y_t-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_y[i],sig_yy[i])
myrho_y-do.call(rho_y,c(list(y=newy),mypar))
rho_y_t[ i-1, ]-myrho_y
}
rho_y_t_scaled-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_y[i],sig_yy[i])
myrho_y-do.call(rho_y,c(list(y=newy),mypar))
rho_y_t_scaled[ i-1, ]-myrho_y/max(myrho_y)
}
# The following 2
Re: [R] Creating Movies with R
See the flag= argument on formatC:
n - 10
formatC(1:n, dig = nchar(n)-1, flag = 0)
# Here is another way:
n - 10
sprintf(paste(%0, nchar(n), .0f, sep = ), 1:n)
On 9/22/06, J.R. Lockwood [EMAIL PROTECTED] wrote:
An alternative that I've used a few times is the jpg() function to
create the sequence of images, and then converting these to an mpeg
movie using mencoder distributed with mplayer. This works on both
windows and linux. I have a pretty self-contained example file
written up that I can send to anyone who is interested. Oddly, the
most challenging part was creating a sequence of file names that would
be correctly ordered - for this I use:
lex - function(N){
## produce vector of N lexicograpically ordered strings
ndig - nchar(N)
substr(formatC((1:N)/10^ndig,digits=ndig,format=f),3,1000)
}
On Fri, 22 Sep 2006, Jeffrey Horner wrote:
Date: Fri, 22 Sep 2006 13:46:52 -0500
From: Jeffrey Horner [EMAIL PROTECTED]
To: Lorenzo Isella [EMAIL PROTECTED], r-help@stat.math.ethz.ch
Subject: Re: [R] Creating Movies with R
If you run R on Linux, then you can run the ImageMagick command called
convert. I place this in an R function to use a sequence of PNG plots as
movie frames:
make.mov.plotcol3d - function(){
unlink(plotcol3d.mpg)
system(convert -delay 10 plotcol3d*.png plotcol3d.mpg)
}
Examples can be seen here:
http://biostat.mc.vanderbilt.edu/JrhRgbColorSpace
Look for the 'Download Movie' links.
Cheers,
Jeff
Lorenzo Isella wrote:
Dear All,
I'd like to know if it is possible to create animations with R.
To be specific, I attach a code I am using for my research to plot
some analytical results in 3D using the lattice package. It is not
necessary to go through the code.
Simply, it plots some 3D density profiles at two different times
selected by the user.
I wonder if it is possible to use the data generated for different
times to create something like an .avi file.
Here is the script:
rm(list=ls())
library(lattice)
# I start defining the analytical functions needed to get the density
as a function of time
expect_position - function(t,lam1,lam2,pos_ini,vel_ini)
{1/(lam1-lam2)*(lam1*exp(lam2*t)-lam2*exp(lam1*t))*pos_ini+
1/(lam1-lam2)*(exp(lam1*t)-exp(lam2*t))*vel_ini
}
sigma_pos-function(t,q,lam1,lam2)
{
q/(lam1-lam2)^2*(
(exp(2*lam1*t)-1)/(2*lam1)-2/(lam1+lam2)*(exp(lam1*t+lam2*t)-1) +
(exp(2*lam2*t)-1)/(2*lam2) )
}
rho_x-function(x,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(x-expect_position)^2/sigma_pos)
}
Now the physical parameters
tau-0.1
beta-1/tau
St-tau ### since I am in dimensionless units and tau is already in
units of 1/|alpha|
D=2e-2
q-2*beta^2*D
### Now the grid in space and time
time-5 # time extent
tsteps-501 # time steps
newtime-seq(0,time,len=tsteps)
Now the things specific for the dynamics along x
lam1- -beta/2*(1+sqrt(1+4*St))
lam2- -beta/2*(1-sqrt(1+4*St))
xmin- -0.5
xmax-0.5
x0-0.1
vx0-x0
nx-101 ## grid intervals along x
newx-seq(xmin,xmax,len=nx) # grid along x
# M1 - do.call(g, c(list(x = newx), mypar))
mypar-c(q,lam1,lam2)
sig_xx-do.call(sigma_pos,c(list(t=newtime),mypar))
mypar-c(lam1,lam2,x0,vx0)
exp_x-do.call(expect_position,c(list(t=newtime),mypar))
#rho_x-function(x,expect_position,sigma_pos)
#NB: at t=0, the density blows up, since I have a delta as the initial
state!
# At any t0, instead, the result is finite.
#for this reason I now redefine time by getting rid of the istant t=0
to work out
# the density
rho_x_t-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_x[i],sig_xx[i])
myrho_x-do.call(rho_x,c(list(x=newx),mypar))
rho_x_t[ i-1, ]-myrho_x
}
### Now I also define a scaled density
rho_x_t_scaled-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar-c(exp_x[i],sig_xx[i])
myrho_x-do.call(rho_x,c(list(x=newx),mypar))
rho_x_t_scaled[ i-1, ]-myrho_x/max(myrho_x)
}
###Now I deal with the dynamics along y
lam1- -beta/2*(1+sqrt(1-4*St))
lam2- -beta/2*(1-sqrt(1-4*St))
ymin- 0
ymax- 1
y0-ymax
vy0- -y0
mypar-c(q,lam1,lam2)
sig_yy-do.call(sigma_pos,c(list(t=newtime),mypar))
mypar-c(lam1,lam2,y0,vy0)
exp_y-do.call(expect_position,c(list(t=newtime),mypar))
# now I introduce the function giving the density along y: this has to
include the BC of zero
# density at wall
rho_y-function(y,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y-expect_position)^2/sigma_pos)-
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y+expect_position)^2/sigma_pos)
}
newy-seq(ymin,ymax,len=nx) # grid along y with the same # of points
as the one along x
rho_y_t-matrix(ncol=nx,nrow=tsteps-1)
for (i in
[R] logistic + neg binomial + ...
Hi Folks, I've just come across a kind of problem which leads me to wonder how to approach it in R. Basically, each a set of items is subjected to a series of impacts until it eventually fails. The force of each impact would depend on covariates X,Y say; but as a result of preceding impacts an item would be expected to have a cumulative frailty such that the probability of failure due to a particular impact would possibly increase according to the series of impacts already survived. Without the cumulative frailty one could envisage something like a logistic model for the probabiliy of failure at each impact, leading to a kind of generalised exponential distribution -- that is, the likelihood for each item would be of the form (1-P[1])*(1-P[2])*...*(1-P[n-1])*P[n] where P[i] could have a logistic model in terms of the values of X[i] and Y[i], and n is the index of the impact at which failure occurred. That is then a solvable problem. Even so, I'm not (so far) finding in the R resources the appropriate analogue of glm for this kind of model. I dare say a prolonged trawl through the various survival resources might lead to something applicable, but ... And then there's the cumulative frailty ... ! Suggestions welcome! With thanks, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 22-Sep-06 Time: 20:25:12 -- 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] Lattice strip labels for two factors
Dear list, My problem is to change the strip text of lattice panels when using two factors. I have a data frame with two factors: df - expand.grid( fact1=c(y,b,r), fact2=c(far,por,lis,set), year=1991:2000, value= NA) df[,value] - sample(1:50, 120, replace=TRUE) I can make simple xyplot and change the text of the factor levels with strip.custom: require(lattice) xyplot( value ~ year | fact1, data=df, type=b, subset= fact2==far, strip = strip.custom(bg=gray.colors(1,0.95), factor.levels=c(yellow, black, red)), layout=c(1,3) ) But how can I change the text of the factor levels when using both factors as in this plot: xyplot( value ~ year | fact1*fact2, data=df, type=b) (fact2 levels text should change to: c(faro,porto,lisbon,setubal)) I read the help for strip.default and the emails archive, tried with which.given but could not find out how to accomplish this. Many thanks, Rafael Duarte -- Rafael Duarte Marine Resources Department - DRM IPIMAR - National Research Institute for Agriculture and Fisheries Av. Brasília, 1449-006 Lisbon - Portugal Tel:+351 21 302 7000 Fax:+351 21 301 5948 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] Lattice strip labels for two factors
Try this: levels(df$fact2) - c(faro,porto,lisbon,setubal) xyplot( value ~ year | fact1*fact2, data=df, type=b) On 9/22/06, Rafael Duarte [EMAIL PROTECTED] wrote: Dear list, My problem is to change the strip text of lattice panels when using two factors. I have a data frame with two factors: df - expand.grid( fact1=c(y,b,r), fact2=c(far,por,lis,set), year=1991:2000, value= NA) df[,value] - sample(1:50, 120, replace=TRUE) I can make simple xyplot and change the text of the factor levels with strip.custom: require(lattice) xyplot( value ~ year | fact1, data=df, type=b, subset= fact2==far, strip = strip.custom(bg=gray.colors(1,0.95), factor.levels=c(yellow, black, red)), layout=c(1,3) ) But how can I change the text of the factor levels when using both factors as in this plot: xyplot( value ~ year | fact1*fact2, data=df, type=b) (fact2 levels text should change to: c(faro,porto,lisbon,setubal)) I read the help for strip.default and the emails archive, tried with which.given but could not find out how to accomplish this. Many thanks, Rafael Duarte -- Rafael Duarte Marine Resources Department - DRM IPIMAR - National Research Institute for Agriculture and Fisheries Av. Brasília, 1449-006 Lisbon - Portugal Tel:+351 21 302 7000 Fax:+351 21 301 5948 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] setMethod
Hello R-help:
I was hoping someone could help me understand a particular function i came
across in a package:
$.myClass - function( x, name ) {
sym = paste( foo, name, sep = _ )
if( is.loaded(sym) )
.Call(sym,x)
}
I understand the paste, and .Call part, but I'm not sure how this function
would get called? What exactly is $. ? is it a regular expression? Would a
user call this method directly, or is it an internal function that gets
called according to a class/type.
something similar to:
plot.lm()
Thanks in advance !
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] legends in a plot
Hi there, I have the following plot. The circles and the line do not cross in the main plot. But in the legend, the circle and the line cross. I am wondering if there is a way to make the legend look like the plot without crossing. I looked around but did not find a way to do that. Is it doable? plot(x,x^1.5, pch=1, lty=1, type='b') legend(1,25, legend=hello,pch=1, lty=1) Thanks and have a great weekend! __ 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] legends in a plot
If the main thing you want is just to ensure that the legend and plot are consistent it would be easier to change the plot than change the legend: x - 1:10 plot(x,x^1.5, pch=1) lines(x, x^1.5,lty=1) legend(1,25, legend=hello,pch=1, lty=1) On 9/22/06, Bingshan Li [EMAIL PROTECTED] wrote: Hi there, I have the following plot. The circles and the line do not cross in the main plot. But in the legend, the circle and the line cross. I am wondering if there is a way to make the legend look like the plot without crossing. I looked around but did not find a way to do that. Is it doable? plot(x,x^1.5, pch=1, lty=1, type='b') legend(1,25, legend=hello,pch=1, lty=1) Thanks and have a great weekend! __ 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] Statitics Textbook - any recommendation?
Iuri Gavronski schrieb: Other text I am trying to find is multivariate data analysis (EFA, cluster, mult regression, MANOVA, etc.) with examples with R. Hi Iuri, for your second answer I would recommend B. Everitt: An R and S-PLUS Companion to Multivariate Analysis. Springer 2005. isbn 1-85233-882-2. Best Wolfgang -- privat: Wolfgang Lindner, Stieglitzweg 6, D-42799 Leichlingen __ 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] proj4R library will not install
I'm hoping someone can help me. I have downloaded the proj4R.zip and under my version of R (2.3.1) I install the package from local zip file. This worked great. I then type library(proj4R) to load the library and I get the error: Error in library(proj4R) : 'proj4R' is not a valid package -- installed 2.0.0? I have read through the install documentation and have downloaded and unpacked PROJ.4 to c:\proj\ so the bin is located at C:\proj\bin. I then set the environmental variables PATH which now looks like : %SystemRoot%\system32;%SystemRoot%;%SystemRoot%\System32\Wbem;C:\Program Files\Intel\DMIX;C:\Program Files\UltraEdit;C:\proj and I created a new user variable PROJ_LIB to c:\proj\nad. I'm not sure if I am missing anything here but I still get the 2.0.0 error. If you can help me in any way I would truly appreciate it. Thanks in advance, Philip Bermingham __ 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] Double integral
Hi all, I need to solve double integrals with no closed solution. Calling x and y the two variables we have x ~ Normal(y*v,1) and y ~Half-Normal(0,1). In fact, given a joint funcion g(x,y), I need evaluate the integral of this function under that random structure. Could anyone suggest me a package or even a suitable method to solve this problem? Thanks all, Caio - Você quer respostas para suas perguntas? Ou você sabe muito e quer compartilhar seu conhecimento? Experimente o Yahoo! Respostas! [[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] Double integral
Hi all, I need to solve double integrals with no closed solution. Calling x and y the two variables we have x ~ Normal(y*v,1) and y ~Half-Normal(0,1). In fact, given a joint funcion g(x,y), I need evaluate the integral of this function under that random structure. Could anyone suggest me a package or even a suitable method to solve this problem? Thanks all, Caio - [[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] proj4R library will not install
On Fri, 2006-09-22 at 17:16 -0400, Philip Bermingham wrote: I'm hoping someone can help me. I have downloaded the proj4R.zip and under my version of R (2.3.1) I install the package from local zip file. This worked great. I then type library(proj4R) to load the library and I get the error: Error in library(proj4R) : 'proj4R' is not a valid package -- installed 2.0.0? I have read through the install documentation and have downloaded and unpacked PROJ.4 to c: \proj\ so the bin is located at C:\proj\bin. I then set the environmental variables PATH which now looks like : %SystemRoot% \system32;%SystemRoot%;%SystemRoot%\System32\Wbem;C:\Program Files \Intel\DMIX;C:\Program Files\UltraEdit;C:\proj and I created a new user variable PROJ_LIB to c:\proj\nad. I'm not sure if I am missing anything here but I still get the 2.0.0 error. If you can help me in any way I would truly appreciate it. Thanks in advance, Philip Bermingham proj4R is not a base or CRAN R package. Some Googling suggests that the R package might be deprecated, as it has not been updated for some time (hence the error msgs) based upon a review of the R related archive files available at: http://spatial.nhh.no/R/Devel/ I would suggest communicating with Roger Bivand (who I have cc'd here) as to the status of the package and any subsequent replacements. 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.
Re: [R] Double integral
Caio Lucidius Naberezny Azevedo said the following on 9/22/2006 4:40 PM: Hi all, I need to solve double integrals with no closed solution. Calling x and y the two variables we have x ~ Normal(y*v,1) and y ~Half-Normal(0,1). In fact, given a joint funcion g(x,y), I need evaluate the integral of this function under that random structure. Could anyone suggest me a package or even a suitable method to solve this problem? Thanks all, Caio - [[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. Have you tried RSiteSearch(double integral) as the posting guide suggests? --sundar __ 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] setMethod
On 9/22/2006 4:02 PM, Sender wrote:
Hello R-help:
I was hoping someone could help me understand a particular function i came
across in a package:
$.myClass - function( x, name ) {
sym = paste( foo, name, sep = _ )
if( is.loaded(sym) )
.Call(sym,x)
}
I understand the paste, and .Call part, but I'm not sure how this function
would get called? What exactly is $. ? is it a regular expression? Would a
user call this method directly, or is it an internal function that gets
called according to a class/type.
It's the implementation of the $ function for myClass objects, so if x
was a myClass object, it would be called when something like
x$name
was used in an expression. It would call an external function named
foo_name, passing x as an argument.
Duncan Murdoch
something similar to:
plot.lm()
Thanks in advance !
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.
[R] two questions associated with heatmap
hi, there:
i have 2 questions associated with heatmap
in heatmap.2{gplot}, there is a bar called raw z-score showing the
coloring legend for each pixel. Where can I find that z-score's
formulae?
question 2,
for example, I have 5 groups and I want to label each group with a
color name from #FF to #FF evenly. so basically i need a
vector like this:
c(#FF, ?, ?, ?, #FF)
the number of groups can be 10 or whatever.
thanks.
--
Weiwei Shi, Ph.D
Research Scientist
GeneGO, Inc.
Did you always know?
No, I did not. But I believed...
---Matrix III
__
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] Adding .R to source file keeps R from reading it?
So...are you trying to modify a contributed package by adding a *.R file to the 'R' subdirectory in the package? One thing to consider besides the previous tip is that the package might be using a NAMESPACE, which lives in the package root directory (one directory up from the 'R' subdirectory). If so, then you must add the name of your function to the argument list of the call to 'export' in the NAMESPACE file e.g. export(sortgenes, pvalues, MyFunction) -Original Message- From: John Tillinghast [mailto:[EMAIL PROTECTED] Sent: Thu 9/21/2006 9:55 PM To: r-help@stat.math.ethz.ch Subject: Re: [R] Adding .R to source file keeps R from reading it? Yes, this was exactly the problem: I was using the unzipped package, not the source. Now it works! -- Forwarded message -- From: Deepayan Sarkar [EMAIL PROTECTED] Date: Sep 21, 2006 2:58 PM Subject: Re: [R] Adding .R to source file keeps R from reading it? To: John Tillinghast [EMAIL PROTECTED] Cc: r-help@stat.math.ethz.ch On 9/21/06, John Tillinghast [EMAIL PROTECTED] wrote: Hi, I'm updating the LMGene package from Bioconductor. Writing R Extensions suggests that all source files (the ones in the R directory) have a .R ending, so I added it to the (one) source file. The next time I installed and ran R, R didn't understand any of the functions. I tried various things and eventually went back to the file and dropped the .R ending, installed, ran R. It worked! For purposes of distributing the package, do I want to leave the name without the .R, or add the .R and change something else? I'm guessing that the source you are working on has been obtained by unzipping the windows binary zip file. Despite appearances, that is not the source code. For the proper source code, download the file that's marked as source. In this case, http://bioconductor.org/packages/1.8/bioc/html/LMGene.html clearly labels the following as Source (and the corresponding zip file as Windows Binary) http://bioconductor.org/packages/1.8/bioc/src/contrib/LMGene_1.0.0.tar.gz This likely answers your other question as well. -Deepayan [[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] install error uder Cygwin
Dear R users, I have a question about R installation under Cygwin. When running ./configure, I can't pass the checking phase due to an error: --with-readline=yes (default) and headers/libs are not available... Anybody who can tell me why this error arises? Thanks a lot. Xiaohui Chen Dept. of Statistics UBC, Canada _ Buy what you want when you want it on Sympatico / MSN Shopping __ 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] install error uder Cygwin
On 9/22/2006 10:29 PM, X.H Chen wrote: Dear R users, I have a question about R installation under Cygwin. When running ./configure, I can't pass the checking phase due to an error: --with-readline=yes (default) and headers/libs are not available... Anybody who can tell me why this error arises? Thanks a lot. We don't support builds under Cygwin. 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] [BioC] two questions associated with heatmap
Weiwei Shi wrote:
hi, there:
i have 2 questions associated with heatmap
in heatmap.2{gplot}, there is a bar called raw z-score showing the
coloring legend for each pixel. Where can I find that z-score's
formulae?
A Z-score is the mean of the group divided by the standard deviation.
question 2,
for example, I have 5 groups and I want to label each group with a
color name from #FF to #FF evenly. so basically i need a
vector like this:
c(#FF, ?, ?, ?, #FF)
You could look at RColorBrewer or geneplotter for some ideas on how to
do this.
If you don't have a copy of the Bioconductor book, it may be a good idea
to pick one up. It is quite helpful for getting started using
bioconductor and gene expression data.
Hope this helps.
Sean
__
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] [BioC] two questions associated with heatmap
Sean Davis wrote:
Weiwei Shi wrote:
hi, there:
i have 2 questions associated with heatmap
in heatmap.2{gplot}, there is a bar called raw z-score showing the
coloring legend for each pixel. Where can I find that z-score's
formulae?
A Z-score is the mean of the group divided by the standard deviation.
It's a bit late--what I meant to say is that the z-score is the (actual
value minus the mean of the group) divided by the standard deviation.
Sean
__
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] Fatal error: unable to restore saved data in .RData --- no package called 'nlme'
I am using Windows XP and R 2.3.1. During my lastest session I updated my packages and now when I try to start R 2.3.1 I get the following error message: Fatal error: unable to restore saved data in .RData in an error window and Error in loadNamespace(name): there is no package called 'nlme' in the R console. I had been using nlme and lme4 when I updated my packages. I did a search on this error but unfortunatley I don't really understand the solution. I would like to continue using R and 'nlme'. Any suggestions, Mike [[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] contrasts in aov
useRs, A no doubt simple question, but I am baffled. Indeed, I think I once knew the answer, but can't recover it. The default contrasts for aov (and lm, and...) are contr.treatment and contr.poly for unordered and ordered factors, respectively. But, how does one invoke the latter? That is, in a data.frame, how does one indicate that a factor is an *ordered* factor such that contr.poly is invoked in the aov or lm call? -- Please avoid sending me Word or PowerPoint attachments. See http://www.gnu.org/philosophy/no-word-attachments.html -Dr. John R. Vokey __ 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] Y-axis on pbinom
With this: plot(seq(from=0,to=10,by=1),1-pbinom (seq(from=0,to=10,by=1),size=6,prob=0.50),pch=15) How do you change the Y-axis from 0 ~ 0.6? I only get 0.0 ~1.0. thx much __ 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.
