[R] ploting the two sets of data side by side
Hello R-users, I am new to R-commands. I have two sets of data: x - c(7, 7 , 8, 9, 15, 17, 18) y - c(7, 8, 9, 15, 17, 19, 20, 20, 25, 23, 22) I have used 'cut' command to seperate them as follows a - cut(x, breaks =c(0,5,10,20,25,30)) b - cut(y, breaks =c(0,5,10,20,25,30)) table(a) a (0,5] (5,10] (10,20] (20,25] (25,30] 0 4 3 0 0 table(b) b (0,5] (5,10] (10,20] (20,25] (25,30] 0 3 5 3 0 Now if I want to a single graph with both sets of data side by side for same range. Can some one help me regarding the above problem. With Regards Subhabrata Pal [EMAIL PROTECTED] [[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
Re: [R] R newbie...
Thank you for all your answers...
I solved my problem thanks to you all !
david
2005/12/6, paul sorenson [EMAIL PROTECTED]:
Return something that can hold more than one value, eg:
calculate - function(x, y) {
list(a=x+y, b=x-y)
}
David Hajage wrote:
Thank you for your answer.
And what if my first function gives 2 results :
calculate - function(x,y)
{
a - x + y
b - x - y
}
How can I use both a and b in a new function ?
--
David
[[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
Re: [R] ploting the two sets of data side by side
try : barplot(do.call(rbind,lapply(list(x,y), function(x) table(cut(x, breaks =c(0,5,10,20,25,30),beside=T) Subhabrata a écrit : Hello R-users, I am new to R-commands. I have two sets of data: x - c(7, 7 , 8, 9, 15, 17, 18) y - c(7, 8, 9, 15, 17, 19, 20, 20, 25, 23, 22) I have used 'cut' command to seperate them as follows a - cut(x, breaks =c(0,5,10,20,25,30)) b - cut(y, breaks =c(0,5,10,20,25,30)) table(a) a (0,5] (5,10] (10,20] (20,25] (25,30] 0 4 3 0 0 table(b) b (0,5] (5,10] (10,20] (20,25] (25,30] 0 3 5 3 0 Now if I want to a single graph with both sets of data side by side for same range. Can some one help me regarding the above problem. With Regards Subhabrata Pal [EMAIL PROTECTED] [[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 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Bandwidth selection for ksmooth( )
Dear R Users, Before running ksmooth( ), a suitable bandwidth selection is needed. I use some functions for this task and receive these results for my data: width.SJ(y,nb=100,method=ste) : 40.25 bcv(y,nb=100) : 40.53 ucv(y): 41.26 bandwidth.nrd(y) : 45.43 After implementing the function ksmooth(x,y, bandwidth= each of abovementioned bandwidths), I have some NAs in output with regard to each bandwidth. But if the bandwidth be equal to 62 or bigger, then the function ksmooth(x,y, bandwidth=62 or bigger) works without NA,i.e., the bandwidth must be at least 62 in order to not have NA in output. That was the first case. In second case, I have to choose bandwidth bigger than the number of observations in order to not have NA in output of the function of ksmooth(). What is my mistake ? Thank you for any help. Amir Safari - Let fate take it's course directly to your email. [[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
[R] Constructing a transition matrix
If you order your factor levels in your vectors in the order you want in the output, then the prop.table(prop()) command will give you what you want. But you have to reorder the factor levels so that the levels commands give the following output: levels(trans$class) [1] seed veg repr levels(trans$fate) [1] seed veg repr dead That means that you need to include seed as a factor level in fate (even if that level is unused). The prop.table(table()) command will then produce a 3 by 4 table. Remove the last row (that contains the proportion dead), and your set. Hans Gardfjell Dept of Ecology and Environmental science Umeå University, Sweden Chris Stubben wrote: Hi again, I almost figured this out, but still need some help on the last part. I can use prop.table to get survival probabilities... A - t(prop.table( table(trans$class, trans$fate),1) ) rep seed veg dead 0.000 0.333 0.000 rep 0.500 0.000 1.000 veg 0.500 0.667 0.000 so now I just need to format the matrix. I thought I could create a matrix of zeroes using size class names, dev- c(seed,veg, rep). A0-matrix(numeric(9), nrow=3, dimnames=list(dev,dev) ) seed veg rep seed 0 0 0 veg 0 0 0 rep 0 0 0 but how do I assign values in A to the corresponding rows and columns in A0? I hope there is an easy solution that I'm overlooking. seed veg rep seed 0 0 0 veg 0.67 0 0.5 rep 0 1 0.5 Thanks, Chris Chris Stubben wrote: / Hi, // // I would like to construct a transition matrix from a data frame with // annual transitions of marked plants. // // plant-c(1:6) // class-c(seed,seed, seed, veg, rep, rep) // fate-c(dead, veg,veg,rep, rep, veg) // // trans-data.frame(plant, class, fate) // // plant class fate // 1 1 seed dead // 2 2 seed veg // 3 3 seed veg // 4 4 veg rep // 5 5 rep rep // 6 6 rep veg // // I have been using sql queries to do this, but I would like to construct // the matrix in R since I plan to resample transitions using // trans[sample(nrow(trans), 6, replace=T), ] // // I know I can get the original size vector using table() // // data.matrix(table(trans$class)) // [,1] // rep 2 // seed 3 // veg 1 // // // but I don't know how to get counts of each class-fate combination where // fate does NOT equal dead // // seed veg = 2 // veg rep = 1 // rep rep = 1 // rep veg = 1 // // // or how to divide the class-fate count by the original class count in the // size vector to get survival probabilities // // seed veg = 2 / 3 seed = 0.67 // veg rep = 1 / 1 veg = 1 // rep rep = 1 / 2 rep = 0.5 // rep veg = 1 / 2 rep = 0.5 // // // or construct the square matrix with rows and columns in the same // developmental sequence like dev- c(seed,veg, rep). // // seed veg rep // seed 0 0 0 // veg 0.67 0 0.5 // rep 0 1 0.5 // // Any help or suggestions would be appreciated. // Thanks, // // // Chris Stubben // // // -- // Los Alamos National Lab // BioScience Division // MS M888 // Los Alamos, NM 87545 // // / __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] summary[[r.squared]] gives strange results
Felix Flory wrote:
I am simulating an ANOVA model and get a strange behavior from the
summary function. To be more specific: please run the following code
and see for yourself: the summary()[[r.squared]] values of two
identical models are quite different!!
## 3 x 3 ANOVA of two factors x and z on outcome y
s.size - 300 # the sample size
p.z - c(0.25, 0.5, 0.25) # the probabilities of factor z
## probabilities of x given z
p.x.z - matrix(c(40/60, 20/120, 6/60,
14/60, 80/120, 14/60,
6/60, 20/120, 40/60), 3, 3, byrow = TRUE)
## the regression coefficients according to the model.matrix X, that
## is computed later
beta - c(140, -60, -50, -80, 120, 100, -40, 60, 50)/40
## the factor z and the factor x (in dependence of z)
z - x - vector(mode = integer, length = s.size)
for(j in 1:s.size) {
z[j] - sample(1:3, 1, prob = p.z)
x[j] - sample(1:3, 1, prob = p.x.z[, z[j]])
}
## constructing the model.matrix X
zf - factor(z)
contrasts(zf) - contr.treatment(nlevels(zf), base = 3)
zm - model.matrix(~ zf)
xf - factor(x)
contrasts(xf) - contr.treatment(nlevels(xf), base = 3)
xm - model.matrix(~ xf)
X - cbind(zm, zm * xm[,2], zm * xm[,3])
## the outcome y
y - X %*% beta + rnorm(s.size, 0, 4)
## the two linear models
lm.v1 - lm(y ~ X -1)
lm.v2 - lm(y ~ zf * xf)
## which are equivalent
anova(lm.v1, lm.v2)
print(sd(model.matrix(lm.v1) %*% coef(lm.v1))^2 / sd(y)^2) -
print(sd(model.matrix(lm.v2) %*% coef(lm.v2))^2 / sd(y)^2)
## so far everything is fine but why are the following r.squared
## values quite different?
print(summary(lm.v1)[[r.squared]]) -
print(summary(lm.v2)[[r.squared]])
Hi, Felix,
The first model fits your data without an intercept and thus has a
different formula of R^2. The justification should be in any intro
regression text. Here is the relevant snippet from summary.lm:
mss - if (attr(z$terms, intercept))
sum((f - mean(f))^2)
else sum(f^2)
rss - sum(r^2)
snip
ans$r.squared - mss/(mss + rss)
Try the following to see a direct comparison of the two methods:
lm.v1 - lm(y ~ X - 1)
lm.v2 - lm(y ~ zf * xf)
lm.v3 - lm(y ~ X[, -1])
summary(lm.v1)$r.squared
summary(lm.v2)$r.squared
summary(lm.v3)$r.squared
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
Re: [R] heatmap aspect ratio
You can change the code of layout: layout(lmat, widths = lwid, heights = lhei, respect = TRUE) layout(lmat, widths = lwid, heights = lhei, respect = FALSE) (best to create a new function my.layout with the modified code) Jacob Michaelson wrote: Hi all, Does anyone know of a fairly easy way to stretch a heatmap vertically? I've got 42 arrays and would like to be able to see as many significant genes as possible (right now I can only get 50 genes with it still being readable). In some comparisons there are several hundred significant genes. I've fiddled with the asp argument, but that doesn't give the results I'm looking for -- only scales the images, not the dendrograms. Is there any way to make the heatmap rectangular rather than square without hacking the heatmap function itself (which is where I'm headed next)? Thanks in advance, Jake __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Dots argument in apply method
Hello everyone,
I'm working on a package using S4 classes and methods and I ran into the
following problem when I tried to create an apply method for objects
of one of my new classes. I've found a way around the problem but I
wonder if I did not paint myself into the corner. I'd like your opinion
about that.
So I have an object myObj of class myClass. I define a new function
.apply.myClass which is a myClass specific version of apply. The
trick is that I would like to have an additional formal argument in
.apply.myClass compared to apply. More precisely we have:
apply(X, MARGIN, FUN, ...)
and I want:
.apply.myClass(x, margin, fun, groups = NULL, ...)
As long as I stay at the function level there is no problem. Life
becomes harder when I want to define an apply method for myClass
objects, method which should call .apply.myClass.
The formal argument groups in the myClass specific apply method will
have to be passed in the dots argument, together with the FUN specific
arguments. Then if the groups argument is provided it will have to be
extracted and the remaining dots argument(s), if any, will have to be
passed as such to .apply.myClass. Here is the way I did it:
## Start by setting a generic apply method
if (!isGeneric(apply))
setGeneric(apply, function(X, MARGIN, FUN, ...)
standardGeneric(apply))
## set apply method for myClass objects
setMethod(apply,
signature(X = myClass,
MARGIN = numeric,
FUN = function),
function(X, MARGIN, FUN, ...) {
.call - match.call(.apply.myClass)
if (is.null(.call$groups)) myGroups - NULL
else myGroups - .call$groups
argList - list(obj = .call$obj,
margin = .call$margin,
fun = .call$fun,
groups = myGroups
)
if(!all(names(.call)[-1] %in% names(formals(.apply.myClass {
## Some dots arguments have been provided
otherNames - (names(.call)[-1])[!(names(.call)[-1] %in%
names(formals(.apply.myClass)))]
remainingDots - lapply(otherNames, function(i) .call[[i]])
names(remainingDots) - otherNames
argList - c(argList,remainingDots)
}
do.call(.apply.myClass, args = argList)
}
)
Does anyone have a quicker solution?
Thanks in advance,
Christophe.
PS: If you want a full example with actual class and .apply.myClass
definitions, here is one:
## define class myClass
setClass(myClass, representation(Data = data.frame, timeRange =
numeric))
## create myObj an instantiation of myClass
myObj - new(myClass,
Data = data.frame(Time = sort(runif(10)),
observation = I(matrix(rnorm(20),nrow=10,ncol=2)),
label = factor(rep(1:2,5),levels = 1:2, labels = c(cat.
1, cat. 2))
),
timeRange = c(0,1)
)
## create function .apply.myClass for myClass objects
.apply.myClass - function(obj, ## object of class myClass
margin, ## a numeric which should be 1 or 2
fun, ## a function
groups = NULL, ## should fun be applied in a
group
## specific manner?
... ## additional arguments passed to fun
) {
## attach the data frame contained in obj
attach([EMAIL PROTECTED])
## make sure to detach it at the end
on.exit(detach([EMAIL PROTECTED]))
## get the variable names
variableNames - names([EMAIL PROTECTED])
## check that one variable is named observation
if (!(observation %in% variableNames))
stop(paste(The slot Data of,
deparse(substitute(obj)),
does not contain an observation variable as it should.
)
)
if (margin == 1) {
## in that case we don't care of the group
myResult - apply(observation, 1, fun, ...)
return(myResult)
} else if (margin == 2) {
if (is.null(groups)) {
## no groups defined
myResult - apply(observation, 2, fun, ...)
return(myResult)
} else {
## groups defined
groups - eval(groups)
X - levels(groups)
dim(X) - c(1,length(X))
myResult - apply(X,
2,
function(i) apply(observation[groups == i,],
2,
fun, ...)
)
return(myResult)
}
} else {
stop(margin should be set to 1 or 2.)
}
}
--
A Master Carpenter has many tools and is expert with most of them.If you
only know how to use a hammer, every problem starts to look like a nail.
Stay away from that trap.
Richard B Johnson.
--
Christophe Pouzat
Laboratoire de Physiologie Cerebrale
CNRS UMR
Re: [R] Plot
On Wed, 2005-12-07 at 18:08 +1100, paul sorenson wrote: Apple Ho wrote: Hello, I have a problem about using the command plot. Suppose I have some points, and one of them is (0,0), how can I show the figure with this point which is at the corner? How close to the corner do you want it? plot(0, 0, xlim=c(0, 1), ylim=c(0,1)) you could also add: grid() By default, R adds +/- 4% to each axis range, based upon the range of the x and y values or the 'xlim' and 'ylim' arguments. This behavior is set by pars 'xaxs' and 'yaxs', which are set to r by default. Setting both to i will provide exact axis ranges. Note the difference between: plot(0:5, 0:5) and plot(0:5, 0:5, xaxs = i, yaxs = i) See ?par for more information. 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
Re: [R] ploting the two sets of data side by side
Or building on that solution but eliminating the do.call and lapply: f - function(x) table(cut(x, breaks = seq(0, 30, 5))) barplot(rbind(f(x), f(y)), beside = TRUE) On 12/7/05, Jacques VESLOT [EMAIL PROTECTED] wrote: try : barplot(do.call(rbind,lapply(list(x,y), function(x) table(cut(x, breaks =c(0,5,10,20,25,30),beside=T) Subhabrata a écrit : Hello R-users, I am new to R-commands. I have two sets of data: x - c(7, 7 , 8, 9, 15, 17, 18) y - c(7, 8, 9, 15, 17, 19, 20, 20, 25, 23, 22) I have used 'cut' command to seperate them as follows a - cut(x, breaks =c(0,5,10,20,25,30)) b - cut(y, breaks =c(0,5,10,20,25,30)) table(a) a (0,5] (5,10] (10,20] (20,25] (25,30] 0 4 3 0 0 table(b) b (0,5] (5,10] (10,20] (20,25] (25,30] 0 3 5 3 0 Now if I want to a single graph with both sets of data side by side for same range. Can some one help me regarding the above problem. With Regards Subhabrata Pal [EMAIL PROTECTED] [[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 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Dots argument in apply method
Why does simply
setMethod(apply,
signature(X = myClass,
MARGIN = numeric,
FUN = function),
function(X, MARGIN, FUN, ...) .apply.myClass(X, MARGIN, FUN, ...))
not do what you want? It works for me in your example, e.g.
apply(myObj, 2, sum, [EMAIL PROTECTED])
gives exactly the same answer as your complicated solution.
I do wonder if you have misunderstood what '...' does.
I also wonder why you chose to overload a basic R function as an S4
generic like this. If you think that thereby existing calls to apply()
will go via your S4 methods then I fear you have overlooked the effects of
namespaces.
A simpler example
setClass(myClass, representation(tt=numeric))
setMethod(lapply, signature(X=myClass), function(X, FUN, ...) FUN([EMAIL
PROTECTED]))
myObj - new(myClass, tt=1:10)
lapply(myObj, sum)
[1] 55
sapply(myObj, sum)
list()
since sapply is calling base::lapply, not the lapply S4 generic.
On Wed, 7 Dec 2005, Christophe Pouzat wrote:
Hello everyone,
I'm working on a package using S4 classes and methods and I ran into the
following problem when I tried to create an apply method for objects
of one of my new classes. I've found a way around the problem but I
wonder if I did not paint myself into the corner. I'd like your opinion
about that.
So I have an object myObj of class myClass. I define a new function
.apply.myClass which is a myClass specific version of apply. The
trick is that I would like to have an additional formal argument in
.apply.myClass compared to apply. More precisely we have:
apply(X, MARGIN, FUN, ...)
and I want:
.apply.myClass(x, margin, fun, groups = NULL, ...)
As long as I stay at the function level there is no problem. Life
becomes harder when I want to define an apply method for myClass
objects, method which should call .apply.myClass.
The formal argument groups in the myClass specific apply method will
have to be passed in the dots argument, together with the FUN specific
arguments. Then if the groups argument is provided it will have to be
extracted and the remaining dots argument(s), if any, will have to be
passed as such to .apply.myClass. Here is the way I did it:
## Start by setting a generic apply method
if (!isGeneric(apply))
setGeneric(apply, function(X, MARGIN, FUN, ...)
standardGeneric(apply))
## set apply method for myClass objects
setMethod(apply,
signature(X = myClass,
MARGIN = numeric,
FUN = function),
function(X, MARGIN, FUN, ...) {
.call - match.call(.apply.myClass)
if (is.null(.call$groups)) myGroups - NULL
else myGroups - .call$groups
argList - list(obj = .call$obj,
margin = .call$margin,
fun = .call$fun,
groups = myGroups
)
if(!all(names(.call)[-1] %in% names(formals(.apply.myClass {
## Some dots arguments have been provided
otherNames - (names(.call)[-1])[!(names(.call)[-1] %in%
names(formals(.apply.myClass)))]
remainingDots - lapply(otherNames, function(i) .call[[i]])
names(remainingDots) - otherNames
argList - c(argList,remainingDots)
}
do.call(.apply.myClass, args = argList)
}
)
Does anyone have a quicker solution?
Thanks in advance,
Christophe.
PS: If you want a full example with actual class and .apply.myClass
definitions, here is one:
## define class myClass
setClass(myClass, representation(Data = data.frame, timeRange =
numeric))
## create myObj an instantiation of myClass
myObj - new(myClass,
Data = data.frame(Time = sort(runif(10)),
observation = I(matrix(rnorm(20),nrow=10,ncol=2)),
label = factor(rep(1:2,5),levels = 1:2, labels = c(cat.
1, cat. 2))
),
timeRange = c(0,1)
)
## create function .apply.myClass for myClass objects
.apply.myClass - function(obj, ## object of class myClass
margin, ## a numeric which should be 1 or 2
fun, ## a function
groups = NULL, ## should fun be applied in a
group
## specific manner?
... ## additional arguments passed to fun
) {
## attach the data frame contained in obj
attach([EMAIL PROTECTED])
## make sure to detach it at the end
on.exit(detach([EMAIL PROTECTED]))
## get the variable names
variableNames - names([EMAIL PROTECTED])
## check that one variable is named observation
if (!(observation %in% variableNames))
stop(paste(The slot Data of,
deparse(substitute(obj)),
does not contain an observation variable as it
Re: [R] Matrix of dummy variables from a factor
Bert, how about when making predictions of contrasts using, for example, estimable() from the gmodels package? I find it very useful. Andrew -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Berton Gunter Sent: Wednesday, 7 December 2005 8:27 AM To: 'Liaw, Andy'; 'Charles H. Franklin'; r-help@stat.math.ethz.ch Subject: Re: [R] Matrix of dummy variables from a factor But note: There are (almost?) no situations in R where the dummy variables coding is needed. The coding is (almost?) always handled properly by the modeling functions themselves. Question: Can someone provide a straightforward example where the dummy variable coding **is** explicitly needed? -- 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 -- Andrew Robinson Senior Lecturer in Statistics Tel: +61-3-8344-9763 Department of Mathematics and StatisticsFax: +61-3-8344-4599 University of Melbourne, VIC 3010 Australia Email: [EMAIL PROTECTED]Website: http://www.ms.unimelb.edu.au __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R is GNU S, not C.... [was how to get or store .....]
Adaikalavan Ramasamy wrote: Yes, it drives me mad too when people use = instead of - for assignment and suppress spaces in an naive attempt for saving space. In fact, I like the - assignment operator, but tend to write code densely myself as that is the way I like to view it. As R formats the source code it displays in the typographic density approved by the gurus, it does seem a tidge authoritarian to insist that everyone adhere to the same coding style in private. After all, I don't whinge about having to scroll up and down like a yoyo when attempting to decipher other people's spaced out code. I intentionally don't use semicolons as line terminators in R to forestall my C reflexes. The greatest difficulty I have with this sort of interference is between R and Tcl-Tk which, combined with the correct titre of fatigue, almost always leads to syntax errors. I did enjoy Ted's defense of the semicolon in R, as it emphasizes the diversity of styles upon which cooperative programming depends. Jim __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] lmer and glmmPQL
On 12/5/05, Cox, Stephen [EMAIL PROTECTED] wrote: I have been looking into both of these approaches to conducting a GLMM, and want to make sure I understand model specification in each. In particular - after looking at Bates' Rnews article and searching through the help archives, I am unclear on the specification of nested factors in lmer. Do the following statements specify the same mode within each approach? m1 = glmmPQL(RICH ~ ZONE, family = poisson, data, random = ~ YEAR | SITE / QUADRAT) m2 = lmer(RICH ~ ZONE +(YEAR|SITE)+ (YEAR|QUADRAT), family = poisson, data) If you want to ensure that QUADRAT is nested within SITE then use the interaction operator explicitly m2 - lmer(RICH ~ ZONE +(YEAR|SITE)+ (YEAR|SITE:QUADRAT), family = poisson, data) For the grouping factors nested versus non-nested depends on the coding. If QUADRAT has a distinct level for each SITE:QUADRAT combination then the nesting will automatically be detected. However, if the nesting is implicit (that is, if levels of QUADRAT are repeated at different SITES) then it is necessary to use the interaction operator. There is no harm in using the interaction operator when the nesting is explicit. As a follow up - what would be the most appropriate model formula (using glmmPQL syntax) to specify both a nested facor and repeated observations? Specifically, I am dealing with experimental data with three factors. ZONE is a fixed effect. Three sites (SITE) are nested within each ZONE. Multiple quadrats within each SITE are measured across multiple years. I want to represent the nesting of SITE within ZONE and allow for repeated observations within each QUADRAT over time (the YEAR | QUADRAT random effect). -- I am assuming that glmmPQL is the best option at this point because of recent discussion on Rhelp about issues associated with the Matrix package used in lmer (i.e., the anova results do not seem to match parameter tests). I believe the anova problems only occur with a binomial response. They are caused by my failure to use the prior.weights appropriately. For a Poisson model this should not be a problem. Any information would be very much appreciated! Regards Stephen __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] How to simplify
Dear list,
I have a list containing parameters (time and X1), and have n
similar data set like
the following:
cal
[[1]]
timeX1
1 0.0 10.006306
2 0.5 9.433443
3 1.0 8.893405
4 2.0 7.904274
5 4.0 6.243807
6 6.0 4.932158
7 8.0 3.896049
8 10.0 3.077604
[[2]]
timeX1
1 0.0 10.015972
2 0.5 9.460064
3 1.0 8.935039
4 2.0 7.970755
5 4.0 6.343151
6 6.0 5.047900
7 8.0 4.017131
8 10.0 3.196856
[[3]]
time X1
1 0.0 9.985741
2 0.5 9.552583
3 1.0 9.138239
4 2.0 8.362664
5 4.0 7.003394
6 6.0 5.865057
7 8.0 4.911747
8 10.0 4.113382
[[4]]
...
[[n]]
...
And I would like to put all X1( when time=0) together, time=0.5,1...
are the same.
then calculate the mean value.
a-list()
b-list()
c-list()
d-list()
e-list()
...
for(i in 1:n){
+ a[[i]]-cal[[i]][1,2]
+ b[[i]]-cal[[i]][2,2]
+ c[[i]]-cal[[i]][3,2]
+ d[[i]]-cal[[i]][4,2]
+ e[[i]]-cal[[i]][5,2]
+ .
}
mean.a-(a[[1]][1]+a[[2]][1]+a[[3]][1]+.)/n
mean.b-(b[[1]][1]+b[[2]][1]+b[[3]][1]+.)/n
mean.c-(c[[1]][1]+c[[2]][1]+c[[3]][1]+.)/n
mean.d-(d[[1]][1]+d[[2]][1]+d[[3]][1]+.)/n
.
xy-c(mean.a,mean.b,mean.c,mean.d,)
But the way I use seem not very smart.
So please give me some hints to the simplify this.
Thanks in advance !!
Sincerely!!
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R is GNU S, not C.... [was how to get or store .....]
Oh Patrick, surely German Capitalization is better! :) Peter Patrick Connolly wrote: On Tue, 06-Dec-2005 at 04:21PM +, Patrick Burns wrote: | I don't put in extraneous ';' because I maybe get a | blister on my little finger. | | I suspect that those who find the semi-colons ugly in | R do not find them ugly in C. Nor in perl, mysql or php. I quite like how R is rather different from those other languages. Somehow that very different look helps to change (me at least) into the appropriate type of thinking. It's somewhat akin to the difference between English and German use of capital letters -- each system makes perfect sense in its own context and are beter not mixed. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] ploting the two sets of data side by side
As usual, Gabor provides an elegant solution. But I hope that, in this case, the OP provided a toy example. Otherwise, I don't see the point of applying cut() to a vector of length 7. Why not just use stripchart()? Peter Ehlers Gabor Grothendieck wrote: Or building on that solution but eliminating the do.call and lapply: f - function(x) table(cut(x, breaks = seq(0, 30, 5))) barplot(rbind(f(x), f(y)), beside = TRUE) On 12/7/05, Jacques VESLOT [EMAIL PROTECTED] wrote: try : barplot(do.call(rbind,lapply(list(x,y), function(x) table(cut(x, breaks =c(0,5,10,20,25,30),beside=T) Subhabrata a écrit : Hello R-users, I am new to R-commands. I have two sets of data: x - c(7, 7 , 8, 9, 15, 17, 18) y - c(7, 8, 9, 15, 17, 19, 20, 20, 25, 23, 22) I have used 'cut' command to seperate them as follows a - cut(x, breaks =c(0,5,10,20,25,30)) b - cut(y, breaks =c(0,5,10,20,25,30)) table(a) a (0,5] (5,10] (10,20] (20,25] (25,30] 0 4 3 0 0 table(b) b (0,5] (5,10] (10,20] (20,25] (25,30] 0 3 5 3 0 Now if I want to a single graph with both sets of data side by side for same range. Can some one help me regarding the above problem. With Regards Subhabrata Pal [EMAIL PROTECTED] [[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 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Dots argument in apply method
Dear Prof. Ripley, Thanks for the clarification. Prof Brian Ripley wrote: Why does simply setMethod(apply, signature(X = myClass, MARGIN = numeric, FUN = function), function(X, MARGIN, FUN, ...) .apply.myClass(X, MARGIN, FUN, ...)) not do what you want? It works for me in your example, e.g. apply(myObj, 2, sum, [EMAIL PROTECTED]) gives exactly the same answer as your complicated solution. I do wonder if you have misunderstood what '...' does. I hope not (entirely). I got a bizarre bug yesterday while working on that method/function issue which lead me (wrongly, because of a too quick diagnostic probably) to think there was a problem with the way the ... was passed from one function to the next. That explains the complicated construct of my previous e-mail. Your solution works (of course) and is definitely simpler (!). I also wonder why you chose to overload a basic R function as an S4 generic like this. If you think that thereby existing calls to apply() will go via your S4 methods then I fear you have overlooked the effects of namespaces. No, I just want the user to be able to call apply with myClass objects without having to bother with the internal structure of these objects (which is slightly more complicated in my real case than on the example I sent). Thanks again, Christophe. -- A Master Carpenter has many tools and is expert with most of them.If you only know how to use a hammer, every problem starts to look like a nail. Stay away from that trap. Richard B Johnson. -- Christophe Pouzat Laboratoire de Physiologie Cerebrale CNRS UMR 8118 UFR biomedicale de l'Universite Paris V 45, rue des Saints Peres 75006 PARIS France tel: +33 (0)1 42 86 38 28 fax: +33 (0)1 42 86 38 30 web: www.biomedicale.univ-paris5.fr/physcerv/C_Pouzat.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to simplify
I suggest that when displaying test data in a post that
you do it like this:
dput(cal)
since then others can simply copy and paste it into
their session.
At any rate, using this test data:
cal - list(A = data.frame(time = 1:3, X1 = 1:3),
B = data.frame(1:3, X1 = 3:5))
Pick off the second element of each list component, turn
the result into a data frame and take the row means:
rowMeans(as.data.frame(lapply(cal, [, 2)))
Since the times are irregular, you might prefer to represent cal
as a multivariate zoo object using the zoo library:
library(zoo)
cal.zoo - do.call(merge, lapply(cal, function(x) zoo(x[,2], x[,1])))
Once you have done this and other operations simplify. For
this one its just:
rowMeans(cal.zoo)
or to plot them all
plot(cal.zoo) # separate plots
plot(call.zoo, plot.type = single) # all on one plot
On 12/7/05, Rhett Eckstein [EMAIL PROTECTED] wrote:
Dear list,
I have a list containing parameters (time and X1), and have n
similar data set like
the following:
cal
[[1]]
timeX1
1 0.0 10.006306
2 0.5 9.433443
3 1.0 8.893405
4 2.0 7.904274
5 4.0 6.243807
6 6.0 4.932158
7 8.0 3.896049
8 10.0 3.077604
[[2]]
timeX1
1 0.0 10.015972
2 0.5 9.460064
3 1.0 8.935039
4 2.0 7.970755
5 4.0 6.343151
6 6.0 5.047900
7 8.0 4.017131
8 10.0 3.196856
[[3]]
time X1
1 0.0 9.985741
2 0.5 9.552583
3 1.0 9.138239
4 2.0 8.362664
5 4.0 7.003394
6 6.0 5.865057
7 8.0 4.911747
8 10.0 4.113382
[[4]]
...
[[n]]
...
And I would like to put all X1( when time=0) together, time=0.5,1...
are the same.
then calculate the mean value.
a-list()
b-list()
c-list()
d-list()
e-list()
...
for(i in 1:n){
+ a[[i]]-cal[[i]][1,2]
+ b[[i]]-cal[[i]][2,2]
+ c[[i]]-cal[[i]][3,2]
+ d[[i]]-cal[[i]][4,2]
+ e[[i]]-cal[[i]][5,2]
+ .
}
mean.a-(a[[1]][1]+a[[2]][1]+a[[3]][1]+.)/n
mean.b-(b[[1]][1]+b[[2]][1]+b[[3]][1]+.)/n
mean.c-(c[[1]][1]+c[[2]][1]+c[[3]][1]+.)/n
mean.d-(d[[1]][1]+d[[2]][1]+d[[3]][1]+.)/n
.
xy-c(mean.a,mean.b,mean.c,mean.d,)
But the way I use seem not very smart.
So please give me some hints to the simplify this.
Thanks in advance !!
Sincerely!!
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to simplify
assuming that you have the same number of measurements in each
sub-data.frame, you could use something like:
cal - lapply(1:10, function(x) data.frame(time = c(0, 0.5, 1, 2, 4,
6, 8, 10), X1 = rnorm(8, 10:3)))
##
rowMeans(as.data.frame(lapply(cal, [, X1)))
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://www.med.kuleuven.be/biostat/
http://www.student.kuleuven.be/~m0390867/dimitris.htm
- Original Message -
From: Rhett Eckstein [EMAIL PROTECTED]
To: r-help@stat.math.ethz.ch
Sent: Wednesday, December 07, 2005 3:38 PM
Subject: [R] How to simplify
Dear list,
I have a list containing parameters (time and X1), and have n
similar data set like
the following:
cal
[[1]]
timeX1
1 0.0 10.006306
2 0.5 9.433443
3 1.0 8.893405
4 2.0 7.904274
5 4.0 6.243807
6 6.0 4.932158
7 8.0 3.896049
8 10.0 3.077604
[[2]]
timeX1
1 0.0 10.015972
2 0.5 9.460064
3 1.0 8.935039
4 2.0 7.970755
5 4.0 6.343151
6 6.0 5.047900
7 8.0 4.017131
8 10.0 3.196856
[[3]]
time X1
1 0.0 9.985741
2 0.5 9.552583
3 1.0 9.138239
4 2.0 8.362664
5 4.0 7.003394
6 6.0 5.865057
7 8.0 4.911747
8 10.0 4.113382
[[4]]
...
[[n]]
...
And I would like to put all X1( when time=0) together,
time=0.5,1...
are the same.
then calculate the mean value.
a-list()
b-list()
c-list()
d-list()
e-list()
...
for(i in 1:n){
+ a[[i]]-cal[[i]][1,2]
+ b[[i]]-cal[[i]][2,2]
+ c[[i]]-cal[[i]][3,2]
+ d[[i]]-cal[[i]][4,2]
+ e[[i]]-cal[[i]][5,2]
+ .
}
mean.a-(a[[1]][1]+a[[2]][1]+a[[3]][1]+.)/n
mean.b-(b[[1]][1]+b[[2]][1]+b[[3]][1]+.)/n
mean.c-(c[[1]][1]+c[[2]][1]+c[[3]][1]+.)/n
mean.d-(d[[1]][1]+d[[2]][1]+d[[3]][1]+.)/n
.
xy-c(mean.a,mean.b,mean.c,mean.d,)
But the way I use seem not very smart.
So please give me some hints to the simplify this.
Thanks in advance !!
Sincerely!!
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
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
Re: [R] lmer and glmmPQL
Thanks for the reply Doug! A follow up question and comment ... 1) If I understand correctly, looking at a simple situation in which SITES are nested in ZONES, the following should be similar. However, despite the same F values, the p-value from lmer is 1/2 the other methods. Why is this true? anova(lmer(RICH ~ ZONE + (1|SITE:ZONE), data)) Analysis of Variance Table Df Sum Sq Mean Sq Denom F value Pr(F) ZONE 4 97.824.5 6610.0 2.1046 0.07753 . # make the nesting explicit data$SinZ = with(data, ZONE:SITE)[drop=TRUE] anova(lme(RICH ~ ZONE, data, random = ~1 | SinZ)) numDF denDF F-value p-value (Intercept) 1 6600 100.38331 .0001 ZONE410 2.10459 0.155 summary(aov(RICH ~ ZONE + Error(SITE:ZONE), data)) Error: SITE:ZONE Df Sum Sq Mean Sq F value Pr(F) ZONE 4 296697417 2.1046 0.155 Residuals 10 352433524 2) I think the anova problems with lmer may also apply to poisson. Compare the following (which includes a covariate). Based on the parameter estimates, the covariate should be significant. However, its anova p-value is .998: lmer(RICH ~ ZONE + lANPP + (1|SITE:ZONE), family = poisson, data) Generalized linear mixed model fit using PQL Formula: RICH ~ ZONE + lANPP + (1 | SITE:ZONE) Data: data Family: poisson(log link) AIC BIClogLik deviance 9700.252 9754.628 -4842.126 9684.252 Random effects: GroupsNameVarianceStd.Dev. SITE:ZONE (Intercept)0.069493 0.26361 # of obs: 6615, groups: SITE:ZONE, 15 Estimated scale (compare to 1) 1.183970 Fixed effects: Estimate Std. Error z value Pr(|z|) (Intercept) 1.5169605 0.1533564 9.8917 2e-16 *** ZONE20.4034169 0.2156956 1.8703 0.06144 . ZONE3 -0.1772011 0.2158736 -0.8209 0.41173 ZONE4 -0.2368290 0.2158431 -1.0972 0.27254 ZONE5 -0.1011186 0.2158114 -0.4686 0.63939 lANPP0.2201926 0.0081857 26.8995 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 anova(lmer(RICH ~ ZONE + lANPP + (1|SITE:ZONE), family = poisson, data)) Analysis of Variance Table DfSum Sq Mean Sq Denom F value Pr(F) ZONE 4 2.809e-05 7.022e-06 6609 4.298e-06 1. lANPP 1 5.229e-06 5.229e-06 6609 3.200e-06 0.9986 Thanks again for any insight you may be able to provide!! -Original Message- From: Douglas Bates [mailto:[EMAIL PROTECTED] Sent: Wednesday, December 07, 2005 8:28 AM To: Cox, Stephen Cc: r-help@stat.math.ethz.ch Subject: Re: [R] lmer and glmmPQL On 12/5/05, Cox, Stephen [EMAIL PROTECTED] wrote: I have been looking into both of these approaches to conducting a GLMM, and want to make sure I understand model specification in each. In particular - after looking at Bates' Rnews article and searching through the help archives, I am unclear on the specification of nested factors in lmer. Do the following statements specify the same mode within each approach? m1 = glmmPQL(RICH ~ ZONE, family = poisson, data, random = ~ YEAR | SITE / QUADRAT) m2 = lmer(RICH ~ ZONE +(YEAR|SITE)+ (YEAR|QUADRAT), family = poisson, data) If you want to ensure that QUADRAT is nested within SITE then use the interaction operator explicitly m2 - lmer(RICH ~ ZONE +(YEAR|SITE)+ (YEAR|SITE:QUADRAT), family = poisson, data) For the grouping factors nested versus non-nested depends on the coding. If QUADRAT has a distinct level for each SITE:QUADRAT combination then the nesting will automatically be detected. However, if the nesting is implicit (that is, if levels of QUADRAT are repeated at different SITES) then it is necessary to use the interaction operator. There is no harm in using the interaction operator when the nesting is explicit. As a follow up - what would be the most appropriate model formula (using glmmPQL syntax) to specify both a nested facor and repeated observations? Specifically, I am dealing with experimental data with three factors. ZONE is a fixed effect. Three sites (SITE) are nested within each ZONE. Multiple quadrats within each SITE are measured across multiple years. I want to represent the nesting of SITE within ZONE and allow for repeated observations within each QUADRAT over time (the YEAR | QUADRAT random effect). -- I am assuming that glmmPQL is the best option at this point because of recent discussion on Rhelp about issues associated with the Matrix package used in lmer (i.e., the anova results do not seem to match parameter tests). I believe the anova problems only occur with a binomial response. They are caused by my failure to use the prior.weights appropriately. For a Poisson model this should not be a problem. Any information would be very much appreciated! Regards Stephen __ R-help@stat.math.ethz.ch mailing
[R] Change labels of x-axes in Plot of stl() function?
Hi all, How can the label of the x-axes in the plot() of a stl.object be adapted? e.g., When plotting: plot(stl(nottem, per)) In the labels of the x-axes is time. How can this be changed to e.g., Time (dekade) ? It does not work with xlab or others anymore Thanks, Jan ___ Ir. Jan Verbesselt Research Associate Biosystems Department ~ M³-BIORES Vital Decosterstraat 102, 3000 Leuven, Belgium Tel: +32-16-329750 Fax: +32-16-329760 http://gloveg.kuleuven.ac.be/ ___ Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm [[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
Re: [R] How to simplify
Hi
changing list to matrix and making summary could do the trick
lll - list(a=cbind(1:10, rnorm(10)), b=cbind(1:10, rnorm(10)))
mat - do.call(rbind, lll)
tapply(mat[,2], mat[,1], mean)
BTW I found a suitable thread with similar question in CRAN search
list summary mean
HTH
Petr
On 7 Dec 2005 at 22:38, Rhett Eckstein wrote:
Date sent: Wed, 7 Dec 2005 22:38:56 +0800
From: Rhett Eckstein [EMAIL PROTECTED]
To: r-help@stat.math.ethz.ch
Subject:[R] How to simplify
Dear list,
I have a list containing parameters (time and X1), and have n
similar data set like
the following:
cal
[[1]]
timeX1
1 0.0 10.006306
2 0.5 9.433443
3 1.0 8.893405
4 2.0 7.904274
5 4.0 6.243807
6 6.0 4.932158
7 8.0 3.896049
8 10.0 3.077604
[[2]]
timeX1
1 0.0 10.015972
2 0.5 9.460064
3 1.0 8.935039
4 2.0 7.970755
5 4.0 6.343151
6 6.0 5.047900
7 8.0 4.017131
8 10.0 3.196856
[[3]]
time X1
1 0.0 9.985741
2 0.5 9.552583
3 1.0 9.138239
4 2.0 8.362664
5 4.0 7.003394
6 6.0 5.865057
7 8.0 4.911747
8 10.0 4.113382
[[4]]
...
[[n]]
...
And I would like to put all X1( when time=0) together, time=0.5,1...
are the same. then calculate the mean value. a-list() b-list()
c-list() d-list() e-list() ... for(i in 1:n){ +
a[[i]]-cal[[i]][1,2] + b[[i]]-cal[[i]][2,2] +
c[[i]]-cal[[i]][3,2] + d[[i]]-cal[[i]][4,2] +
e[[i]]-cal[[i]][5,2] + . }
mean.a-(a[[1]][1]+a[[2]][1]+a[[3]][1]+.)/n
mean.b-(b[[1]][1]+b[[2]][1]+b[[3]][1]+.)/n
mean.c-(c[[1]][1]+c[[2]][1]+c[[3]][1]+.)/n
mean.d-(d[[1]][1]+d[[2]][1]+d[[3]][1]+.)/n .
xy-c(mean.a,mean.b,mean.c,mean.d,) But the way I use seem
not very smart. So please give me some hints to the simplify this.
Thanks in advance !! Sincerely!!
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
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
Re: [R] Change labels of x-axes in Plot of stl() function?
If you look through the output of:
stats:::plot.stl
you see right near the end that time is hard coded in the call to mtext.
However, we could temporarily redefine mtext so that it works as you
wish and then redefine plot.stl so that it looks within the environment
of our function to find mtext (rather than looking in the stats package):
plot.stl - function(..., xlab = time) {
mtext - function(text, ...) graphics::mtext(xlab, ...)
plot.stl - stats:::plot.stl
environment(plot.stl) - environment()
plot.stl(...)
}
# test it
example(stl)
plot.stl(stmd, xlab = X)
On 12/7/05, Jan Verbesselt [EMAIL PROTECTED] wrote:
Hi all,
How can the label of the x-axes in the plot() of a stl.object be adapted?
e.g.,
When plotting: plot(stl(nottem, per))
In the labels of the x-axes is time. How can this be changed to e.g.,
Time (dekade) ?
It does not work with xlab or others anymore…
Thanks,
Jan
___
Ir. Jan Verbesselt
Research Associate
Biosystems Department ~ M³-BIORES
Vital Decosterstraat 102, 3000 Leuven, Belgium
Tel: +32-16-329750 Fax: +32-16-329760
http://gloveg.kuleuven.ac.be/
___
Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
[[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
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Closed form for regression splines - solution
Greetings,
The question was about getting closed form equations from a bs() representation
of a curve so enabling publication of a fitted curve for uise outside R. In
case anyone else is interested, the solution lies in exploiting the fact that
we can get accurate derivatives at the breakpoints and from these reconstruct
the individual polynomial pieces.
Prototype function below - maybe the basis for a useful addition to the splines
library?
Steve.
###
require(splines)
# pieces
#
# Prototype function to compute the coefficients of the individual
# polynomial pieces making up a cubic spline
# Works for bs() and ns() with intercept=FALSE
#
# Based on the observation that the derivatives at the breakpoints can be
# reliably computed using splineDesign and these define the functions.
# See De Boor SIAM J Num Anal 14 441-472 (1977)
#
# Arguments:
# call: spline function call used in modelling function (eg bs(x,4)) could
#be taken from attr(fit$terms,predvars)
# data: frame to evaluate call in - usually data argument of modelling
#function call
# coef: spline function coefficients (no intercept assumed!)
#
# Notes:
# Currently assumes intercept=FALSE used in bs and ns call - the usual for #
modelling use
# - relatively trivial to modify to allow this, but messy.
#
# Future:
# Needs bulletproofing and deeper testing and code could be more elegant
# Would like to add a text representation or maybe an R function.
#
# Value:
# A list with components
# npieces: no of pieces
# degree: degree of fitted polynomial pieces
# cutpoints: the ranges of the pieces:
# piece i covers cutpoint[i] to cutpoint[i+1]
# pars: A list with one vector per piece giving coefficients of each
# piece for 1,(x-c),(x-c)^2
# where c=cutpoint[i] for the ith piece.
#
# Steve Roberts Dec 2005
# [EMAIL PROTECTED]
pieces - function(call,coef,data=NULL)
{
basis - eval(call,list(data,sys.parent))
degree - attr(basis,degree)
#knot positions as in bs and ns
if (class(basis)[1]==bs)
Aknots - sort(c(rep(attr(basis,Boundary.knots), degree+1),
attr(basis,knots)))
else if (class(basis)[1]==ns)
Aknots - sort(c(rep(attr(basis,Boundary.knots), 4),
attr(basis,knots)))
else
stop( paste(Unrecognised spline type, class(basis)[1]))
cutpoints - sort(unique(c(min(attr(basis,Boundary.knots)),
attr(basis,knots
npieces-length(cutpoints)
#evaluate derivatives and collect in array dd[order,piece]
dd-matrix(NA,degree+1,length(cutpoints))
if (class(basis)[1]==bs)
{
for (j in 0:degree)
dd[j+1,] - splineDesign(Aknots,cutpoints,ord=degree+1,outer=T,
deriv=rep(j,npieces))[,-1,drop=F] %*% coef
}
else if (class(basis)[1]==ns)
{
for (j in 0:degree)
{#Code from ns()
sdes - splineDesign(Aknots, cutpoints, 4,outer=T,
deriv=rep(j,npieces))[, -1, drop = FALSE]
const - splineDesign(Aknots, attr(basis,Boundary.knots), 4,
c(2, 2))[, -1, drop = FALSE]
qr.const - qr(t(const))
sdes - as.matrix((t(qr.qty(qr.const, t(sdes[, -(1:2)])
dd[j+1,] - sdes %*% coef
}
}
#add upper boundary knot to output list for convenience
cutpoints - c(cutpoints,max(attr(basis,Boundary.knots) ))
#create output coefficient lists
pars - list()
for ( i in 1:npieces)
{
pars[[i]] - dd[,i]/factorial(0:degree)
names(pars[[i]]) - c(1,paste(x^,1:degree,sep=))
names(pars)[i] - paste(cutpoints[i],to,cutpoints[i+1])
}
list(npieces=npieces,degree=degree,cutpoints=cutpoints,pars=pars)
}
# example
x - seq(10,30,len=10)
y - sin(x/2)+1
fit.bs - lm( y~bs(x,5) )
fit.ns - lm( y~ns(x,5) )
coef.bs - coef(fit.bs)[-1]
coef.ns - coef(fit.ns)[-1]
pieces.bs - pieces(attr(fit.bs$terms,predvars)[[3]],coef.bs)
pieces.ns - pieces(attr(fit.ns$terms,predvars)[[3]],coef.ns)
#verify in plot
plot(y~x)
xx - seq(min(x),max(x),len=3000)
lines(predict(fit.ns,new=data.frame(x=xx))~xx, lty=2)
lines(predict(fit.bs,new=data.frame(x=xx))~xx, lty=1)
abline(v=pieces.bs$cutpoints)
abline(v=pieces.ns$cutpoints,lty=2)
for (i in 1:pieces.bs$npieces)
{
xx - seq(pieces.bs$cutpoints[i],pieces.bs$cutpoints[i+1],len=1000)
mm - matrix(NA,length(xx),pieces.bs$degree+1)
for (j in 0:pieces.bs$degree) mm[,j+1] - (xx-pieces.bs$cutpoints[i])^j
yy - coef(fit.bs)[1] + mm %*% pieces.bs$pars[[i]]
lines(yy~xx, col=i+1)
}
for (i in 1:pieces.ns$npieces)
{
xx - seq(pieces.ns$cutpoints[i],pieces.ns$cutpoints[i+1],len=1000)
mm - matrix(NA,length(xx),pieces.ns$degree+1)
for (j in 0:pieces.ns$degree) mm[,j+1] -
[R] A question on colors for plotog groupedData
I have a groupedData object named data. When I use plot(data) I obtain a trellis plot for each group with a grey bakground. How can I change the background to white? I tried with par(bg=white) but got no change at all. I would appreciate any suggestion. -- Albert Sorribas Grup de Bioestadística i Biomatematica Departament de Ciències Mèdiques Bàsiques Universitat de Lledia tel: +34 973 702 406 FAX: +34 973 702 426 Home page: http://www.udl.es/Biomath/Group __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R is GNU S, not C.... [was how to get or store .....]
Well, this has been an interesting thread. I guess my own perspective
is warped, having never been a C programmer. My native languages are
FORTRAN, python, and R, all of which accept (or demand) a linefeed as a
terminator, rather than a semicolon, and two of which are very
particular about whitespace.
Accepting for a moment Ted's argument about wanting to compact his code,
the problem as I understand it is that ; is a statement separator, not
a statement terminator, so that Ted really does need all those ; in
between his statements on a single line, but that the last one implies a
NULL statement on every line. Given that the ; facilitates compacting
lines in vi (or vim), how about when the code is compacted you try
1,$s/;$//
which will remove all the final trailing ; leaving the other necessary
; separators.
Dave R.
(Ted Harding) wrote:
On 06-Dec-05 Martin Maechler wrote:
[But really, I'm more concerned and quite bit disappointed by
the diehard ; lovers]
Martin Maechler
Well, while not die-hard, I will put in my own little reason
for often using ; at the end of lines which don't need them.
Basically, this is done to protect me from myself (so in fact
is quite a strong reason).
I tend to develop extended R code in a side-window, using
a text editor (vim) in that window, and cutpasting the
chunks of R code from that window into the R window.
This usually means that I have a lot of short lines,
since it is easier when developing code to work with the
commands one per line, as they are easier to find and
less likely to be corrected erroneously.
Finally, when when I am content that the code does the job
I then put several short lines into one longer one.
For example (a function to do with sampling with probability
proportional to weights); first, as written line-by-line:
myfunction - function(X,n1,n2,n3,WTS){
N1-n1;
N2-n1+n2;
N3-n1+n2+n3;
# first selection
pii-WTS/sum(WTS);
alpha-N2;
Pi-alpha*pii;
r-runif(N3);
ix-sort(which(r=Pi));
# second selection
ix0-(1:N3);
ix3-ix0[-ix];
ix20-ix0[ix];
W-WTS[ix];
pii-W/sum(W);
Pi-N1*pii;
r-runif(length(Pi));
ix10-sort(which(r=Pi));
ix1-ix20[ix10];
ix2-ix20[-ix10];
# return the results
list(X1=X[ix1],X2=X[ix2],X3=X[ix3],ix1=ix1,ix2=ix2,ix3=ix3)
}
Having got that function right, with 'vim' in command mode
successive lines are readily brought up to the current line
by simply pressing J, which is very fast. This, in the
above case, then results in
MARselect-function(X,n1,n2,n3,WTS){
N1-n1; N2-n1+n2; N3-n1+n2+n3;
# first selection
pii-WTS/sum(WTS); alpha-N2; Pi-alpha*pii;
r-runif(N3); ix-sort(which(r=Pi));
# second selection
ix0-(1:N3); ix3-ix0[-ix]; ix20-ix0[ix];
W-WTS[ix]; pii-W/sum(W); Pi-N1*pii;
r-runif(length(Pi)); ix10-sort(which(r=Pi));
ix1-ix20[ix10]; ix2-ix20[-ix10];
# return the results
list(X1=X[ix1],X2=X[ix2],X3=X[ix3],ix1=ix1,ix2=ix2,ix3=ix3)
}
The greater readability of the first relative to the second is
obvious. The compactness of the second relative to the first
is evident. Obtaining the second from the first by repeated J
is very quick.
BUT -- if I had not put the ; at the ends of the lines in the
string-out version (which is easy to do as you type in the line
in the first place), then it would be much more trouble to get
the second version, and very easy to get it wrong!
Also, being long used to programming in C and octave/matlab,
putting ; at the end of a command is an easy reflex, and of
course does no harm at all to an R command.
Not that I'm trying to encourage others to do the same as I
do -- as I said, it's a self-protective habit -- but equally
if people (e.g. me) may find it useful I don't think it should
be discouraged either -- especially on aesthetic grounds!
Just my little bit ...
Best wishes,
Ted.
E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 094 0861
Date: 06-Dec-05 Time: 19:02:23
-- 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
--
David W. Roberts office 406-994-4548
Professor and Head FAX 406-994-3190
Department of Ecology email [EMAIL PROTECTED]
Montana State University
Bozeman, MT 59717-3460
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] A question on colors for plotog groupedData
On 12/7/05, Albert Sorribas [EMAIL PROTECTED] wrote: I have a groupedData object named data. When I use plot(data) I obtain a trellis plot for each group with a grey bakground. How can I change the background to white? I tried with par(bg=white) but got no change at all. I would appreciate any suggestion. Read ?trellis.device (after attaching the lattice package). Deepayan __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Time series influenced by half-time, intake and treatment...
Have you read Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer)? The latter part about nonlinear modeling with mixed effects sounds like it could help you a lot. 1. Consistent with that, I might start by averaging over all 15 people, then making plots and from the plots decide how to model everything else. 2. Then I'd try to fit that model to each person one at a time, as suggested by Pinheiro and Bates. 3. With the output from step 2, you can then use function nlme in package nlme. 4. I case this does not bring sufficient enlightenment and you would like more help from this listserve, I suggest you read the posting guide! www.R-project.org/posting-guide.html. I believe that questions more consistent with that posting guide tend to get more useful replies quicker. hope this helps. Kåre Edvardsen wrote: Hi! First of all: I'm a newbie to both statistics and R, so please be patient with me... I do however, like R because I've been programming (pascal, IDL, perl, C etc) and designing models since -92, but never related to statistics. Ok, here we go: I've got a set of 15 people, all of them observed over 10 weeks (10 analysed blood samples) with - let us kall it the A-value - influenced by: 1) a half time of the A-value of ~3 weeks 2) intake through diet (constant in time for each individ, but a big variation between the individs) 3) a treatment with a quick response, which may influence the A-value if sufficient dose is given. Problem: I do not know the limit for sufficient in 3). I also think there is a possibility the limit is individual, and thus may vary a lot between the individs. How do I approach such a problem? I know it's Friday, but I could not during the week figure out a model telling me if the treatment were significant or not... Anyone out there willing to guide me? Ked [[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 -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R is GNU S, not C.... [was how to get or store .....]
On 12/7/05, Dave Roberts [EMAIL PROTECTED] wrote:
Well, this has been an interesting thread. I guess my own perspective
is warped, having never been a C programmer. My native languages are
FORTRAN, python, and R, all of which accept (or demand) a linefeed as a
terminator, rather than a semicolon, and two of which are very
particular about whitespace.
Accepting for a moment Ted's argument about wanting to compact his code,
the problem as I understand it is that ; is a statement separator, not
a statement terminator, so that Ted really does need all those ; in
between his statements on a single line, but that the last one implies a
NULL statement on every line. Given that the ; facilitates compacting
lines in vi (or vim), how about when the code is compacted you try
1,$s/;$//
which will remove all the final trailing ; leaving the other necessary
; separators.
But watch out for this:
x - abc;
def
Dave R.
(Ted Harding) wrote:
On 06-Dec-05 Martin Maechler wrote:
[But really, I'm more concerned and quite bit disappointed by
the diehard ; lovers]
Martin Maechler
Well, while not die-hard, I will put in my own little reason
for often using ; at the end of lines which don't need them.
Basically, this is done to protect me from myself (so in fact
is quite a strong reason).
I tend to develop extended R code in a side-window, using
a text editor (vim) in that window, and cutpasting the
chunks of R code from that window into the R window.
This usually means that I have a lot of short lines,
since it is easier when developing code to work with the
commands one per line, as they are easier to find and
less likely to be corrected erroneously.
Finally, when when I am content that the code does the job
I then put several short lines into one longer one.
For example (a function to do with sampling with probability
proportional to weights); first, as written line-by-line:
myfunction - function(X,n1,n2,n3,WTS){
N1-n1;
N2-n1+n2;
N3-n1+n2+n3;
# first selection
pii-WTS/sum(WTS);
alpha-N2;
Pi-alpha*pii;
r-runif(N3);
ix-sort(which(r=Pi));
# second selection
ix0-(1:N3);
ix3-ix0[-ix];
ix20-ix0[ix];
W-WTS[ix];
pii-W/sum(W);
Pi-N1*pii;
r-runif(length(Pi));
ix10-sort(which(r=Pi));
ix1-ix20[ix10];
ix2-ix20[-ix10];
# return the results
list(X1=X[ix1],X2=X[ix2],X3=X[ix3],ix1=ix1,ix2=ix2,ix3=ix3)
}
Having got that function right, with 'vim' in command mode
successive lines are readily brought up to the current line
by simply pressing J, which is very fast. This, in the
above case, then results in
MARselect-function(X,n1,n2,n3,WTS){
N1-n1; N2-n1+n2; N3-n1+n2+n3;
# first selection
pii-WTS/sum(WTS); alpha-N2; Pi-alpha*pii;
r-runif(N3); ix-sort(which(r=Pi));
# second selection
ix0-(1:N3); ix3-ix0[-ix]; ix20-ix0[ix];
W-WTS[ix]; pii-W/sum(W); Pi-N1*pii;
r-runif(length(Pi)); ix10-sort(which(r=Pi));
ix1-ix20[ix10]; ix2-ix20[-ix10];
# return the results
list(X1=X[ix1],X2=X[ix2],X3=X[ix3],ix1=ix1,ix2=ix2,ix3=ix3)
}
The greater readability of the first relative to the second is
obvious. The compactness of the second relative to the first
is evident. Obtaining the second from the first by repeated J
is very quick.
BUT -- if I had not put the ; at the ends of the lines in the
string-out version (which is easy to do as you type in the line
in the first place), then it would be much more trouble to get
the second version, and very easy to get it wrong!
Also, being long used to programming in C and octave/matlab,
putting ; at the end of a command is an easy reflex, and of
course does no harm at all to an R command.
Not that I'm trying to encourage others to do the same as I
do -- as I said, it's a self-protective habit -- but equally
if people (e.g. me) may find it useful I don't think it should
be discouraged either -- especially on aesthetic grounds!
Just my little bit ...
Best wishes,
Ted.
E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 094 0861
Date: 06-Dec-05 Time: 19:02:23
-- 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
--
David W. Roberts office 406-994-4548
Professor and Head FAX 406-994-3190
Department of Ecology email [EMAIL PROTECTED]
Montana State University
Bozeman, MT 59717-3460
[R] organizing plot drawing routines; creating complex expressions
Hi,
My general goal is to find a coding strategy to efficiently store and
retrieve drawing routines for different plots.
This is my approach so far. In a single text file I store multiple drawing
routines where each routine draws a different plot.
A drawing routine may look like this:
panel.mean - function(means,...){
dots - list(...)
dots$x - means$data$trial.bin
dots$y - means$data$Freq
dots$group-
means$data$expected.hierarchyTo
dots$typ - l
dots$lty - 3
dots$lwd - 1
do.call(panel.superpose,dots)
}
arg - list()
arg$formula - formula(Freq~trial.bin|subj)
arg$group - quote(expected.hierarchyTo)
arg$data - atm.errors.intrusion.none[[frequency~hdt*prac*subj]]
arg$as.table - TRUE
arg$layout - c(3,4)
arg$col - c(green,blue,red)
arg$lwd - 2
arg$typ - b
arg$pch - 17
arg$lty - 1
arg$xlab - Trial Bin
arg$ylab - Frequency of Error Type
arg$main - ATM Experiment: Frequency of Intrusion-none Button
Presses During Error Trials\nAs A Function of Hierarchy Traversal Distance
To Current Expected Button Press
arg$scales - list(rot=c(90,90))
arg$sub - date()
arg$means -
list(data=gsummary(subset(arg$data,select=c(expected.hierarchyTo,trial.bi
n,Freq)),form=~expected.hierarchyTo/trial.bin,mean))
arg$par.settings -
list(background=list(col=white),strip.background=list(col=white))
arg$panel -
function(x,y,means,...){panel.mean(means,...);panel.superpose(x,y,...)}
print(do.call(xyplot,arg))
key -
draw.key(list(background=white,border=TRUE,text=list(lab=levels(arg$data$e
xpected.hierarchyTo)),lines=list(col=arg$col,pch=arg$pch,lty=arg$lty,lwd=arg
$lwd)))
vp-viewport(x=0.9,y=grobHeight(key)+unit(0.07,npc),just=c(right,top),
width=grobWidth(key),height=grobHeight(key))
pushViewport(vp)
grid.draw(key)
As you can see each routine is a sequence of commands required to draw a
single plot. Each routine is separated in the text file by comments. If I
wanted all these plots drawn I could source the file, but often I want only
one. Each time I want to draw one of these I 1) open the text file, 2)
scroll to the code for the desired plot, then copy and paste the code into
the R console. Is there a more efficient way to store and retrieve these
drawing routines?
For example, I am exploring the following idea. Save the plot routines as
expressions in a list object. The result is a list, l, with x components,
where each component contains the routine to draw a plot as an expression.
To draw a single plot all I need to do is type eval(l[[x]]).
While this approach has succeeded, it has some drawbacks. One drawback comes
about because I also want to be able to copy and paste directly from the
routine as it is written into the R console. If I use the expression command
to convert the routine to an expression, then each statement must be
separated by commas. Using that syntax I cannot copy and paste command1
through command3 (see below) to the console because the commas will be
included and cause a syntax error.
expression(
command1,
command2,
command3,
ect.
)
Another option is to call parse(file=) and then run each command in
sequence without commas. The problem here is that parse does not work on
multi-line functions at the command line (see below). Why doesn't this work?
parse(file=)
?function(x){x-9}
expression(function(x) {
x - 9
})
parse(file=)
?function(x){x-9;x-8}
Error in parse(file, n, text, prompt) : parse error
parse(file=)
?function(x){x-9
Error in parse(file, n, text, prompt) : parse error
Using parse also requires knowing how many statements will be processed
ahead of time. If I go this route I would like some way to automate that as
well.
Assuming the parse problem described above cannot be fixed, there are two
other possible options.
1) Rather than have multiple drawing routines in a single file, have each
drawing routine in it's own file. In that case the syntax in the file is
just as it would be at the command line, which is the way I want it.
Furthermore, calling prase with a file now works--that is, it can handle
multiple line functions with ease. The only real draw back to this is that
having separate files for each and every plot may be overwhelming. I'd
rather have a series of related plots in the same text file than the same
folder. But, I could get over that.
2) Convert the drawing routines to functions without arguments--not
expressions. In this scheme each component of the list would contain a
function with the drawing
Re: [R] Change labels of x-axes in Plot of stl() function?
Thanks a lot!
However, its not working perfectly yet. The function plot.stl now also
changes the labels of Data, Seasonal, Trend, and Remainder to the text
defined for xlab. How could this be fine-tuned?
Regards,
Jan
-Original Message-
From: Gabor Grothendieck [mailto:[EMAIL PROTECTED]
Sent: Wednesday, December 07, 2005 5:07 PM
To: Jan Verbesselt
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] Change labels of x-axes in Plot of stl() function?
If you look through the output of:
stats:::plot.stl
you see right near the end that time is hard coded in the call to mtext.
However, we could temporarily redefine mtext so that it works as you
wish and then redefine plot.stl so that it looks within the environment
of our function to find mtext (rather than looking in the stats package):
plot.stl - function(..., xlab = time) {
mtext - function(text, ...) graphics::mtext(xlab, ...)
plot.stl - stats:::plot.stl
environment(plot.stl) - environment()
plot.stl(...)
}
# test it
example(stl)
plot.stl(stmd, xlab = X)
On 12/7/05, Jan Verbesselt [EMAIL PROTECTED] wrote:
Hi all,
How can the label of the x-axes in the plot() of a stl.object be
adapted?
e.g.,
When plotting: plot(stl(nottem, per))
In the labels of the x-axes is time. How can this be changed to e.g.,
Time (dekade) ?
It does not work with xlab or others anymore
Thanks,
Jan
___
Ir. Jan Verbesselt
Research Associate
Biosystems Department ~ M³-BIORES
Vital Decosterstraat 102, 3000 Leuven, Belgium
Tel: +32-16-329750 Fax: +32-16-329760
http://gloveg.kuleuven.ac.be/
___
Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
[[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
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
Re: [R] Change labels of x-axes in Plot of stl() function?
One other thought. This can be arguably done more compactly
using the proto package. We define a proto object consisting
of three components:
- plot.stl which is just a copy of the corresponding routine
in the stats package,
- our redefined mtext and
- our desired x label.
and then run plot.stl in the scope of the proto object. There
are two simplifications here:
1. we don't need to do explicit manipulations of environments
since proto automatically resets the environments of component
functions. In particular, the environment of plot.stl is automatically
reset to point to the proto object so that its scope is no longer
the stats package.
2. we don't need to redefine the argument list of plot.stl since
with its environment automatically redefined as just explained,
plot.stl will look for xlab in the proto object so we can just put it
there rather than pass it in a new, different, argument list
library(proto)
STL - proto(plot.stl = stats:::plot.stl,
mtext = function(text, ...) graphics::mtext(xlab, ...),
xlab = time (decade))
# test. First line computes stmd to be used as test data.
example(stl)
with(STL, plot.stl(stmd))
# note that to change the xlab we can do this
STL$xlab - New xlab
with(STL, plot.stl(stmd))
On 12/7/05, Gabor Grothendieck [EMAIL PROTECTED] wrote:
If you look through the output of:
stats:::plot.stl
you see right near the end that time is hard coded in the call to mtext.
However, we could temporarily redefine mtext so that it works as you
wish and then redefine plot.stl so that it looks within the environment
of our function to find mtext (rather than looking in the stats package):
plot.stl - function(..., xlab = time) {
mtext - function(text, ...) graphics::mtext(xlab, ...)
plot.stl - stats:::plot.stl
environment(plot.stl) - environment()
plot.stl(...)
}
# test it
example(stl)
plot.stl(stmd, xlab = X)
On 12/7/05, Jan Verbesselt [EMAIL PROTECTED] wrote:
Hi all,
How can the label of the x-axes in the plot() of a stl.object be adapted?
e.g.,
When plotting: plot(stl(nottem, per))
In the labels of the x-axes is time. How can this be changed to e.g.,
Time (dekade) ?
It does not work with xlab or others anymore…
Thanks,
Jan
___
Ir. Jan Verbesselt
Research Associate
Biosystems Department ~ M³-BIORES
Vital Decosterstraat 102, 3000 Leuven, Belgium
Tel: +32-16-329750 Fax: +32-16-329760
http://gloveg.kuleuven.ac.be/
___
Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
[[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
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] R programming job in Boston
Hello, I'm a senior statistician at Affinnova, a market research firm outside of Boston, MA. We are looking to hire a strong R programmer. Here is the description of the job: Come join a growing company in the market innovation field. Affinnova is changing innovation and product development by applying Interactive Genetic Algorithms to bridge the chasm between the worlds of the designer, the marketer, and the consumer. We work with large Consumer Products Goods companies and retail clients. See our website www.affinnova.com for more information on this job under the Company - Careers section. We are seeking statisticians with R or S-Plus experience. The Statistician is a vital member of a highly talented cross-functional team that delivers market research projects to our clients. Responsibilities * Work with state-of-the-art data collection, analysis, and visualization strategies for product and service optimization. * Collaborate company-wide with project managers, sales, product managers, etc. and the client to identify business initiatives, and solve business problems * Help develop and design new statistical tools Skills and Experience * MS or PhD in statistics, biostatistics, statistical genetics, econometrics, or other quantitative fields * 3+ years experience in research design and statistical modeling in a business environment (marketing, market research, etc.) * Modeling experience in linear and nonlinear regression, logistic regression, cluster analysis or other multivariate methods required * Excellent written, oral, and presentation skills * Experience with latent-class models * Experience with conjoint and discrete choice models * Experience with Bayesian models and simulation models * Proficiency with R, S-Plus a must, further experience with SAS, Stata, Matlab, and SPSS a plus. * Knowledge of multi-dimensional visualization techniques. * Ability to translate business problems into technical solutions, and translate results into visuals the client can understand * Success at working on cross-functional teams to meet a common goal * Experience communicating with internal and/or external customers * Creative and independent thinker. Strong initiative * Experience working in Consumer Packaged Goods and/or other consumer-oriented industries * Experience with survey design and delivery, Exposure to genetic and evolutionary computation, Proficiency with Excel (graphical functions and VBA), Experience with Java, SQL, VBA all a plus. For information on how to apply, please see our website. Cheers, Oliver Will, PhD Senior Statistician Affinnova, Inc 52 Second Ave 2nd Floor South Waltham, MA 02451 USA Phone: 781-464-4751 Fax: 781-464-4701 Cell: 978-844-1176 www.affinnova.com [[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
Re: [R] Change labels of x-axes in Plot of stl() function?
This should fix that problem:
plot.stl - function(..., xlab = time) {
mtext - function(text, ...)
graphics::mtext(if (text == time) xlab else text, ...)
plot.stl - stats:::plot.stl
environment(plot.stl) - environment()
plot.stl(...)
}
Also for the proto solution try this:
library(proto)
STL - proto(plot.stl = stats:::plot.stl,
mtext = function(text, ...)
graphics::mtext(if (text == time) xlab else text, ...),
xlab = time (decade))
# test. First line computers stmd to be used as test data.
example(stl)
with(STL, plot.stl(stmd))
On 12/7/05, Jan Verbesselt [EMAIL PROTECTED] wrote:
Thanks a lot!
However, it's not working perfectly yet. The function plot.stl now also
changes the labels of Data, Seasonal, Trend, and Remainder to the text
defined for xlab. How could this be fine-tuned?
Regards,
Jan
-Original Message-
From: Gabor Grothendieck [mailto:[EMAIL PROTECTED]
Sent: Wednesday, December 07, 2005 5:07 PM
To: Jan Verbesselt
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] Change labels of x-axes in Plot of stl() function?
If you look through the output of:
stats:::plot.stl
you see right near the end that time is hard coded in the call to mtext.
However, we could temporarily redefine mtext so that it works as you
wish and then redefine plot.stl so that it looks within the environment
of our function to find mtext (rather than looking in the stats package):
plot.stl - function(..., xlab = time) {
mtext - function(text, ...) graphics::mtext(xlab, ...)
plot.stl - stats:::plot.stl
environment(plot.stl) - environment()
plot.stl(...)
}
# test it
example(stl)
plot.stl(stmd, xlab = X)
On 12/7/05, Jan Verbesselt [EMAIL PROTECTED] wrote:
Hi all,
How can the label of the x-axes in the plot() of a stl.object be
adapted?
e.g.,
When plotting: plot(stl(nottem, per))
In the labels of the x-axes is time. How can this be changed to e.g.,
Time (dekade) ?
It does not work with xlab or others anymore…
Thanks,
Jan
___
Ir. Jan Verbesselt
Research Associate
Biosystems Department ~ M³-BIORES
Vital Decosterstraat 102, 3000 Leuven, Belgium
Tel: +32-16-329750 Fax: +32-16-329760
http://gloveg.kuleuven.ac.be/
___
Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
[[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
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
[R] concatenate data frame
hi all Here is a small part of my code: tab_tmp-tab[1:(no[off_set[i-1]+1]+(no[off_set[i]+1]-no[off_set[i-1]+1])),length(tab)]; tab_tmp1-tab[(no[off_set[i-1]+1]+(no[off_set[i]+1]-no[off_set[i-1]+1])):length(TotalFillTimeHours),length(tab)]; tab-c(tab_tmp,tab_tmp1); attach(tab); Here is the output: Error in attach(tab) : attach only works for lists and data frames Execution halted How do i concatenate them in order to keep the data frame structure? thks for your help guillaume __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Maintaining factors when copying from one data frame to another
Greetings all: OK, this is bugging the @[EMAIL PROTECTED] out of me. I know the answer is simple and straightforward but for the life of me I cannot find it in the documentation, in the archives, or in my notes (because I know I've encountered this in the past). My problem is: I have a data frame with columns A, B, C, D, and E. A, B, and E are factors and C and D are numeric. I need a new data frame with just A, C, and D. When I copy these columns to a new data frame newDF - data.frame(cbind(oldDF$A, oldDF$C, oldDF$D)) all the factor data comes out as levels rather than the original factors. How do I preserve the factors when I copy from one data frame to another? Thanks vary much, Kurt Wollenberg, Ph.D. Tufts Center for Vision Research Tufts-New England Medical Center 750 Washington St #450 Boston, MA 02111 Office: 617-636-9028 Fax: 617-636-8945 email: kurt.wollenberg at gmail dot com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] ttda on R 2.1.1: error
Do you know where I can download R 2.1.0 ?? I would like to use ttda. Thanks Lucie [[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
Re: [R] Constructing a transition matrix
Hi Peter, Thanks for pointing out the set functions. I can use setdiff to find missing rows setdiff(dev, rownames(A)) [1] seed and intersect to find common rows d1- intersect(dev, rownames(A) ) [1] veg rep I was trying to use a negative index like A[-1,] to remove the dead row, but d1 is a better solution. Now I can add the missing seed row and get a square matrix. rbind( seed=numeric(3), A[d1,] )[dev,dev] Another post by Hans Gardfjell suggested reordering factor levels before using prop.table(table()) and this solution works great! trans$class - ordered(trans$class, levels=dev) trans$fate - ordered(trans$fate, levels=c(dev,dead) ) A - t(prop.table(table(trans$class, trans$fate),1))[-4,] seed veg rep seed 0.000 0 0.0 veg 0.667 0 0.5 rep 0.000 1 0.5 Thanks for the help, Chris Peter Dalgaard wrote: Are you looking for something like d1 - setdiff(dev,seed) A0[d1,dev] - A[d1,dev] ? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Maintaining factors when copying from one data frame to another
Does newDF - oldDF[,c(A,C,D)] work? Kurt Wollenberg wrote: Greetings all: OK, this is bugging the @[EMAIL PROTECTED] out of me. I know the answer is simple and straightforward but for the life of me I cannot find it in the documentation, in the archives, or in my notes (because I know I've encountered this in the past). My problem is: I have a data frame with columns A, B, C, D, and E. A, B, and E are factors and C and D are numeric. I need a new data frame with just A, C, and D. When I copy these columns to a new data frame newDF - data.frame(cbind(oldDF$A, oldDF$C, oldDF$D)) all the factor data comes out as levels rather than the original factors. How do I preserve the factors when I copy from one data frame to another? Thanks vary much, Kurt Wollenberg, Ph.D. Tufts Center for Vision Research Tufts-New England Medical Center 750 Washington St #450 Boston, MA 02111 Office: 617-636-9028 Fax: 617-636-8945 email: kurt.wollenberg at gmail dot com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Kevin E. Thorpe Biostatistician/Trialist, Knowledge Translation Program Assistant Professor, Department of Public Health Sciences Faculty of Medicine, University of Toronto email: [EMAIL PROTECTED] Tel: 416.946.8081 Fax: 416.946.3297 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Maintaining factors when copying from one data frame to another
newDF - oldDF[,c(A,C,D)] HTH Kurt Wollenberg wrote: Greetings all: OK, this is bugging the @[EMAIL PROTECTED] out of me. I know the answer is simple and straightforward but for the life of me I cannot find it in the documentation, in the archives, or in my notes (because I know I've encountered this in the past). My problem is: I have a data frame with columns A, B, C, D, and E. A, B, and E are factors and C and D are numeric. I need a new data frame with just A, C, and D. When I copy these columns to a new data frame newDF - data.frame(cbind(oldDF$A, oldDF$C, oldDF$D)) all the factor data comes out as levels rather than the original factors. How do I preserve the factors when I copy from one data frame to another? Thanks vary much, Kurt Wollenberg, Ph.D. Tufts Center for Vision Research Tufts-New England Medical Center 750 Washington St #450 Boston, MA 02111 Office: 617-636-9028 Fax: 617-636-8945 email: kurt.wollenberg at gmail dot com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Ferdinand Alimadhi Programmer / Analyst Harvard University The Institute for Quantitative Social Science (617) 496-0187 [EMAIL PROTECTED] www.iq.harvard.edu [[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
Re: [R] Maintaining factors when copying from one data frame to another
Kurt Wollenberg wrote: Greetings all: OK, this is bugging the @[EMAIL PROTECTED] out of me. I know the answer is simple and straightforward but for the life of me I cannot find it in the documentation, in the archives, or in my notes (because I know I've encountered this in the past). My problem is: I have a data frame with columns A, B, C, D, and E. A, B, and E are factors and C and D are numeric. I need a new data frame with just A, C, and D. When I copy these columns to a new data frame newDF - data.frame(cbind(oldDF$A, oldDF$C, oldDF$D)) all the factor data comes out as levels rather than the original factors. How do I preserve the factors when I copy from one data frame to another? The following may be better: newDF - subset(oldDF, select = c(A, C, D)) HTH, Tobias __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Maintaining factors when copying from one data frame to another
On Wed, 7 Dec 2005, Kurt Wollenberg wrote: Greetings all: OK, this is bugging the @[EMAIL PROTECTED] out of me. I know the answer is simple and straightforward but for the life of me I cannot find it in the documentation, in the archives, or in my notes (because I know I've encountered this in the past). My problem is: I have a data frame with columns A, B, C, D, and E. A, B, and E are factors and C and D are numeric. I need a new data frame with just A, C, and D. When I copy these columns to a new data frame newDF - data.frame(cbind(oldDF$A, oldDF$C, oldDF$D)) set.seed(1) oldDF - data.frame(A=rnorm(5), B=rnorm(5), C=factor(letters[1:5]), D=rnorm(5)) str(oldDF) newDF - data.frame(cbind(oldDF$A, oldDF$C, oldDF$D)) str(newDF) # cbind forces to numeric as a matrix, so avoid it here newDFa - data.frame(Aa=oldDF$A, Ca=oldDF$C, Da=oldDF$D) str(newDFa) all the factor data comes out as levels rather than the original factors. How do I preserve the factors when I copy from one data frame to another? Thanks vary much, Kurt Wollenberg, Ph.D. Tufts Center for Vision Research Tufts-New England Medical Center 750 Washington St #450 Boston, MA 02111 Office: 617-636-9028 Fax: 617-636-8945 email: kurt.wollenberg at gmail dot com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R is GNU S, not C.... [was how to get or store .....]
On Tue, Dec 06, 2005 at 04:21:01PM +, Patrick Burns wrote: I don't put in extraneous ';' because I maybe get a blister on my little finger. I suspect that those who find the semi-colons ugly in R do not find them ugly in C. I think the reason there would be a visceral reaction in R but not in C is that there is a danger when using them in R that they really mean something. We get questions on R-help often enough about why code like: if(x 0) y - 4 else y - 4.5e23 doesn't work. If people habitually used semi-colons, those sorts of questions would probably multiply. I don't understand. It would seem to me that in if (x 0) y - 4; else y - 4.5e23; it's pretty obvious that the if statement is terminated by the semicolon at the end of the first line and that therefore, the else on the next line is erroneous because it is not associated with any if. At least, the version above fails consistently, i.e. regardless of context. On the other hand, I've studied the R Language Definition for quite some time before fully understanding why if (x 0) y - 4 else y - 4.5e23 works inside of a function (or other enclosing block) while it does not work interactively (or at the top level of a script). Best regards, Jan -- +- Jan T. Kim ---+ | email: [EMAIL PROTECTED] | | WWW: http://www.cmp.uea.ac.uk/people/jtk | *-= hierarchical systems are for files, not for humans =-* __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Are minbucket and minsplit rpart options working as expected?
Dear r-list: I am using rpart to build a tree on a dataset. First I obtain a perhaps too large tree: arbol.bsvg.02 - rpart(formula, data = bsvg, subset=grp.entr, control=rpart.control(cp=0.001)) arbol.bsvg.02 n= 10 node), split, n, loss, yval, (yprob) * denotes terminal node 1) root 10 6657 0 (0.93343000 0.06657000) 2) meses_antiguedad_svg=10.5 73899 3658 0 (0.95050001 0.0494) 4) eor_n1_gns 1.5 63968 2807 0 (0.95611868 0.04388132) 8) tarifa_gas=31,32,33,34 63842 2771 0 (0.95659597 0.04340403) * 9) tarifa_gas=NO 126 36 0 (0.71428571 0.28571429) 18) tipo_mercado=ESP,N/A 90 10 0 (0.8889 0.) * 19) tipo_mercado=NE ,SAH,SAV 36 10 1 (0.2778 0.7222) * 5) eor_n1_gns=1.5 9931 851 0 (0.91430873 0.08569127) 10) sn_calef=0.5 8390 546 0 (0.93492253 0.06507747) * 11) sn_calef 0.5 1541 305 0 (0.80207657 0.19792343) 22) tarifa_gas=31,NO 1134 141 0 (0.87566138 0.12433862) * 23) tarifa_gas=32 407 164 0 (0.59705160 0.40294840) 46) cons_gas_delta_1 6997 196 51 0 (0.73979592 0.26020408) * 47) cons_gas_delta_1=6997 211 98 1 (0.46445498 0.53554502) 94) meses_antiguedad_svg=23.5 134 54 0 (0.59701493 0.40298507) 188) altitud 312 61 16 0 (0.73770492 0.26229508) * 189) altitud=312 73 35 1 (0.47945205 0.52054795) 378) back_office=1.5 39 12 0 (0.69230769 0.30769231) * 379) back_office 1.5 348 1 (0.23529412 0.76470588) * 95) meses_antiguedad_svg 23.5 77 18 1 (0.23376623 0.76623377) * 3) meses_antiguedad_svg 10.5 26101 2999 0 (0.88510019 0.11489981) 6) sn_calef=0.5 20129 1853 0 (0.90794376 0.09205624) * 7) sn_calef 0.5 5972 1146 0 (0.80810449 0.19189551) 14) tarifa_gas=31 4406 664 0 (0.84929641 0.15070359) * 15) tarifa_gas=32,NO 1566 482 0 (0.69220945 0.30779055) 30) eor_n1_gns 0.5 1168 306 0 (0.73801370 0.26198630) * 31) eor_n1_gns=0.5 398 176 0 (0.55778894 0.44221106) 62) back_office=0.5 148 35 0 (0.76351351 0.23648649) * 63) back_office 0.5 250 109 1 (0.4360 0.5640) * So I decide not to consider branches with less than 1000 observations, a 1% of the original number of observations. Therefore, according to the rpart.control help pages, I set minbucket=1000. However, arbol.bsvg.02 n= 10 node), split, n, loss, yval, (yprob) * denotes terminal node 1) root 10 6657 0 (0.9334300 0.0665700) * And I get an empty tree. But there were branches in the original tree with more than 1000 observations. Something similar happens if I set minsplit (or both minbucket and minsplit) to a similar value: I end up with the same root, branch-less tree. Am I misreading something? Can anybody cast a light on the correct usage of the minbucket (and/or minsplit) for me? Sincerely, Carlos J. Gil Bellosta http://www.datanalytics.com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] concatenate data frame
Guillaume, I assume that 'tab' is a data frame and that, for some unspecified reason, you want to get two subsets of the last column of tab, overlapping one case, and coercing the final result to a data frame. If that is correct, then as.data.frame(c(tab_tmp, tab_tmp1)) will give you a data frame. Alternatively, check out the 'drop =' argument to '[.data.frame' and then rbind your pieces. Peter Ehlers [EMAIL PROTECTED] wrote: hi all Here is a small part of my code: tab_tmp-tab[1:(no[off_set[i-1]+1]+(no[off_set[i]+1]-no[off_set[i-1]+1])),length(tab)]; tab_tmp1-tab[(no[off_set[i-1]+1]+(no[off_set[i]+1]-no[off_set[i-1]+1])):length(TotalFillTimeHours),length(tab)]; tab-c(tab_tmp,tab_tmp1); attach(tab); Here is the output: Error in attach(tab) : attach only works for lists and data frames Execution halted How do i concatenate them in order to keep the data frame structure? thks for your help guillaume __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] KMO sampling adequacy and SPSS -- partial solution
Dear colleagues,
I've been searching for information on the Kaiser-Meyer-Olkin (KMO)
Measure of Sampling Adequacy (MSA). This statistic is generated in
SPSS and is often used to determine if a dataset is appropriate for
factor analysis -- it's true utility seems quite low, but it seems to
come up in stats classes a lot. It did in mine, and a glance through
the R-help archives suggests I'm not alone.
I finally found a reference describing the calculation, and wrote the
following R function to perform it. Note that the function depends on
a partial correlation function found in library(corpcor).
kmo.test - function(df){
###
## Calculate the Kaiser-Meyer-Olkin Measure of Sampling Adequacy.
## Input should be a data frame or matrix, output is the KMO statistic.
## Formula derived from Hutcheson et al, 1999,
## The multivariate social scientist, page 224, ISBN 0761952012
## see http://www2.chass.ncsu.edu/garson/pa765/hutcheson.htm
###
cor.sq = cor(df)^2
cor.sumsq = (sum(cor.sq)-dim(cor.sq)[1])/2
library(corpcor)
pcor.sq = cor2pcor(cor(df))^2
pcor.sumsq = (sum(pcor.sq)-dim(pcor.sq)[1])/2
kmo = sus.cor.ss/(sus.cor.ss+sus.pcor.ss)
return(kmo)
}
Also, for those trying to reproduce the SPSS factor analysis output,
(-1 * cor2pcor(cor(yourDataFrame))) will produce the anti-image
correlation matrix. Unfortunately, the most useful property of that
matrix in SPSS is that the diagonals represent the individual MSA
values -- I haven't found a way to derive those yet. Still working on
that, any suggestions appreciated.
--Ash.
-
Ashish Ranpura
Institute of Cognitive Neuroscience
University College London
17 Queen Square
London WC1N 3AR
tel: +44 (20) 7679 1126
web: http://www.icn.ucl.ac.uk
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] KMO sampling adequacy and SPSS -- partial solution
Sorry, there was an error in that function, a hangover from a
previous session. The corrected function is:
kmo.test - function(df){
###
## Calculate the Kaiser-Meyer-Olkin Measure of Sampling Adequacy.
## Input should be a data frame or matrix, output is the KMO statistic.
## Formula derived from Hutcheson et al, 1999,
## The multivariate social scientist, page 224, ISBN 0761952012
## see http://www2.chass.ncsu.edu/garson/pa765/hutcheson.htm
###
cor.sq = cor(df)^2
cor.sumsq = (sum(cor.sq)-dim(cor.sq)[1])/2
library(corpcor)
pcor.sq = cor2pcor(cor(df))^2
pcor.sumsq = (sum(pcor.sq)-dim(pcor.sq)[1])/2
kmo = cor.sumsq/(cor.sumsq+pcor.sumsq)
return(kmo)
}
--Ashish.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] ttda on R 2.1.1: error
Lucie RYBARCZYK [EMAIL PROTECTED] writes: Do you know where I can download R 2.1.0 ?? I would like to use ttda. Which platform? Both source and Windows binaries of historic releases are on CRAN and are not particularly hard to find. Any particular reason why this cannot work with R-2.2.0 (or 2.2.1 beta for that matter)? Thanks Lucie -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] contrasts for lm
I would like estimate a number of contrasts from a one-way ANOVA model. I see that the lm command has a contrasts option, but I can't figure out how to use it! Any help that can be offered would be greatly apreciated. Here is my model statement: Model-lm(log2PM~P+T+P*T) where P has 16 levels, T(treatment) has 12 levels and I am interested in looking at different treatment comparisons. Thanks! Ann __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] error message from building R in Mac
Following the instruction in R website and downloading gcc and g77, I am trying to configure and build R in my Mac laptop, but got some error message that I do not know how to resolve. Do any of you know how to solve this problem? After type ./configure, I got the following message. R is now configured for powerpc-apple-darwin8.3.0 Source directory: . Installation directory:/Library/Frameworks C compiler:gcc -g -O2 C++ compiler: g++ -g -O2 Fortran compiler: gfortran -g -O2 Interfaces supported: X11, aqua, tcltk External libraries:readline, BLAS(generic), LAPACK(in blas) Additional capabilities: iconv, MBCS, NLS Options enabled: framework, R profiling Recommended packages: yes But after typing make, I got the following error message. make[1]: Nothing to be done for `R'. make[1]: Nothing to be done for `R'. make[2]: Nothing to be done for `R'. creating src/scripts/R.fe config.status: creating src/include/config.h config.status: src/include/config.h is unchanged Rmath.h is unchanged make[3]: Nothing to be done for `R'. make[4]: `libbz2.a' is up to date. make[4]: `libpcre.a' is up to date. make[4]: `libz.a' is up to date. make[3]: Nothing to be done for `R'. make[3]: `stamp-lo' is up to date. make[3]: `stamp-lo' is up to date. make[3]: `stamp-lo' is up to date. /Users/fang-yichiou/tmp/R-2.2.0/lib/libR.dylib is unchanged /Users/fang-yichiou/tmp/R-2.2.0/bin/exec/R is unchanged make[4]: `R_X11.so' is up to date. make[4]: `internet.so' is up to date. make[4]: `lapack.so' is up to date. make[4]: `vfonts.so' is up to date. building system startup profile building package 'base' all.R is unchanged building package 'tools' all.R is unchanged make[5]: `Makedeps' is up to date. ../../../../library/tools/libs/tools.so is unchanged building package 'utils' all.R is unchanged building package 'grDevices' all.R is unchanged ../../../library/grDevices/R/grDevices is unchanged make[5]: `Makedeps' is up to date. ../../../../library/grDevices/libs/grDevices.so is unchanged Error in solve.default(rgb) : lapack routines cannot be loaded In addition: Warning message: unable to load shared library '/Users/fang-yichiou/tmp/R-2.2.0/ modules/lapack.so': dlopen(/Users/fang-yichiou/tmp/R-2.2.0/modules/lapack.so, 6): Symbol not found: _rcblas_zdotc_sub__ Referenced from: /Users/fang-yichiou/tmp/R-2.2.0/modules/lapack.so Expected in: flat namespace Error: unable to load R code in package 'grDevices' Execution halted make[3]: *** [all] Error 1 make[2]: *** [R] Error 1 make[1]: *** [R] Error 1 make: *** [R] Error 1 Thanks for your advice in advance. Fang-Yi __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Change labels of x-axes in Plot of stl() function?
I noticed that we can combine the function and proto
approaches by placing the proto in the function with
these advantages:
1. the function body can be reduced to just two statements
2. no explicit manipulation of environments via
environment(...) is required (as proto does that itself
automatically)
In our new solution, the function signature and the
redefinition of mtext are the same as in our prior function
solution but the remaining lines in the function solution
(that manipulate environments) are replaced with the single
'with' command as shown.
As before, placing stats:::plot.stl in the proto results in
(1) a copy of that function being placed in the proto object and
(2) the environment of that copy being set to the proto object itself
The parent of that proto object defaults to its lexical environment
so all we need to do in order to ensure that xlab and the new mtext
are accessible from the copy of plot.stl are to ensure that they
are in that parent environment -- they need not be in the proto
object itself. They will be accessed via inheritance anyways. Thus
the solution reduces to just:
library(proto)
plot.stl - function(..., xlab = time) {
mtext - function(text, ...)
graphics::mtext(if (text == time) xlab else text, ...)
with( proto(plot.stl = stats:::plot.stl), plot.stl(...) )
}
# test
example(stl) # defines stdm for use in the test
plot.stl(stdm, xlab = My xlab)
On 12/7/05, Gabor Grothendieck [EMAIL PROTECTED] wrote:
This should fix that problem:
plot.stl - function(..., xlab = time) {
mtext - function(text, ...)
graphics::mtext(if (text == time) xlab else text, ...)
plot.stl - stats:::plot.stl
environment(plot.stl) - environment()
plot.stl(...)
}
Also for the proto solution try this:
library(proto)
STL - proto(plot.stl = stats:::plot.stl,
mtext = function(text, ...)
graphics::mtext(if (text == time) xlab else text, ...),
xlab = time (decade))
# test. First line computers stmd to be used as test data.
example(stl)
with(STL, plot.stl(stmd))
On 12/7/05, Jan Verbesselt [EMAIL PROTECTED] wrote:
Thanks a lot!
However, it's not working perfectly yet. The function plot.stl now also
changes the labels of Data, Seasonal, Trend, and Remainder to the text
defined for xlab. How could this be fine-tuned?
Regards,
Jan
-Original Message-
From: Gabor Grothendieck [mailto:[EMAIL PROTECTED]
Sent: Wednesday, December 07, 2005 5:07 PM
To: Jan Verbesselt
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] Change labels of x-axes in Plot of stl() function?
If you look through the output of:
stats:::plot.stl
you see right near the end that time is hard coded in the call to mtext.
However, we could temporarily redefine mtext so that it works as you
wish and then redefine plot.stl so that it looks within the environment
of our function to find mtext (rather than looking in the stats package):
plot.stl - function(..., xlab = time) {
mtext - function(text, ...) graphics::mtext(xlab, ...)
plot.stl - stats:::plot.stl
environment(plot.stl) - environment()
plot.stl(...)
}
# test it
example(stl)
plot.stl(stmd, xlab = X)
On 12/7/05, Jan Verbesselt [EMAIL PROTECTED] wrote:
Hi all,
How can the label of the x-axes in the plot() of a stl.object be
adapted?
e.g.,
When plotting: plot(stl(nottem, per))
In the labels of the x-axes is time. How can this be changed to e.g.,
Time (dekade) ?
It does not work with xlab or others anymore…
Thanks,
Jan
___
Ir. Jan Verbesselt
Research Associate
Biosystems Department ~ M³-BIORES
Vital Decosterstraat 102, 3000 Leuven, Belgium
Tel: +32-16-329750 Fax: +32-16-329760
http://gloveg.kuleuven.ac.be/
___
Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
[[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
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
Re: [R] R is GNU S, not C.... [was how to get or store .....]
On 12/7/05, Jan T. Kim [EMAIL PROTECTED] wrote:
On Tue, Dec 06, 2005 at 04:21:01PM +, Patrick Burns wrote:
I don't put in extraneous ';' because I maybe get a
blister on my little finger.
I suspect that those who find the semi-colons ugly in
R do not find them ugly in C. I think the reason there
would be a visceral reaction in R but not in C is that
there is a danger when using them in R that they really
mean something.
We get questions on R-help often enough about why
code like:
if(x 0) y - 4
else y - 4.5e23
doesn't work.
If people habitually used semi-colons, those sorts of
questions would probably multiply.
I don't understand. It would seem to me that in
if (x 0) y - 4;
else y - 4.5e23;
it's pretty obvious that the if statement is terminated by the
semicolon at the end of the first line and that therefore, the else
on the next line is erroneous because it is not associated with any
if.
Why is it obvious?
if (x 0) y = 4;
else y = 4.5e23;
is perfectly legal C.
At least, the version above fails consistently, i.e. regardless of
context.
It fails for unrelated reasons.
if (x 0) { y - 4; } else { y - 4.5e23; }
doesn't fail.
-Deepayan
On the other hand, I've studied the R Language Definition for
quite some time before fully understanding why
if (x 0) y - 4
else y - 4.5e23
works inside of a function (or other enclosing block) while it does not
work interactively (or at the top level of a script).
Best regards, Jan
--
+- Jan T. Kim ---+
| email: [EMAIL PROTECTED] |
| WWW: http://www.cmp.uea.ac.uk/people/jtk |
*-= hierarchical systems are for files, not for humans =-*
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] mle.stepwise versus step/stepAIC
Hello, I have a question pertaining to the stepwise regression which I am trying to perform. I have a data set in which I have 14 predictor variables accompanying my response variable. I am not sure what the difference is between the function mle.stepwise found in the wle package and the functions step or stepAIC? When would one use mle.stepwise versus step/stepAIC? Other than the references in the R-software, is there anything which discusses the use of mle.stepwise? Thank you kindly in advance for any insight offered. Jim Brindle __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Ancova and lme use
Mon cher M. MENICACCI: It looks to me like you ultimately want to use lmer in library(lme4) [which also requires library(Matrix)]. For documentation, I suggest you start with Doug Bates (2005) Fitting Linear Mixed Models in R, R News, vol. 5/1: 27-30 (available from www.r-project.org - Newsletter). After install.packages(lme4), I suggest you read Implementation.pdf in ~R\library\lme4\doc. I also suggest you get Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer). For me, this book was essential documentation for lme, the previous implementation of lmer. Studying that book might help you understand lmer. Also, I encourage you to use the extensive random number generation capabilities in R (including the nlme and lme4 packages) to produce simulated data like you expect to collect and try to analyze the simulated data. You should simulate both what you expect to see and the null hypothesis as well. If you encounter difficulties doing that, please submit another question to this listserve. Before submitting another post, I suggest you help yourself by reading the posting guide! www.R-project.org/posting-guide.html. Anecdotal evidence suggests that posts that are more consistent with this posting guide generally get more useful replies quicker. bon chance. spencer graves [EMAIL PROTECTED] wrote: Dear R-users, We expect to develop statistic procedures and environnement for the computational analysis of our experimental datas. To provide a proof of concept, we plan to implement a test for a given experiment. Its design split data into 10 groups (including a control one) with 2 mesures for each (ref at t0 and response at t1). We aim to compare each group response with control response (group 1) using a multiple comparison procedure (Dunnett test). Before achieving this, we have to normalize our data : response values cannot be compared if base line isn't corrected. Covariance analysis seems to represent the best way to do this. But how to perform this by using R ? Actually, we have identify some R functions of interest regarding this matter (lme(), lm() and glm()). For example we plan to do as describe : glm(response~baseline) and then simtest(response_corrected~group, type=Dunnett, ttype=logical) If a mixed model seems to better fit our experiment, we have some problems on using the lme function : lme(response~baseline) returns an error (Invalid formula for groups). So : Are fitted values represent the corrected response ? Is it relevant to perform these tests in our design ? And how to use lme in a glm like way ? If someone could bring us your its precious knowledge to validate our analytical protocol and to express its point of view on implementation strategy ? Best regards. Alexandre MENICACCI Bioinformatics - FOURNIER PHARMA 50, rue de Dijon - 21121 Daix - FRANCE [EMAIL PROTECTED] tél : 03.80.44.76.17 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Tidal Time Series Analysis in R
I haven't seen any replies, so I will offer a few thoughts: If this were my problem, I'd start using lm with terms that predict the tide in terms of the positions of the sun and moon. If other things were thought to be important, I'd include them also in the lm model. However, I would not believe the answer from lm until the autocorrelation function (acf) of the residuals was plausibly white noise. If for some reason I couldn't get there, I'd use arima with its xreg argument to describe the plausible deterministic effects plus some ARMA model identified in the residuals. For documentation, I suggest you start with Venables and Ripley (2002) Modern Applied Statistics with S, 4th ed. (Springer), especially ch. 14. There may be better books on R available today, but among the books I've used, this is the premier general reference on R, and ch. 14 on time series should help you get started. I'm not familiar with hoa and Rwave. They may ultimately be the best tools for your need, but I don't know that. If it were my project, I would likely try lm and arima as discussed in Venables and Ripley first. Then I'd look at hoa, Rwave and anything else that might look promising. Should you desire other help from this listserve, I suggest you PLEASE do read the posting guide! www.R-project.org/posting-guide.html. Anecdotal evidence suggests that people who follow more closely that guide tend to get more useful replies quicker. hope this helps. spencer graves William H. Asquith wrote: I am looking at using R to analyze time series data containing a tidal component. I need to remove the tidal signal to extract the time series of the phenomena I seek to study. A browse of R-project search engines has not been too fruitful? I've found 'hoa' and 'Rwave', but need further help getting started. THANKS. -wa __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Accounting for within family correlation in genetic analysis
I am hoping for help with a genetic analysis. I am trying to perform an analysis of the relation between genes at a given locus (rs2304795) and a phenotypic trait (zerotg). Multiple subjects are recruited from each family (and so share a part of their genome and are correlated). Family groups are identified by the variable FAMILY. For each subject multiple measurements are made of the trait of interest (zerotg) over the course of several hours (time). The data within each subject are of course correlated and over time each subject's response is curvilinear. I know how to deal with the within subject correlation using a growth-curve analysis: zerotgQuad0ML.lme - lme (zerotg ~ time + time^2 + rs2304795 + rs2304795 * time + rs2304795 * time^2, data = GC, random = ~ 1 + time + time^2 | HAPI, na.action = na.omit, method = REML) but I don't know how I can modify the model to account for the within-family correlation. I would appreciate any suggestions for modifications to the code above that would allow me to account for the within-family correlation of observations. Thanks, John John Sorkin M.D., Ph.D. Chief, Biostatistics and Informatics Baltimore VA Medical Center GRECC and University of Maryland School of Medicine Claude Pepper OAIC University of Maryland School of Medicine Division of Gerontology Baltimore VA Medical Center 10 North Greene Street GRECC (BT/18/GR) Baltimore, MD 21201-1524 410-605-7119 - NOTE NEW EMAIL ADDRESS: [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
[R] read.table error
Hey, Once again I ask for some quick help. Here is some code: ovendata- read.table(ovens.dat,header=TRUE) attach(ovendata) print(ovendata) Here is the .dat file: DOne Two Three FourFiveSeven Eight 1130254 252 375 384 252 375 876 127 250 250 384 386 251 378 875 Here is the R Console output: ovendata- read.table(ovens.dat,header=TRUE) Warning message: incomplete final line found by readTableHeader on 'ovens.dat' attach(ovendata) The following object(s) are masked from ovendata ( position 3 ) : D Eight Five Four One Seven Three Two The following object(s) are masked from ovendata ( position 4 ) : D Eight Five Four One Seven Three Two The following object(s) are masked from ovendata ( position 5 ) : Eight Five Four One Seven Three Two The following object(s) are masked from package:stats : D print(ovendata) D One Two Three Four Five Seven Eight 1 1130 254 252 375 384 252 375 876 2 127 250 250 384 386 251 378 875 I've never seen anything like theis before. What's going on? Eric __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Warnings about user error (was read.table error)
I see no error here, let alone an error in read.table as claimed in your subject line. The posting guide does specifically ask `Use an informative subject line'. Please distinguish warnings about _your_ usage from errors in R. The first warning is that R fixed up an error in your file: it is missing a newline at the end of the last line (we can't see that in your listing). The remaining warnings come from attach() and say you have already repeatedly attach()ed ovendata. Learn to use detach() to match attach(). Also, in attaching ovendata you mask the function D in package stats, which is probably OK as you are not using it, and your D is a not a function. On Wed, 7 Dec 2005, Eric C. Jennings wrote: Hey, Once again I ask for some quick help. Here is some code: ovendata- read.table(ovens.dat,header=TRUE) attach(ovendata) print(ovendata) Here is the .dat file: DOne Two Three FourFiveSeven Eight 1130254 252 375 384 252 375 876 127 250 250 384 386 251 378 875 Here is the R Console output: ovendata- read.table(ovens.dat,header=TRUE) Warning message: incomplete final line found by readTableHeader on 'ovens.dat' attach(ovendata) The following object(s) are masked from ovendata ( position 3 ) : D Eight Five Four One Seven Three Two The following object(s) are masked from ovendata ( position 4 ) : D Eight Five Four One Seven Three Two The following object(s) are masked from ovendata ( position 5 ) : Eight Five Four One Seven Three Two The following object(s) are masked from package:stats : D print(ovendata) D One Two Three Four Five Seven Eight 1 1130 254 252 375 384 252 375 876 2 127 250 250 384 386 251 378 875 I've never seen anything like theis before. What's going on? Eric __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] mle.stepwise versus step/stepAIC
You could ask the author of the contributed package 'wle', which is not part of the `R-software'. The documentation is minimal, but the references are all to classical stepwise methods such as those in package leaps, that is they select columns in the model matrix and _not_ terms. On Wed, 7 Dec 2005, Jim Brindle wrote: Hello, I have a question pertaining to the stepwise regression which I am trying to perform. I have a data set in which I have 14 predictor variables accompanying my response variable. I am not sure what the difference is between the function mle.stepwise found in the wle package and the functions step or stepAIC? When would one use mle.stepwise versus step/stepAIC? Other than the references in the R-software, is there anything which discusses the use of mle.stepwise? Thank you kindly in advance for any insight offered. Jim Brindle -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R is GNU S, not C.... [was how to get or store .....]
DeepS == Deepayan Sarkar [EMAIL PROTECTED]
on Wed, 7 Dec 2005 19:15:52 -0600 writes:
DeepS On 12/7/05, Jan T. Kim [EMAIL PROTECTED] wrote:
On Tue, Dec 06, 2005 at 04:21:01PM +, Patrick Burns wrote:
I don't put in extraneous ';' because I maybe get a
blister on my little finger.
I suspect that those who find the semi-colons ugly in
R do not find them ugly in C. I think the reason there
would be a visceral reaction in R but not in C is that
there is a danger when using them in R that they really
mean something.
We get questions on R-help often enough about why
code like:
if(x 0) y - 4
else y - 4.5e23
doesn't work.
If people habitually used semi-colons, those sorts of
questions would probably multiply.
I don't understand. It would seem to me that in
if (x 0) y - 4;
else y - 4.5e23;
it's pretty obvious that the if statement is terminated by the
semicolon at the end of the first line and that therefore, the else
on the next line is erroneous because it is not associated with any
if.
DeepS Why is it obvious?
DeepSif (x 0) y = 4;
DeepSelse y = 4.5e23;
DeepS is perfectly legal C.
Very good! I'm glad we got here, namely back to the original
subject, and now have seen a nice example of why it can be a bad
idea to maltreat the S language by extraneous ; ..
At least, the version above fails consistently, i.e. regardless of
context.
DeepS It fails for unrelated reasons.
DeepS if (x 0) { y - 4; } else { y - 4.5e23; }
DeepS doesn't fail.
DeepS -Deepayan
Thank you, Deepayan, Brian,
and all the others in the knowning who hadn't lost patience with
this thread yet... ;-)
Martin
On the other hand, I've studied the R Language Definition for
quite some time before fully understanding why
if (x 0) y - 4
else y - 4.5e23
works inside of a function (or other enclosing block) while it does not
work interactively (or at the top level of a script).
Best regards, Jan
--
+- Jan T. Kim ---+
| email: [EMAIL PROTECTED] |
| WWW: http://www.cmp.uea.ac.uk/people/jtk |
*-= hierarchical systems are for files, not for humans =-*
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
