Re: [R] String in data frame
Please read `An Introduction to R': R does not have `strings' or `namelists' whereas it does have `lists' but d$num is not one. To communicate, you need a common language with your readers. You are not speaking R, and your `examples' are not R output. I suspect d$num is a numeric vector and d$name is a factor. Please look up those concepts and the help for c(), noting that `NA' is a logical vector. In particular you need to understand The default method combines its arguments to form a vector. All arguments are coerced to a common type which is the type of the returned value. On Fri, 15 Apr 2005, Cuichang Zhao wrote: hello, how can take the string in the data frame. right now i have a table that create as a data frame and stored in the file called data.xls and now i want to read data frame as a table in my another r program, i used the following command: the first column of the data frame is just one number called num, but the second one a list of string, called name. d - read.table(data.xls, header = T); what i get for d is a table of 2 columns: num and name: the one for d$num is just a normal list without any levels with it, and this is what i want. but the second d$name is a list with a levels with it. the list i have for d$name is: d$name = (bal bal bal bal bal bal), levlels:bal, however, when i want to have for following code, something different happens: namelist - NA; a - d$name[1]; #this will outputs a = bal, lelvels:bal namelist - c(namelist, a); #this does not outptu (NA, bal), instead it outputs (NA, 1), if i keep #adding a to the namelist, it keeps adding 1 to the namelist instead of bal. However, i want to add bal to the namelist, not 1, so how i can do this? [[alternative HTML version deleted]] PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html PLEASE do as we ask (and not sent HTML mail). -- 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] example on front page doesn't work in R 2.0.1
Hubert Feyrer wrote: On Fri, 15 Apr 2005, Sundar Dorai-Raj wrote: You snipped too much: ... This is where it says ade4, etc. is not found. If you would like to install these packages, it's rather easy: install.packages(c(ade4, RColorBrewer, pixmap)) Ah, thanks! i'll try to see how I can add the packages - using Unix, I cannot write to the R package's directory. I'll see how to work around this, but thank you very much for the help! Read the docs how to create your own library of packages. It's easy! Uwe Ligges Maybe it would be good to not require any special packages for an example offered on the frontpage of www.r-project.org though... :) - Hubert __ R-help@stat.math.ethz.ch 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] Pearson corelation and p-value for matrix
Dear Andy,
That's clearly much better -- and illustrates an effective strategy for
vectorizing (or matricizing) a computation. I think I'll add this to my
list of programming examples. It might be a little dangerous to pass ...
through to cor(), since someone could specify type=spearman, for example.
Thanks,
John
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
-Original Message-
From: Liaw, Andy [mailto:[EMAIL PROTECTED]
Sent: Friday, April 15, 2005 9:51 PM
To: 'John Fox'; [EMAIL PROTECTED]
Cc: 'R-Help'; 'Dren Scott'
Subject: RE: [R] Pearson corelation and p-value for matrix
We can be a bit sneaky and `borrow' code from cor.test.default:
cor.pval - function(x, alternative=two-sided, ...) {
corMat - cor(x, ...)
n - nrow(x)
df - n - 2
STATISTIC - sqrt(df) * corMat / sqrt(1 - corMat^2)
p - pt(STATISTIC, df)
p - if (alternative == less) {
p
} else if (alternative == greater) {
1 - p
} else 2 * pmin(p, 1 - p)
p
}
The test:
system.time(c1 - cor.pvals(X), gcFirst=TRUE)
[1] 13.19 0.01 13.58NANA
system.time(c2 - cor.pvalues(X), gcFirst=TRUE)
[1] 6.22 0.00 6.42 NA NA
system.time(c3 - cor.pval(X), gcFirst=TRUE)
[1] 0.07 0.00 0.07 NA NA
Cheers,
Andy
From: John Fox
Dear Mark,
I think that the reflex of trying to avoid loops in R is often
mistaken, and so I decided to try to time the two approaches (on a
3GHz Windows XP system).
I discovered, first, that there is a bug in your function -- you
appear to have indexed rows instead of columns; fixing that:
cor.pvals - function(mat)
{
cols - expand.grid(1:ncol(mat), 1:ncol(mat))
matrix(apply(cols, 1,
function(x) cor.test(mat[, x[1]], mat[,
x[2]])$p.value),
ncol = ncol(mat))
}
My function is cor.pvalues and yours cor.pvals. This is for a data
matrix with 1000 observations on 100 variables:
R - diag(100)
R[upper.tri(R)] - R[lower.tri(R)] - .5
library(mvtnorm)
X - rmvnorm(1000, sigma=R)
dim(X)
[1] 1000 100
system.time(cor.pvalues(X))
[1] 5.53 0.00 5.53 NA NA
system.time(cor.pvals(X))
[1] 12.66 0.00 12.66NANA
I frankly didn't expect the advantage of my approach to be
this large,
but there it is.
Regards,
John
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
-Original Message-
From: Marc Schwartz [mailto:[EMAIL PROTECTED]
Sent: Friday, April 15, 2005 7:08 PM
To: John Fox
Cc: 'Dren Scott'; R-Help
Subject: RE: [R] Pearson corelation and p-value for matrix
Here is what might be a slightly more efficient way to
get to John's
question:
cor.pvals - function(mat)
{
rows - expand.grid(1:nrow(mat), 1:nrow(mat))
matrix(apply(rows, 1,
function(x) cor.test(mat[x[1], ], mat[x[2],
])$p.value),
ncol = nrow(mat))
}
HTH,
Marc Schwartz
On Fri, 2005-04-15 at 18:26 -0400, John Fox wrote:
Dear Dren,
How about the following?
cor.pvalues - function(X){
nc - ncol(X)
res - matrix(0, nc, nc)
for (i in 2:nc){
for (j in 1:(i - 1)){
res[i, j] - res[j, i] - cor.test(X[,i],
X[,j])$p.value
}
}
res
}
What one then does with all of those non-independent test
is another
question, I guess.
I hope this helps,
John
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of
Dren Scott
Sent: Friday, April 15, 2005 4:33 PM
To: r-help@stat.math.ethz.ch
Subject: [R] Pearson corelation and p-value for matrix
Hi,
I was trying to evaluate the pearson correlation and
the p-values
for an nxm matrix, where each row represents a vector.
One way to do
it would be to iterate through each row, and find its
correlation
value( and the p-value) with respect to the other rows.
Is there
some function by which I can use the matrix as input?
Ideally, the
output would be an nxn matrix, containing the p-values
between the
respective vectors.
I have tried cor.test for the iterations, but couldn't find a
function that would take the matrix as input.
Thanks for the help.
Dren
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
Re: [R] Getting subsets of a data frame
Thanks, it's interesting reading. I also noticed that sw[, 1, drop = TRUE] is a vector (coerces to the lowest dimension) but sw[1, , drop = TRUE] is a one-row data frame (does not convert it into a list or vector) FS On 4/16/05, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote: You should look at ?[ and look very carefully at the drop argument. For your example sw[, 1] is the first component of the data frame, but sw[, 1, drop = FALSE] is a data frame consisting of just the first component, as mathematically fastidious people would expect. This is a convention, and like most arbitrary conventions it can be very useful most of the time, but some of the time it can be a very nasty trap. Caveat emptor. Bill Venables. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Fernando Saldanha Sent: Saturday, 16 April 2005 1:07 PM To: Submissions to R help Subject: [R] Getting subsets of a data frame I was reading in the Reference Manual about Extract.data.frame. There is a list of examples of expressions using [ and [[, with the outcomes. I was puzzled by the fact that, if sw is a data frame, then sw[, 1:3] is also a data frame, but sw[, 1] is just a vector. Since R has no scalars, it must be the case that 1 and 1:1 are the same: 1 == 1:1 [1] TRUE Then why isn't sw[,1] = sw[, 1:1] a data frame? FS __ R-help@stat.math.ethz.ch 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] Pearson corelation and p-value for matrix
Dear Mark,
-Original Message-
From: Marc Schwartz [mailto:[EMAIL PROTECTED]
Sent: Friday, April 15, 2005 9:41 PM
To: John Fox
Cc: 'R-Help'; 'Dren Scott'
Subject: RE: [R] Pearson corelation and p-value for matrix
John,
Interesting test. Thanks for pointing that out.
You are right, there is a knee-jerk reaction to avoid loops,
especially nested loops.
On the indexing of rows, I did that because Dren had
indicated in his initial post:
I was trying to evaluate the pearson correlation and the p-values
for an nxm matrix, where each row represents a vector.
One way to do it would be to iterate through each row, and find its
correlation value( and the p-value) with respect to the other rows.
So I ran the correlations by row, rather than by column.
That's the second time yesterday that I responded to a posting without
reading it carefully enough -- a good lesson for me. I guess that Dren could
just apply my solution to the transpose of his matrix -- i.e.,
cor.pvalues(t(X)).
Sorry,
John
Thanks again. Good lesson.
Marc
On Fri, 2005-04-15 at 21:36 -0400, John Fox wrote:
Dear Mark,
I think that the reflex of trying to avoid loops in R is often
mistaken, and so I decided to try to time the two approaches (on a
3GHz Windows XP system).
I discovered, first, that there is a bug in your function -- you
appear to have indexed rows instead of columns; fixing that:
cor.pvals - function(mat)
{
cols - expand.grid(1:ncol(mat), 1:ncol(mat))
matrix(apply(cols, 1,
function(x) cor.test(mat[, x[1]], mat[,
x[2]])$p.value),
ncol = ncol(mat))
}
My function is cor.pvalues and yours cor.pvals. This is for a data
matrix with 1000 observations on 100 variables:
R - diag(100)
R[upper.tri(R)] - R[lower.tri(R)] - .5
library(mvtnorm)
X - rmvnorm(1000, sigma=R)
dim(X)
[1] 1000 100
system.time(cor.pvalues(X))
[1] 5.53 0.00 5.53 NA NA
system.time(cor.pvals(X))
[1] 12.66 0.00 12.66NANA
I frankly didn't expect the advantage of my approach to be
this large,
but there it is.
Regards,
John
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
-Original Message-
From: Marc Schwartz [mailto:[EMAIL PROTECTED]
Sent: Friday, April 15, 2005 7:08 PM
To: John Fox
Cc: 'Dren Scott'; R-Help
Subject: RE: [R] Pearson corelation and p-value for matrix
Here is what might be a slightly more efficient way to
get to John's
question:
cor.pvals - function(mat)
{
rows - expand.grid(1:nrow(mat), 1:nrow(mat))
matrix(apply(rows, 1,
function(x) cor.test(mat[x[1], ], mat[x[2],
])$p.value),
ncol = nrow(mat))
}
HTH,
Marc Schwartz
On Fri, 2005-04-15 at 18:26 -0400, John Fox wrote:
Dear Dren,
How about the following?
cor.pvalues - function(X){
nc - ncol(X)
res - matrix(0, nc, nc)
for (i in 2:nc){
for (j in 1:(i - 1)){
res[i, j] - res[j, i] - cor.test(X[,i],
X[,j])$p.value
}
}
res
}
What one then does with all of those non-independent test
is another
question, I guess.
I hope this helps,
John
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Dren
Scott
Sent: Friday, April 15, 2005 4:33 PM
To: r-help@stat.math.ethz.ch
Subject: [R] Pearson corelation and p-value for matrix
Hi,
I was trying to evaluate the pearson correlation and the
p-values for an nxm matrix, where each row represents
a vector.
One way to do
it would be to iterate through each row, and find its
correlation value( and the p-value) with respect to the other
rows. Is there some function by which I can use the
matrix as input?
Ideally, the
output would be an nxn matrix, containing the p-values
between the
respective vectors.
I have tried cor.test for the iterations, but couldn't find a
function that would take the matrix as input.
Thanks for the help.
Dren
__
R-help@stat.math.ethz.ch 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] function corresponding to map of perl
I defined map as follows:
map - function(x, y, fun) {
mymat - matrix( c(x,y), c(length(x), 2) )
tmat - t(mymat)
oldmat - tmat
result - apply(tmat, 2, function(x) {fun(x[1], x[2])})
}
It seems to work (see below). Of course you can turn it into a one-liner.
a-c(1,2,3)
b-c(4,5,6)
mysum - function(x, y) {x + y}
map - function(x, y, fun) {
+ mymat - matrix( c(x,y), c(length(x), 2) )
+ tmat - t(mymat)
+ oldmat - tmat
+ result - apply(tmat, 2, function(x) {fun(x[1], x[2])})
+ }
(test - map(a, b, mysum))
[1] 5 7 9
__
R-help@stat.math.ethz.ch 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] Getting subsets of a data frame
Because a data frame can hold different data types (even matrices) in different variables, one row of it can not be converted to a vector in general (where all elements need to be of the same type). Andy From: Fernando Saldanha Thanks, it's interesting reading. I also noticed that sw[, 1, drop = TRUE] is a vector (coerces to the lowest dimension) but sw[1, , drop = TRUE] is a one-row data frame (does not convert it into a list or vector) FS On 4/16/05, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote: You should look at ?[ and look very carefully at the drop argument. For your example sw[, 1] is the first component of the data frame, but sw[, 1, drop = FALSE] is a data frame consisting of just the first component, as mathematically fastidious people would expect. This is a convention, and like most arbitrary conventions it can be very useful most of the time, but some of the time it can be a very nasty trap. Caveat emptor. Bill Venables. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Fernando Saldanha Sent: Saturday, 16 April 2005 1:07 PM To: Submissions to R help Subject: [R] Getting subsets of a data frame I was reading in the Reference Manual about Extract.data.frame. There is a list of examples of expressions using [ and [[, with the outcomes. I was puzzled by the fact that, if sw is a data frame, then sw[, 1:3] is also a data frame, but sw[, 1] is just a vector. Since R has no scalars, it must be the case that 1 and 1:1 are the same: 1 == 1:1 [1] TRUE Then why isn't sw[,1] = sw[, 1:1] a data frame? FS __ R-help@stat.math.ethz.ch 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] Pearson corelation and p-value for matrix
From: John Fox
Dear Andy,
That's clearly much better -- and illustrates an effective
strategy for
vectorizing (or matricizing) a computation. I think I'll
add this to my
list of programming examples. It might be a little dangerous
to pass ...
through to cor(), since someone could specify
type=spearman, for example.
Ah, yes, the ... isn't likely to help here! Also, it will
only work correctly if there are no NA's, for example (or
else the degree of freedom would be wrong).
Best,
Andy
Thanks,
John
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
-Original Message-
From: Liaw, Andy [mailto:[EMAIL PROTECTED]
Sent: Friday, April 15, 2005 9:51 PM
To: 'John Fox'; [EMAIL PROTECTED]
Cc: 'R-Help'; 'Dren Scott'
Subject: RE: [R] Pearson corelation and p-value for matrix
We can be a bit sneaky and `borrow' code from cor.test.default:
cor.pval - function(x, alternative=two-sided, ...) {
corMat - cor(x, ...)
n - nrow(x)
df - n - 2
STATISTIC - sqrt(df) * corMat / sqrt(1 - corMat^2)
p - pt(STATISTIC, df)
p - if (alternative == less) {
p
} else if (alternative == greater) {
1 - p
} else 2 * pmin(p, 1 - p)
p
}
The test:
system.time(c1 - cor.pvals(X), gcFirst=TRUE)
[1] 13.19 0.01 13.58NANA
system.time(c2 - cor.pvalues(X), gcFirst=TRUE)
[1] 6.22 0.00 6.42 NA NA
system.time(c3 - cor.pval(X), gcFirst=TRUE)
[1] 0.07 0.00 0.07 NA NA
Cheers,
Andy
From: John Fox
Dear Mark,
I think that the reflex of trying to avoid loops in R is often
mistaken, and so I decided to try to time the two
approaches (on a
3GHz Windows XP system).
I discovered, first, that there is a bug in your function -- you
appear to have indexed rows instead of columns; fixing that:
cor.pvals - function(mat)
{
cols - expand.grid(1:ncol(mat), 1:ncol(mat))
matrix(apply(cols, 1,
function(x) cor.test(mat[, x[1]], mat[,
x[2]])$p.value),
ncol = ncol(mat))
}
My function is cor.pvalues and yours cor.pvals. This is
for a data
matrix with 1000 observations on 100 variables:
R - diag(100)
R[upper.tri(R)] - R[lower.tri(R)] - .5
library(mvtnorm)
X - rmvnorm(1000, sigma=R)
dim(X)
[1] 1000 100
system.time(cor.pvalues(X))
[1] 5.53 0.00 5.53 NA NA
system.time(cor.pvals(X))
[1] 12.66 0.00 12.66NANA
I frankly didn't expect the advantage of my approach to be
this large,
but there it is.
Regards,
John
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
-Original Message-
From: Marc Schwartz [mailto:[EMAIL PROTECTED]
Sent: Friday, April 15, 2005 7:08 PM
To: John Fox
Cc: 'Dren Scott'; R-Help
Subject: RE: [R] Pearson corelation and p-value for matrix
Here is what might be a slightly more efficient way to
get to John's
question:
cor.pvals - function(mat)
{
rows - expand.grid(1:nrow(mat), 1:nrow(mat))
matrix(apply(rows, 1,
function(x) cor.test(mat[x[1], ], mat[x[2],
])$p.value),
ncol = nrow(mat))
}
HTH,
Marc Schwartz
On Fri, 2005-04-15 at 18:26 -0400, John Fox wrote:
Dear Dren,
How about the following?
cor.pvalues - function(X){
nc - ncol(X)
res - matrix(0, nc, nc)
for (i in 2:nc){
for (j in 1:(i - 1)){
res[i, j] - res[j, i] - cor.test(X[,i],
X[,j])$p.value
}
}
res
}
What one then does with all of those non-independent test
is another
question, I guess.
I hope this helps,
John
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of
Dren Scott
Sent: Friday, April 15, 2005 4:33 PM
To: r-help@stat.math.ethz.ch
Subject: [R] Pearson corelation and p-value for matrix
Hi,
I was trying to evaluate the pearson correlation and
the p-values
for an nxm matrix, where each row represents a vector.
One way to do
it would be to iterate through each row, and find its
correlation
value( and the p-value) with respect to the other rows.
Is there
some function by which I can use the matrix as input?
Ideally, the
output would be an nxn matrix, containing the p-values
between the
respective vectors.
Re: [R] How to get predictions, plots, etc. from lmer{lme4}
Michael Kubovy wrote: Kindly send a cc to me when replying to the list. I'm having trouble using lmer beyond a first step. My data: some(exp1B) sub ba amplitude a b c d 2 1 1.00 1.5 65 63 4 8 414 1.15 0.0 92 41 3 4 434 1.15 3.0 88 48 2 2 636 1.00 3.0 50 72 9 9 778 1.15 0.0 112 25 2 1 89 10 1.15 0.0 37 33 36 34 126 13 1.15 1.5 80 50 6 4 140 14 1.15 4.5 12 115 6 7 145 20 1.00 0.0 73 65 0 2 147 20 1.00 3.0 63 72 2 3 etc. The output: summary(exp1B.both.cont.lmer) Linear mixed-effects model fit by maximum likelihood Formula: cbind(b, a) ~ ba + amplitude + (1 | sub) Data: exp1B AIC BIClogLik MLdeviance REMLdeviance 768.0086 783.2579 -379.0043 758.0086 766.3104 Random effects: Groups NameVariance Std.Dev. sub (Intercept) 0.18787 0.43344 Residual 6.36355 2.52261 # of obs: 156, groups: sub, 13 Fixed effects: Estimate Std. Error DF t value Pr(|t|) (Intercept) 4.118806 0.397780 153 10.3545 2.2e-16 ba -4.205518 0.330010 153 -12.7436 2.2e-16 amplitude 0.137958 0.023754 153 5.8076 3.547e-08 This is just what I need. But I also need predicted values, plots, etc, and can't figure out how to proceed. Have I overlooked a more extended document than the rather terse (for me at least) help page? I regret to say no. Methods for generalized linear mixed models fit by lmer are still rather rudimentary. __ R-help@stat.math.ethz.ch 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] Almost seasonal decomposition
Dear ts-friends, I have an almost seasonal signal. It's human respiration pressure, where the respiration signal is regular but non-sinusoidal, with a slightly drifting period. Overlaid on the signal are non-periodic twitches (think weak hickups) with amplitudes 2-5 times higher than the respiration signal. We have used adaptive filtering in a similar case before, but no twitch-free reference is available in the present case. I want to do something like stl, allowing some freedom to drift for the seasonal component. Robust processing is a must, the hickup should not affect processing. Robust autoregressive processing in package wle looked like a good choice, but I could not get the model to converge with my data set. Any ideas Dieter Menne __ R-help@stat.math.ethz.ch 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] can test the if relationship is significant in cancor?
ronggui wrote:
i have try hard to find the answer by google,but i can not find any solution.
so i wan to ask:
1,can we test the if canonical relationship is significant after using cancor?
One reference is T. W. Anderson: An Introduction to Multivariate
Statistical Analysis, second edition, pages 497-498.
2,if it can,how?
Following the reference above:
cancor.test - function(obj, N){
# obj is object returned from cancor
# N is sample size, which is not contained in the cancor object!
p1 - NROW(obj$xcoef)
p2 - NROW(obj$ycoef)
p - p1 + p2
r - length(obj$cor)
# Calculating Bartlett modification of minus twice log likelihood:
bartlett - -(N-0.5*(p+3))*sum( log( 1-obj$cor^2))
# which is approximately chi-squared with p1p2 degrees of freedom:
list(bartlett=bartlett, p.value=pchisq(bartlett, df=p1*p2,
lower.tail=FALSE))
}
This tests if ALLl the canonical correlations are zero. Anybody knows
how good this approximation is,
and how dependent on multivariate normality?
Kjetil
3,if not,is it under-developed or there is not need to do it?or there is no
good way to do it?
i hope my question is not too silly.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
--
Kjetil Halvorsen.
Peace is the most effective weapon of mass construction.
-- Mahdi Elmandjra
__
R-help@stat.math.ethz.ch 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] Define local function
Peter Dalgaard wrote:
Ales Ziberna [EMAIL PROTECTED] writes:
I am also very interested how this could be done, possibly in such a
way that this would be incorporated in the function itself and there
wouldn't be a need to write environment(f) - NULL before calling a
function, as is proposed in the reply below and in a thread a few days
ago!
Why worry about where the line occurs? In R, there is little
distinction between functions and data, so you could say that the
function definition isn't done until you're finished modifying it.
Notice BTW, that environment(f) - NULL may have unexpected
consequences. What it really means is that the lexical scope of f
becomes the base package. This interpretation of NULL may change in
the future, since it is somewhat illogical and it has a couple of
undesirable consequences that there's no way to specify a truly empty
environment. So
a) if you're calling a function outside of the base package, you get
the effect of
f - function(){mean(rnorm(10))}
environment(f)-NULL
f()
Error in mean(rnorm(10)) : couldn't find function rnorm
b) even if it does work now, it may be broken by a future change to R.
Notice that *all* functions contain unbound variables in the form of
functions so if we get an empty NULL environment, even - may stop
working.
A more likely thing to want to do is to see the search path only after
the global environment, which you can do with
environment(f) - as.environment(2)
(the 2 says to use the 2nd environment on the search path). This might
not be perfect, since you can use attach() to modify the search path,
but it would likely be good enough to flush out simple programming bugs.
Duncan Murdoch
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to get predictions, plots, etc. from lmer{lme4}
Any comments on the following strategy: (1) Take the log-odds of the frequencies of b and a (adding 1/6, as Tukey recommended) and run lme. (2) If the estimates from lme are in line with the estimates I got from lmer, then use the results from lme. On Apr 16, 2005, at 9:16 AM, Douglas Bates wrote: Michael Kubovy wrote: I'm having trouble using lmer beyond a first step. My data: some(exp1B) sub ba amplitude a b c d 2 1 1.00 1.5 65 63 4 8 414 1.15 0.0 92 41 3 4 434 1.15 3.0 88 48 2 2 636 1.00 3.0 50 72 9 9 778 1.15 0.0 112 25 2 1 89 10 1.15 0.0 37 33 36 34 126 13 1.15 1.5 80 50 6 4 140 14 1.15 4.5 12 115 6 7 145 20 1.00 0.0 73 65 0 2 147 20 1.00 3.0 63 72 2 3 etc. The output: summary(exp1B.both.cont.lmer) Linear mixed-effects model fit by maximum likelihood Formula: cbind(b, a) ~ ba + amplitude + (1 | sub) Data: exp1B AIC BIClogLik MLdeviance REMLdeviance 768.0086 783.2579 -379.0043 758.0086 766.3104 Random effects: Groups NameVariance Std.Dev. sub (Intercept) 0.18787 0.43344 Residual 6.36355 2.52261 # of obs: 156, groups: sub, 13 Fixed effects: Estimate Std. Error DF t value Pr(|t|) (Intercept) 4.118806 0.397780 153 10.3545 2.2e-16 ba -4.205518 0.330010 153 -12.7436 2.2e-16 amplitude 0.137958 0.023754 153 5.8076 3.547e-08 This is just what I need. But I also need predicted values, plots, etc, and can't figure out how to proceed. Have I overlooked a more extended document than the rather terse (for me at least) help page? I regret to say no. Methods for generalized linear mixed models fit by lmer are still rather rudimentary. _ Professor Michael Kubovy University of Virginia Department of Psychology USPS: P.O.Box 400400 Charlottesville, VA 22904-4400 Parcels:Room 102Gilmer Hall McCormick Road Charlottesville, VA 22903 Office: B011+1-434-982-4729 Lab:B019+1-434-982-4751 Fax:+1-434-982-4766 WWW:http://www.people.virginia.edu/~mk9y/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] [R-pkgs] g.data version 1.6, upgrade in response to changes in R-2.1.0
Version 1.6 of the g.data package is available on CRAN. The g.data package is used to create and maintain delayed-data packages (DDP's). Data stored in a DDP are available on demand, but do not take up memory until requested. You attach a DDP with g.data.attach(), then read from it and assign to it in a manner similar to S-Plus, except that you must run g.data.save() to actually commit to disk. This new version uses delayedAssign (new to R-2.1.0) instead of delay (deprecated in R-2.1.0). If it detects you are using R-2.1.0 or higher, it will automatically upgrade older DDP's. Version 1.6 of g.data will still work fine under earlier versions of R too. Note, however, that once you upgrade an older DDP (by g.data.attach'ing it under R-2.1.0), it will no longer be readable by earlier versions of R (which do not understand delayedAssign), so users running multiple versions of R may want to avoid this upgrade. Please contact me if this presents a problem. -- David Brahm ([EMAIL PROTECTED]) [[alternative HTML version deleted]] ___ R-packages mailing list [EMAIL PROTECTED] https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] [R-pkgs] bayesm: a package for Bayesian infererence for Marketing/Micro-Econometrics
We are pleased to announce the release of version 0.0 of bayesm on CRAN. bayesm covers many important models used in marketing and micro-econometrics applications. The package includes: Bayes Regression (univariate or multivariate dep var) Multinomial Logit Multinomial and Multivariate Probit Multivariate Mixtures of Normals Hierarchical Linear Models with a normal prior and covariates Hierarchical Multinomial Logits with mixture of normals prior Bayesian analysis of choice-based conjoint data Bayesian treatment of linear instrumental variables models Analyis of Multivariate Ordinal survey data with scale usage heterogeneity after installing the package, use help.search(mcmc) to display the key functions. bayem implements the models discussed in our book, Bayesian Statistics and Marketing To view a draft of the book visit http://gsbwww.uchicago.edu/fac/peter.rossi/research/bsm.html We would very much appreciate any comments/errors/suggestions for future development. Peter E. Rossi Joseph T. Lewis Professor of Marketing and Statistics Editor, Quantitative Marketing and Economics Rm 360, Graduate School of Business, U of Chicago 5807 S. Woodlawn Ave, Chicago IL 60637 Tel: (773) 702-7513 | Fax: (773) 834-2081 [EMAIL PROTECTED] WWW: http://ChicagoGsb.edu/fac/peter.rossi SSRN: http://ssrn.com/author=22862 QME: http://www.kluweronline.com/issn/1570-7156 ___ R-packages mailing list [EMAIL PROTECTED] https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] function corresponding to map of perl
From: Fernando Saldanha
I defined map as follows:
map - function(x, y, fun) {
mymat - matrix( c(x,y), c(length(x), 2) )
tmat - t(mymat)
oldmat - tmat
result - apply(tmat, 2, function(x) {fun(x[1], x[2])})
}
It seems to work (see below). Of course you can turn it into
a one-liner.
a-c(1,2,3)
b-c(4,5,6)
mysum - function(x, y) {x + y}
map - function(x, y, fun) {
+ mymat - matrix( c(x,y), c(length(x), 2) )
+ tmat - t(mymat)
+ oldmat - tmat
+ result - apply(tmat, 2, function(x) {fun(x[1], x[2])})
+ }
(test - map(a, b, mysum))
[1] 5 7 9
Maybe you're re-inventing mapply()?
mapply(mysum, a, b)
[1] 5 7 9
Andy
__
R-help@stat.math.ethz.ch 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] HOWTO compare univariate binomial glm lrm models which are not nested
Thanks a lot for the input! I forgot to add family=binomial, for a binomial glm. Now the AIC's are positive! I was planning to use AIC (from the binomial glm) and c-index (lrm) to compare and rank different uni-variate (one continue explanatory variable) logistic models to evaluate the 'performance' of the different explanatory variables in the different models. What is the best technique to compare these lrm.models, which are not nested? I found in literature that ranking based on different parameters (goodness of fit and predictability) that these can be used to compare uni-variate models. Thanks in advance, Regards, Jan- ___ ir. Jan Verbesselt Research Associate Lab of Geomatics Engineering K.U. Leuven Vital Decosterstraat 102. B-3000 Leuven Belgium Tel: +32-16-329750 Fax: +32-16-329760 http://gloveg.kuleuven.ac.be/ ___ -Original Message- From: Prof Brian Ripley [mailto:[EMAIL PROTECTED] Sent: Friday, April 15, 2005 5:06 PM To: Jan Verbesselt Cc: r-help@stat.math.ethz.ch Subject: Re: [R] negetative AIC values: How to compare models with negative AIC's AICs (like log-likelihoods) can be positive or negative. However, you fitted a Gaussian and not a binomial glm (as lrm does if m.arson is binary). For a discrete response with the usual dominating measure (counting measure) the log-likelihood is negative and hence the AIC is positive, but not in general (and it is matter of convention even there). In any case, Akaike only suggested comparing AIC for nested models, no one suggests comparing continuous and discrete models. On Fri, 15 Apr 2005, Jan Verbesselt wrote: Dear, When fitting the following model knots - 5 lrm.NDWI - lrm(m.arson ~ rcs(NDWI,knots) I obtain the following result: Logistic Regression Model lrm(formula = m.arson ~ rcs(NDWI, knots)) Frequencies of Responses 0 1 666 35 Obs Max Deriv Model L.R. d.f. P C Dxy Gamma Tau-a R2 Brier 701 5e-07 34.49 4 0 0.777 0.553 0.563 0.053 0.147 0.045 Coef S.E.Wald Z P Intercept -4.627 3.188 -1.45 0.1467 NDWI 5.333 20.724 0.26 0.7969 NDWI'6.832 74.201 0.09 0.9266 NDWI'' 10.469 183.915 0.06 0.9546 NDWI''' -190.566 254.590 -0.75 0.4541 When analysing the glm fit of the same model Call: glm(formula = m.arson ~ rcs(NDWI, knots), x = T, y = T) Coefficients: (Intercept) rcs(NDWI, knots)NDWIrcs(NDWI, knots)NDWI' rcs(NDWI, knots)NDWI'' rcs(NDWI, knots)NDWI''' 0.02067 0.08441 -0.54307 3.99550-17.38573 Degrees of Freedom: 700 Total (i.e. Null); 696 Residual Null Deviance: 33.25 Residual Deviance: 31.76AIC: -167.7 A negative AIC occurs! How can the negative AIC from different models be compared with each other? Is this result logical? Is the lowest AIC still correct? -- 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
[R] Re: HOWTO compare univariate binomial glm lrm models which are not nested
Compare them by `goodness for purpose': you have not told us the purpose. Please do read some of the extensive literature on model comparison. On Sat, 16 Apr 2005, Jan Verbesselt wrote: Thanks a lot for the input! I forgot to add family=binomial, for a binomial glm. Now the AIC's are positive! I was planning to use AIC (from the binomial glm) and c-index (lrm) to compare and rank different uni-variate (one continue explanatory variable) logistic models to evaluate the 'performance' of the different explanatory variables in the different models. What is the best technique to compare these lrm.models, which are not nested? I found in literature that ranking based on different parameters (goodness of fit and predictability) that these can be used to compare uni-variate models. Thanks in advance, Regards, Jan- ___ ir. Jan Verbesselt Research Associate Lab of Geomatics Engineering K.U. Leuven Vital Decosterstraat 102. B-3000 Leuven Belgium Tel: +32-16-329750 Fax: +32-16-329760 http://gloveg.kuleuven.ac.be/ ___ -Original Message- From: Prof Brian Ripley [mailto:[EMAIL PROTECTED] Sent: Friday, April 15, 2005 5:06 PM To: Jan Verbesselt Cc: r-help@stat.math.ethz.ch Subject: Re: [R] negetative AIC values: How to compare models with negative AIC's AICs (like log-likelihoods) can be positive or negative. However, you fitted a Gaussian and not a binomial glm (as lrm does if m.arson is binary). For a discrete response with the usual dominating measure (counting measure) the log-likelihood is negative and hence the AIC is positive, but not in general (and it is matter of convention even there). In any case, Akaike only suggested comparing AIC for nested models, no one suggests comparing continuous and discrete models. On Fri, 15 Apr 2005, Jan Verbesselt wrote: Dear, When fitting the following model knots - 5 lrm.NDWI - lrm(m.arson ~ rcs(NDWI,knots) I obtain the following result: Logistic Regression Model lrm(formula = m.arson ~ rcs(NDWI, knots)) Frequencies of Responses 0 1 666 35 Obs Max Deriv Model L.R. d.f. P C Dxy Gamma Tau-a R2 Brier 701 5e-07 34.49 4 0 0.777 0.553 0.563 0.053 0.147 0.045 Coef S.E.Wald Z P Intercept -4.627 3.188 -1.45 0.1467 NDWI 5.333 20.724 0.26 0.7969 NDWI'6.832 74.201 0.09 0.9266 NDWI'' 10.469 183.915 0.06 0.9546 NDWI''' -190.566 254.590 -0.75 0.4541 When analysing the glm fit of the same model Call: glm(formula = m.arson ~ rcs(NDWI, knots), x = T, y = T) Coefficients: (Intercept) rcs(NDWI, knots)NDWIrcs(NDWI, knots)NDWI' rcs(NDWI, knots)NDWI'' rcs(NDWI, knots)NDWI''' 0.02067 0.08441 -0.54307 3.99550-17.38573 Degrees of Freedom: 700 Total (i.e. Null); 696 Residual Null Deviance: 33.25 Residual Deviance: 31.76AIC: -167.7 A negative AIC occurs! How can the negative AIC from different models be compared with each other? Is this result logical? Is the lowest AIC still correct? -- 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 -- 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] Getting subsets of a data frame
Perhaps Fernando will also note that is documented in ?[.data.frame, a slightly more appropriate reference than Bill's. It would be a good idea to read a good account of R's indexing: Bill Venables and I know of a couple you will find in the R FAQ. On Sat, 16 Apr 2005, Liaw, Andy wrote: Because a data frame can hold different data types (even matrices) in different variables, one row of it can not be converted to a vector in general (where all elements need to be of the same type). Andy From: Fernando Saldanha Thanks, it's interesting reading. I also noticed that sw[, 1, drop = TRUE] is a vector (coerces to the lowest dimension) but sw[1, , drop = TRUE] is a one-row data frame (does not convert it into a list or vector) FS On 4/16/05, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote: You should look at ?[ and look very carefully at the drop argument. For your example sw[, 1] is the first component of the data frame, but sw[, 1, drop = FALSE] is a data frame consisting of just the first component, as mathematically fastidious people would expect. This is a convention, and like most arbitrary conventions it can be very useful most of the time, but some of the time it can be a very nasty trap. Caveat emptor. Bill Venables. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Fernando Saldanha Sent: Saturday, 16 April 2005 1:07 PM To: Submissions to R help Subject: [R] Getting subsets of a data frame I was reading in the Reference Manual about Extract.data.frame. There is a list of examples of expressions using [ and [[, with the outcomes. I was puzzled by the fact that, if sw is a data frame, then sw[, 1:3] is also a data frame, but sw[, 1] is just a vector. Since R has no scalars, it must be the case that 1 and 1:1 are the same: 1 == 1:1 [1] TRUE Then why isn't sw[,1] = sw[, 1:1] a data frame? FS __ R-help@stat.math.ethz.ch 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 -- 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
[R] [R-pkgs] Version 1.0-0 of the Rcmdr package
Dear list members, I've just uploaded version 1.0-0 of the Rcmdr package to CRAN. For people who haven't seen the package before, the R Commander provides a basic-statistics graphical user interface to R, based on the tcltk package. The new version incorporates a number of improvements to the R Commander interface, as documented in the CHANGES file distributed with the package, and I feel that the package is now sufficiently polished and stable to warrant finally bumping the version number to 1.0-0. More details are available at http://socserv.socsci.mcmaster.ca/jfox/Misc/Rcmdr/index.html and in the manual distributed with the package and accessible through the R Commander Help - Introduction to the R Commander menu. As usual, comments, suggestions, etc., are appreciated. John John Fox Department of Sociology McMaster University Hamilton, Ontario Canada L8S 4M4 905-525-9140x23604 http://socserv.mcmaster.ca/jfox ___ R-packages mailing list [EMAIL PROTECTED] https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Getting subsets of a data frame
On Sat, 16 Apr 2005, Prof Brian Ripley wrote: Perhaps Fernando will also note that is documented in ?[.data.frame, a slightly more appropriate reference than Bill's. It would be a good idea to read a good account of R's indexing: Bill Venables and I know of a couple you will find in the R FAQ. BTW, sw - swiss sw[1,,drop=TRUE] *is* a list (not as claimed, but as documented) sw[1, ] is a data frame sw[, 1] is a numeric vector. I should have pointed out that [.data.frame is in the See Also of Bill's reference. BTW to Andy: a list is a vector, and Kurt and I recently have been trying to correct documentation that means `atomic vector' when it says `vector'. (Long ago lists in R were pairlists and not vectors.) is.vector(list(a=1)) [1] TRUE On Sat, 16 Apr 2005, Liaw, Andy wrote: Because a data frame can hold different data types (even matrices) in different variables, one row of it can not be converted to a vector in general (where all elements need to be of the same type). Andy From: Fernando Saldanha Thanks, it's interesting reading. I also noticed that sw[, 1, drop = TRUE] is a vector (coerces to the lowest dimension) but sw[1, , drop = TRUE] is a one-row data frame (does not convert it into a list or vector) FS On 4/16/05, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote: You should look at ?[ and look very carefully at the drop argument. For your example sw[, 1] is the first component of the data frame, but sw[, 1, drop = FALSE] is a data frame consisting of just the first component, as mathematically fastidious people would expect. This is a convention, and like most arbitrary conventions it can be very useful most of the time, but some of the time it can be a very nasty trap. Caveat emptor. Bill Venables. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Fernando Saldanha Sent: Saturday, 16 April 2005 1:07 PM To: Submissions to R help Subject: [R] Getting subsets of a data frame I was reading in the Reference Manual about Extract.data.frame. There is a list of examples of expressions using [ and [[, with the outcomes. I was puzzled by the fact that, if sw is a data frame, then sw[, 1:3] is also a data frame, but sw[, 1] is just a vector. Since R has no scalars, it must be the case that 1 and 1:1 are the same: 1 == 1:1 [1] TRUE Then why isn't sw[,1] = sw[, 1:1] a data frame? FS __ R-help@stat.math.ethz.ch 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 -- 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 -- 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
[R] Re: Inverse of the Laplace Transform/Gaver Stehfest algorithm
Tolga Uzuner wrote:
Hi there,
Is there an implementation of the Gaveh Stehfest algorithm in R
somewhere ? Or some other inversion ?
Thanks,
Tolga
Well, at least here is Zakian's algorithm, for anyone who needs it:
Zakian-function(Fs,t){
# Fs is the function to be inverted and evaluated at t
a = c(12.83767675+1.666063445i,
12.22613209+5.012718792i,10.93430308+8.409673116i,
8.776434715+11.92185389i,5.225453361+15.72952905i)
K = c(-36902.08210+196990.4257i,
61277.02524-95408.62551i,-28916.56288+18169.18531i,
+4655.361138-1.901528642i,-118.7414011-141.3036911i)
ssum = 0.0
# Zakian's method does not work for t=0. Check that out.
if(t == 0){
print(ERROR: Inverse transform can not be calculated for t=0)
print(WARNING: Routine zakian() exiting.)
return(Error)}
# The for-loop below is the heart of Zakian's Inversion Algorithm.
for(j in 1:5){ssum = ssum + Re(K[j]*Fs(a[j]/t))}
return (2.0*ssum/t)
}
# InvLap(1/(s-1))=exp(t)
# check if Zakian(function(s){1/(s-1)},1)==exp(1)
lapfunc-function(s){1.0/(s-1.0)}
# Function Zakian(functobeinverted,t) is invoked.
Zakian(lapfunc,1.0)
__
R-help@stat.math.ethz.ch 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] Re: Inverse of the Laplace Transform/Gaver Stehfest algorithm
Tolga Uzuner wrote:
Tolga Uzuner wrote:
Hi there,
Is there an implementation of the Gaveh Stehfest algorithm in R
somewhere ? Or some other inversion ?
Thanks,
Tolga
Well, at least here is Zakian's algorithm, for anyone who needs it:
Zakian-function(Fs,t){
# Fs is the function to be inverted and evaluated at t
a = c(12.83767675+1.666063445i,
12.22613209+5.012718792i,10.93430308+8.409673116i,
8.776434715+11.92185389i,5.225453361+15.72952905i)
K = c(-36902.08210+196990.4257i,
61277.02524-95408.62551i,-28916.56288+18169.18531i,
+4655.361138-1.901528642i,-118.7414011-141.3036911i)
ssum = 0.0
# Zakian's method does not work for t=0. Check that out.
if(t == 0){
print(ERROR: Inverse transform can not be calculated for t=0)
print(WARNING: Routine zakian() exiting.)
return(Error)}
# The for-loop below is the heart of Zakian's Inversion Algorithm.
for(j in 1:5){ssum = ssum + Re(K[j]*Fs(a[j]/t))}
return (2.0*ssum/t)
}
# InvLap(1/(s-1))=exp(t)
# check if Zakian(function(s){1/(s-1)},1)==exp(1)
lapfunc-function(s){1.0/(s-1.0)}
# Function Zakian(functobeinverted,t) is invoked.
Zakian(lapfunc,1.0)
By the way, I am told This significance of this specification is that
Zakians Algorithm is very accurate for overdamped and slightly
underdamped systems. But it is not accurate for systems with prolonged
oscillations.
for those considering using it.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Element-wise multiplication
Dear members, The code I am writing heavily use element-wise multiplication of matrix and vectors, e.g. X, is nxm matrix e, is nx1 matrix Doing Z=X*e[,], I obtain a nxm matrix, Z, where each column of X is multiplied (element-wise) by e. Is this the best way to achieve the result I want? By best way, I mean the fastest way. By profiling my code, 45% of the time is spent by * and even a 30% speedup in obtaining Z would greatly benefit the total speed of the code. Aby suggestion is greatly appreciated. |Giuseppe Ragusa |University of California, San Diego Please avoid sending me Word or PowerPoint attachments. See http://www.gnu.org/philosophy/no-word-attachments.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] Element-wise multiplication
On 4/16/05, GiusVa [EMAIL PROTECTED] wrote: Dear members, The code I am writing heavily use element-wise multiplication of matrix and vectors, e.g. X, is nxm matrix e, is nx1 matrix Doing Z=X*e[,], I obtain a nxm matrix, Z, where each column of X is multiplied (element-wise) by e. Is this the best way to achieve the result I want? By best way, I mean the fastest way. By profiling my code, 45% of the time is spent by * and even a 30% speedup in obtaining Z would greatly benefit the total speed of the code. Aby suggestion is greatly appreciated. |Giuseppe Ragusa |University of California, San Diego Try this: mm - matrix(1, 1000, 1000) ee - matrix(1:100,nc=1) system.time(mm*ee[,], TRUE) [1] 0.26 0.02 0.28 NA NA system.time(mm*c(ee), TRUE) [1] 0.07 0.00 0.07 NA NA # if you have to do it multiple times with same ee: cee - c(ee) system.time(mm*cee, TRUE) [1] 0.06 0.00 0.06 NA NA __ R-help@stat.math.ethz.ch 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] Getting subsets of a data frame
I am reading as fast as I can! Just started with R five days ago. I found the following in the documentation: Although the default for 'drop' is 'TRUE', the default behaviour when only one _row_ is left is equivalent to specifying 'drop = FALSE'. To drop from a data frame to a list, 'drop = FALSE' has to (sic) specified explicitly. I think the exception mentioned in the first sentence is the reason for my confusion. I also think the second sentence is wrong and should have 'TRUE' instead of 'FALSE'. While it is true that a data frame is a list, it is not a list of numbers, but rather a list of columns, which, if I understand correctly, can be either vectors or matrices. So regardless of the value assigned to 'drop' the returned object is a list. When I asked why isn't sw[1, ] a list? I should have asked instead why isn't sw[1, ] a list of vectors? I did some experiments with a data frame a, where the columns are vectors (no matrix columns): is.data.frame(a) # just checking [1] TRUE a1- a[3, ] (is.data.frame(a1)) [1] TRUE (did not sop being a data frame) (is.list(a1)) [1] TRUE (but it is a list) a2- a[3, , drop=T] (is.data.frame(a2)) [1] FALSE (no longer a data frame) (is.list(a2)) [1] TRUE (but it is a list) a3- a[3, , drop=F] (is.data.frame(a3)) [1] TRUE(still a data frame) (is.list(a3)) [1] TRUE(but it is a list) I also tried: a2[1] $dates.num [1] 477032400 a3[1] dates.num 3 477032400 (notice the row name) attributes(a3[1]) $names [1] dates.num $class [1] data.frame $row.names [1] 3 attributes(a2[1]) $names [1] dates.num FS On 4/16/05, Prof Brian Ripley [EMAIL PROTECTED] wrote: On Sat, 16 Apr 2005, Prof Brian Ripley wrote: Perhaps Fernando will also note that is documented in ?[.data.frame, a slightly more appropriate reference than Bill's. It would be a good idea to read a good account of R's indexing: Bill Venables and I know of a couple you will find in the R FAQ. BTW, sw - swiss sw[1,,drop=TRUE] *is* a list (not as claimed, but as documented) sw[1, ] is a data frame sw[, 1] is a numeric vector. I should have pointed out that [.data.frame is in the See Also of Bill's reference. BTW to Andy: a list is a vector, and Kurt and I recently have been trying to correct documentation that means `atomic vector' when it says `vector'. (Long ago lists in R were pairlists and not vectors.) is.vector(list(a=1)) [1] TRUE On Sat, 16 Apr 2005, Liaw, Andy wrote: Because a data frame can hold different data types (even matrices) in different variables, one row of it can not be converted to a vector in general (where all elements need to be of the same type). Andy From: Fernando Saldanha Thanks, it's interesting reading. I also noticed that sw[, 1, drop = TRUE] is a vector (coerces to the lowest dimension) but sw[1, , drop = TRUE] is a one-row data frame (does not convert it into a list or vector) FS On 4/16/05, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote: You should look at ?[ and look very carefully at the drop argument. For your example sw[, 1] is the first component of the data frame, but sw[, 1, drop = FALSE] is a data frame consisting of just the first component, as mathematically fastidious people would expect. This is a convention, and like most arbitrary conventions it can be very useful most of the time, but some of the time it can be a very nasty trap. Caveat emptor. Bill Venables. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Fernando Saldanha Sent: Saturday, 16 April 2005 1:07 PM To: Submissions to R help Subject: [R] Getting subsets of a data frame I was reading in the Reference Manual about Extract.data.frame. There is a list of examples of expressions using [ and [[, with the outcomes. I was puzzled by the fact that, if sw is a data frame, then sw[, 1:3] is also a data frame, but sw[, 1] is just a vector. Since R has no scalars, it must be the case that 1 and 1:1 are the same: 1 == 1:1 [1] TRUE Then why isn't sw[,1] = sw[, 1:1] a data frame? FS __ R-help@stat.math.ethz.ch 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!
[R] help on extract variance components from the fitted model by lm
Hey, all: Do we have a convenient command(s) to extract the variance components from a fitted model by lm (actually it's a nexted model)? e.g.: using the following codes we could get MSA,MSB(A) and MSE. How to get the variance component estimates by command in R rather than calculations by hand? A-as.vector(rep(c(rep(1,5), rep(2,5), rep(3,5), rep(4,5), rep(5,5)),2)) B-as.vector(rep(c(rep(c(1,2,3,4,5),5)),2)) y-as.vector(c(15.5,15.2,14.2,14.3,15.8,6.2,7.2,6.6,6.2,5.6,15.4,13.9,13.4,12.5,13.2,10.9,12.5,12.3,11.0,12.3,7.5,6.7,7.2,7.6,6.3,14.9,15.2,14.2,14.3,16.4,7,8.4,7.8,7.6,7.4,14.4,13.3,14.8,14.1,15,11.3,12.7,11.7,12,13.3,6.7,7.3,6,7.6,7.1)) lm1-lm(y~factor(A)/factor(B)) anova(lm1) Thanks a lot! And have a good weekend! Regards, Wenqing __ R-help@stat.math.ethz.ch 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] Ack! Odd correlation matrix problem
Spencer Graves wrote: Does the following answer the question: cor(B, use=complete.obs) ** snip ** cor(B, use=pairwise.complete.obs) Yep. That's exactly the issue. I had thought the reference to casewise deletion in the help for complete.obs was referring solely to the two variables involved, not the entire dataset. The documentation might be a little clearer on this point for those just starting out in statistics, although I suppose it's only an issue if you're working with correlation matrices, which might imply you're really *not* just starting out in statistics, and should know better. -- Chris __ R-help@stat.math.ethz.ch 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] help on extract variance components from the fitted model by lm
On Sun, 2005-04-17 at 02:38 +0800, wenqing li wrote: Hey, all: Do we have a convenient command(s) to extract the variance components from a fitted model by lm (actually it's a nexted model)? e.g.: using the following codes we could get MSA,MSB(A) and MSE. How to get the variance component estimates by command in R rather than calculations by hand? A-as.vector(rep(c(rep(1,5), rep(2,5), rep(3,5), rep(4,5), rep (5,5)),2)) B-as.vector(rep(c(rep(c(1,2,3,4,5),5)),2)) y-as.vector(c (15.5,15.2,14.2,14.3,15.8,6.2,7.2,6.6,6.2,5.6,15.4,13.9,13.4,12.5,13.2, 10.9,12.5,12.3,11.0,12.3,7.5,6.7,7.2,7.6,6.3,14.9,15.2,14.2,14.3,16.4,7,8.4, 7.8,7.6,7.4,14.4,13.3,14.8,14.1,15,11.3,12.7,11.7,12,13.3,6.7,7.3,6,7.6,7.1)) lm1-lm(y~factor(A)/factor(B)) anova(lm1) Thanks a lot! And have a good weekend! Regards, Wenqing First, the use of as.vector() above is redundant: A - as.vector(rep(c(rep(1,5), rep(2,5), rep(3,5), rep(4,5), rep(5,5)),2)) str(A) num [1:50] 1 1 1 1 1 2 2 2 2 2 ... A1 - rep(c(rep(1,5), rep(2,5), rep(3,5), rep(4,5), rep(5,5)),2) str(A1) num [1:50] 1 1 1 1 1 2 2 2 2 2 ... all.equal(A, A1) [1] TRUE On your question, a couple of options: # First create a data frame df - data.frame(factor(A), factor(B), y) library(nlme) # Set your factors as the nested random effects lme1 - lme(y ~ 1, random = ~ 1 | A / B, data = df) summary(lme1) Linear mixed-effects model fit by REML Data: df AIC BIClogLik 147.4539 155.0211 -69.72693 Random effects: Formula: ~1 | A (Intercept) StdDev:3.797784 Formula: ~1 | B %in% A (Intercept) Residual StdDev: 0.3768259 0.6957023 Fixed effects: y ~ 1 Value Std.Error DF t-value p-value (Intercept)11 1.702937 25 6.45943 0 Standardized Within-Group Residuals: Min Q1Med Q3Max -1.7696293 -0.7042397 0.1521976 0.5708701 1.5679386 Number of Observations: 50 Number of Groups: A B %in% A 5 25 Also: library(ape) varcomp(lme1) A B Within 14.4231663 0.1419978 0.4840016 attr(,class) [1] varcomp Note in both cases, lme() is used, not lm(). These are referenced in Section 10.2 (pg 279) of MASS4 by Venables Ripley and in Section 4.2.3 (pg 167) of MEMSS by Pinheiro and Bates. 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] Ack! Odd correlation matrix problem
Spencer Graves wrote: Does the following answer the question: cor(B, use=complete.obs) ** snip ** cor(B, use=pairwise.complete.obs) Yep. That's exactly the issue. I had thought the reference to casewise deletion in the help for complete.obs was referring solely to the two variables involved, not the entire dataset. The documentation might be a little clearer on this point for those just starting out in statistics, although I suppose it's only an issue if you're working with correlation matrices, which might imply you're really *not* just starting out in statistics, and should know better. -- Chris __ R-help@stat.math.ethz.ch 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] Generating a binomial random variable correlated with a
Ted: Thank you for your help. All I want is a binomial random variable that is correlated with a normal random variable with specified correlation. By linear I mean the ordinary Pearson correlation. I tried the following two methods, in each case the resulting correlation is substantially less than the one specified. Method I: Generate two correlated normals (using Cholesky decomposition method) and dichotomize one (0/1) to get the binomial. Method II: Generate two correlated variables, one binomial and one normal using the Cholesky decomposition methods. Here is how I did: X - rnorm(100) Y - rnorm(100) r- 0.7 Y1 - X*r+Y*sqrt(1-r**2) cor(X,Y1) # Correlated normals using Cholesky decomposition cor(X0.84,Y1) # Method I ## X1 - rbinom(100,1,0.5) Y2 - X1*r+Y*sqrt(1-r**2) cor(X1,Y2); # Method II I would like to thank Ben from whom I received the following response: Are you computing the correlation between the continuous variable and the dichotomized variable with the formula for the biserial correlation? If not, that is probably the root of your problem. I looked at the biserial correlation which is a special case of Pearson correlation between a continuous and binomial random variable. I don't know how I can use it to generate the data. Any idea? Regards, Ashraf -Original Message- From: Ted Harding [mailto:[EMAIL PROTECTED] Sent: Saturday, April 16, 2005 3:22 AM To: Ashraf Chaudhary Cc: r-help@stat.math.ethz.ch Subject: RE: [R] Generating a binomial random variable correlated with a On 15-Apr-05 Ashraf Chaudhary wrote: Hi, I am posting this problem again (with some additional detail) as I am stuck and could not get it resolved as yet. I tried to look up in alternative sources but with no success. Here it is: I need to generate a binomial (binary 0/1) random variable linearly correlated with a normal random variable with a specified correlation. Off course, the correlation coefficient would not be same at each run because of randomness. If I generate two correlated normals with specified correlation and dichotomize one, the correlation of a normal and the binomial random variable would not be the same as specified. I greatly appreciate your help. Ashraf Hello Ashraf, I do not know what you mean by a binomial random variable linearly correlated with a normal random variable. You can certainly (and indeed your dichotomy method is one way) generate a binomial and a normal which are correlated. But apparently this gives a result which is not the same as specified: however, I cannot see in your description a specification which would violated by the result of doing so. You cannot expect a binomial variable to be such that, for instance, its expectation conditional on the value of a normal variable would be a linear function of the normal variable, since this would allow a situation where the expectation was greater than 1 or less than 0. But I wonder what else you could possibly mean by linearly correlated. Please therefore be more explicit about the specification of your problem! Trying to help, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 16-Apr-05 Time: 08:21:42 -- 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
Re: [R] How to create a vector with one, two, three, ...?
Sorry, I didn't get the question clear. What I meant is to create a character vector with length 200: one, two, three, ..., two hundred On 4/15/05, Federico Calboli [EMAIL PROTECTED] wrote: On Fri, 2005-04-15 at 14:30 -0400, Frank Duan wrote: Hi R people, I met a naive prolem. Could anyone give me a hint how to create such a vector with entries: one, two, three, ...? rvect - c(one, two, three) rvect [1] one two three Is it what you want? F -- Federico C. F. Calboli Department of Epidemiology and Public Health Imperial College, St Mary's Campus Norfolk Place, London W2 1PG Tel +44 (0)20 7594 1602 Fax (+44) 020 7594 3193 f.calboli [.a.t] imperial.ac.uk f.calboli [.a.t] gmail.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] How to create a vector with one, two, three, ...?
I may be wrong, but I am unaware of anyone that has created a number to text function in R. If you search Google: http://www.google.com/search?q=numbers+into+words There are various program examples, from VB to JavaScript to PHP, etc. It shouldn't be too hard to convert one of them to R. Most have fairly common constructs in terms of parsing and converting the numbers. Some of them handle decimals and currency formats as well. HTH, Marc Schwartz On Sat, 2005-04-16 at 18:01 -0400, Frank Duan wrote: Sorry, I didn't get the question clear. What I meant is to create a character vector with length 200: one, two, three, ..., two hundred On 4/15/05, Federico Calboli [EMAIL PROTECTED] wrote: On Fri, 2005-04-15 at 14:30 -0400, Frank Duan wrote: Hi R people, I met a naive prolem. Could anyone give me a hint how to create such a vector with entries: one, two, three, ...? rvect - c(one, two, three) rvect [1] one two three Is it what you want? F __ R-help@stat.math.ethz.ch 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] Unlogical phenomenon caused by locator
I was wondering if someone could help me!
Here is what I don't understand!
Why the text in the function isn't showing up before the locator action?
###
plot(1:10)
toto-function(){
cat(Why this text isn't showing up before the point selection )
point-locator(2)
return(point)
}
toto()
P.S.: I am using the most recent version of R for mac
Thanks a lot!
--
Sébastien Durand
Maîtrise en biologie
Université de Montréal
(514) 343-6864
Université du Québec à Montréal
(514) 987-3000 (1572#)
__
R-help@stat.math.ethz.ch 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] nls segmented model with unknown joint points
Hello, I am interested in fitting a segmented model with unknown joint points in nls and perhaps eventually in nlme. I can fit this model in sas (see below, joint points to be estimated are a41 and a41), but am unsure how to specify this in the nlm function. I would really appreciate any suggestions or example code. Thanks a lot. -andy proc nlin data=Stems.Trees; params b41=-3 b42=1.5 b43=-1.5 b44=50 a41=0.75 a42=0.1; term1 = (b41*(x - 1) + b42*(x**2 -1)); if (a41 - x) = 0 then term2 = (b43*(a41 - x)**2); else term2 = 0; if (a42 - x) =0 then term3 = (b44*(a42 - x)**2); else term3 = 0; model y = term1+term2+term3; run; __ R-help@stat.math.ethz.ch 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 create a vector with one, two, three, ...?
Dear Frank,
This was an interesting exercise. Here's a solution:
numbers2words - function(x){
helper - function(x){
digits - rev(strsplit(as.character(x), )[[1]])
nDigits - length(digits)
if (nDigits == 1) as.vector(ones[digits])
else if (nDigits == 2)
if (x = 19) as.vector(teens[digits[1]])
else trim(paste(tens[digits[2]],
Recall(as.numeric(digits[1]
else if (nDigits == 3) trim(paste(ones[digits[3]], hundred,
Recall(makeNumber(digits[2:1]
else {
nSuffix - ((nDigits + 2) %/% 3) - 1
if (nSuffix length(suffixes)) stop(paste(x, is too large!))
trim(paste(Recall(makeNumber(digits[
nDigits:(3*nSuffix + 1)])),
suffixes[nSuffix],
Recall(makeNumber(digits[(3*nSuffix):1]
}
}
trim - function(text){
gsub(^\ , , gsub(\ *$, , text))
}
makeNumber - function(...) as.numeric(paste(..., collapse=))
opts - options(scipen=100)
on.exit(options(opts))
ones - c(, one, two, three, four, five, six, seven,
eight, nine)
names(ones) - 0:9
teens - c(ten, eleven, twelve, thirteen, fourteen, fifteen,
sixteen, seventeen, eighteen, nineteen)
names(teens) - 0:9
tens - c(twenty, thirty, forty, fifty, sixty, seventy,
eighty,
ninety)
names(tens) - 2:9
x - round(x)
suffixes - c(thousand, million, billion, trillion)
if (length(x) 1) return(sapply(x, helper))
helper(x)
}
For example:
numbers2words(5673424350)
[1] fifty six trillion seven hundred thirty four billion two hundred
million four thousand three hundred fifty
numbers2words(c(567342, 604))
[1] five billion six hundred seventy three million four hundred twenty
thousand
[2] six hundred four
numbers2words(21:30)
[1] twenty one twenty two twenty three twenty four twenty
five
[6] twenty six twenty seven twenty eight twenty nine thirty
Note that if you want, you could go beyond trillions by adding to suffixes.
I hope that this does what you want,
John
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Frank Duan
Sent: Saturday, April 16, 2005 5:02 PM
To: [EMAIL PROTECTED]
Cc: r-help
Subject: Re: [R] How to create a vector with one, two,
three, ...?
Sorry, I didn't get the question clear. What I meant is to
create a character vector with length 200:
one, two, three, ..., two hundred
On 4/15/05, Federico Calboli [EMAIL PROTECTED] wrote:
On Fri, 2005-04-15 at 14:30 -0400, Frank Duan wrote:
Hi R people,
I met a naive prolem. Could anyone give me a hint how to
create such
a vector with entries: one, two, three, ...?
rvect - c(one, two, three)
rvect
[1] one two three
Is it what you want?
F
--
Federico C. F. Calboli
Department of Epidemiology and Public Health Imperial College, St
Mary's Campus Norfolk Place, London W2 1PG
Tel +44 (0)20 7594 1602 Fax (+44) 020 7594 3193
f.calboli [.a.t] imperial.ac.uk
f.calboli [.a.t] gmail.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-help@stat.math.ethz.ch 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] nls segmented model with unknown joint points
This is how I'd write the formula for use with nls/nlme: y ~ b41*(x - 1) + b42*(x^2 - 1) + ifelse((a41 - x) = 0, b43*(a41 - x)^2, 0) + ifelse((a42 - x) = 0, b44*(a42 - x)^2, 0) This is a direct translation from your funny foreign-looking code below that probably makes it clear what's going on. A more swish R form might be y ~ b41*(x - 1) + b42*(x^2 - 1) + b43*pmax(a41 - x, 0)^2 + b44*pmax(a42 - x, 0)^2 You mention nlm, too. Here you would use a function rather than a formula, but the idea is the same. V. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of andy Sent: Sunday, 17 April 2005 1:09 PM To: r-help@stat.math.ethz.ch Subject: [R] nls segmented model with unknown joint points Hello, I am interested in fitting a segmented model with unknown joint points in nls and perhaps eventually in nlme. I can fit this model in sas (see below, joint points to be estimated are a41 and a41), but am unsure how to specify this in the nlm function. I would really appreciate any suggestions or example code. Thanks a lot. -andy proc nlin data=Stems.Trees; params b41=-3 b42=1.5 b43=-1.5 b44=50 a41=0.75 a42=0.1; term1 = (b41*(x - 1) + b42*(x**2 -1)); if (a41 - x) = 0 then term2 = (b43*(a41 - x)**2); else term2 = 0; if (a42 - x) =0 then term3 = (b44*(a42 - x)**2); else term3 = 0; model y = term1+term2+term3; run; __ R-help@stat.math.ethz.ch 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] nls segmented model with unknown joint points
Thank you very much Dr. Venables, I'll give this a try. Regards- andy On Sun, 2005-04-17 at 13:36 +1000, [EMAIL PROTECTED] wrote: This is how I'd write the formula for use with nls/nlme: y ~ b41*(x - 1) + b42*(x^2 - 1) + ifelse((a41 - x) = 0, b43*(a41 - x)^2, 0) + ifelse((a42 - x) = 0, b44*(a42 - x)^2, 0) This is a direct translation from your funny foreign-looking code below that probably makes it clear what's going on. A more swish R form might be y ~ b41*(x - 1) + b42*(x^2 - 1) + b43*pmax(a41 - x, 0)^2 + b44*pmax(a42 - x, 0)^2 You mention nlm, too. Here you would use a function rather than a formula, but the idea is the same. V. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of andy Sent: Sunday, 17 April 2005 1:09 PM To: r-help@stat.math.ethz.ch Subject: [R] nls segmented model with unknown joint points Hello, I am interested in fitting a segmented model with unknown joint points in nls and perhaps eventually in nlme. I can fit this model in sas (see below, joint points to be estimated are a41 and a41), but am unsure how to specify this in the nlm function. I would really appreciate any suggestions or example code. Thanks a lot. -andy proc nlin data=Stems.Trees; params b41=-3 b42=1.5 b43=-1.5 b44=50 a41=0.75 a42=0.1; term1 = (b41*(x - 1) + b42*(x**2 -1)); if (a41 - x) = 0 then term2 = (b43*(a41 - x)**2); else term2 = 0; if (a42 - x) =0 then term3 = (b44*(a42 - x)**2); else term3 = 0; model y = term1+term2+term3; run; __ R-help@stat.math.ethz.ch 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
