Re: [R] Histograms of lots of variables
Please, experiment a little bit :-)
I have trouble using Win metafile graphics, so I just converted it to
jpeg. The code works - there was a copy/paste error. Try this with your
'total' dataframe - it works with my 'total' given there:
total - list(a=rnorm(1000),b=rnorm(1000),c=(10*runif(1000)-5))
lapply(1:3, function(.ind){
jpeg(filename=paste(file., .ind, .jpg, sep=''))
hist(total[[.ind]], main=paste('Histogram of',names(total)[.ind], sep=
), xlab=)
dev.off()
})
You can add the breaks argument inside the call of hist() as you would
normally do when plotting a histogram. However, consider it carefully -
it will most probably result into an error if your breaks won't fit all
of your data.
Petr
[EMAIL PROTECTED] napsal(a):
Thanks again Jim, I appreciate your time. I've been trying to debug the
code, but am running into a wall. I'm getting a syntax error after the line
containing the hist function. Here's the R session, any ideas?
Also, I'd like to be able to have each histogram use the same x-axis breaks
(0, 1, 2-3, 4-10, 11-20, 20) using the breaks argument, but it's not clear
to me where I'd put that in your code.
Thanks! Mark
total-read.csv(total_ec.csv) #creates a dataframe called total from a
.csv file
lapply(10:20, function(.ind){
+win.metafile(filename=paste(file., .ind, .wmf, sep=''))
+hist(total[[.ind]], main=colnames(total)[.ind], xlab=)
+ sep=''), xlab=)
Error: syntax error in:
hist(total[[.ind]], main=colnames(total)[.ind], xlab=)
sep='')
dev.off()
Error in dev.off() : cannot shut down device 1 (the null device)
})
On 4/16/07, jim holtman [EMAIL PROTECTED] wrote:
This will put the colnames on the histograms:
lapply(3:20, function(.ind){
win.metafile(filename=paste(file., .ind, .wmf, sep=''))
hist(dataframe[[.ind]], main=colnames(dataframe)[.ind], xlab=)
dev.off()
})
On 4/16/07, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote:
Hi Jim, thanks for your help. That looks like it might work, but a
couple of
things...
1) The resulting 73 (in my case) histograms will be named by their
variable
number, not by the variable name contained in the first row. Any way to
include the variable name in the resulting histogram?
2) How can I avoid having the resulting output.wmf files overwrite
each
other? I'd like to keep all of them, not just the last one.
Any ideas?
Thanks, Mark
On 4/16/07, jim holtman [EMAIL PROTECTED] wrote:
try:
win.metafile(file=output.wmf)
lapply(dataframe[3:20], hist)
dev.off()
On 4/16/07, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote:
Hi R-helpers,
I wish to produce frequency histograms of all of the variables in my
dataframe (except some identifying variables).
I have tried
hist(dataframe[,3:20])
to produce histograms of the 3rd through 20th variables in my
dataframe,
but
R doesn't like that.
Could anyone provide a suggestion?
Also, once I produce the histograms, I'd like to save them as
graphic
files
on my computer. How would I do that using code (rather than
Right-click
|
Save as metafile, which would be tedious for dozens of histograms).
Thanks, Mark
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Jim Holtman
Cincinnati, OH
+1 513 646 9390
What is the problem you are trying to solve?
--
Jim Holtman
Cincinnati, OH
+1 513 646 9390
What is the problem you are trying to solve?
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Petr Klasterecky
Dept. of Probability and Statistics
Charles University in Prague
Czech Republic
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Please help.... i know its trivial
Patrick Callier napsal(a): rownames is the command you want, I think. rownames(A) - seq(1:8) Yes. More generally rownames(A) - 1:length(A[,1]) And keep in mind that nobody is going to help you next time, if you won't give a clear and informative SUBJECT in your message. Here it should have been something like 'wrong row names in subset' or 'dataframe subset question' etc, certainly not 'please help me, it's trivial'... Btw, there is something like the posting guide for this list (see footer). Petr -- Petr Klasterecky Dept. of Probability and Statistics Charles University in Prague Czech Republic __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Setting rownames (was Please help.... i know its trivial)
On Tue, 17 Apr 2007, Petr Klasterecky wrote: Patrick Callier napsal(a): rownames is the command you want, I think. rownames(A) - seq(1:8) Yes. More generally rownames(A) - 1:length(A[,1]) Better rownames(A) - seq_len(nrow(A)) which works correctly if the length is 0 or there are zero columns, and is self-explanatory. Whenever you use 1:n you need to consider what you want if n = 0. And keep in mind that nobody is going to help you next time, if you won't give a clear and informative SUBJECT in your message. Here it should have been something like 'wrong row names in subset' or 'dataframe subset question' etc, certainly not 'please help me, it's trivial'... Btw, there is something like the posting guide for this list (see footer). And please change the subject line to something appropriate when replying. Petr -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Greek symbols in xtable rows
Thank You for answer. I Tried code that You provided in reply to my question
on both PC's with Linux and Widows OS.
On linux box output is:
xtable(diag(c($\\sigma_1^2$,$\\sigma_2^2$)))
% latex table generated in R 2.4.0 by xtable 1.4-3 package
% Tue Apr 17 09:18:31 2007
\begin{table}[ht]
\begin{center}
\begin{tabular}{rll}
\hline
1 2 \\
\hline
1 \$$\backslash$sigma\_1\verb|^|2\$ 0 \\
2 0 \$$\backslash$sigma\_2\verb|^|2\$ \\
\hline
\end{tabular}
\end{center}
\end{table}
System information:
sessionInfo ()
R version 2.4.0 (2006-10-03)
i686-redhat-linux-gnu
locale:
LC_CTYPE=lv_LV.UTF-8;LC_NUMERIC=C;LC_TIME=lv_LV.UTF-8;LC_COLLATE=lv_LV.UTF-8;LC_MONETARY=lv_LV.UTF-8;LC_MESSAGES=lv_LV.UTF-8;LC_PAPER=lv_LV.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=lv_LV.UTF-8;LC_IDENTIFICATION=C
attached base packages:
[1] tcltk methods stats graphics grDevices utils
[7] datasets base
other attached packages:
xtable pls ellipse
1.4-3 1.2-1 0.3-2
version
platform i686-redhat-linux-gnu
arch i686
os linux-gnu
system i686, linux-gnu
status
major 2
minor 4.0
year 2006
month 10
day03
svn rev39566
language R
version.string R version 2.4.0 (2006-10-03)
And there is an output form my windows PC:
mat - diag (c($\\sigma_1^2,$\\sigma_2^2$))
xtable (mat)
% latex table generated in R 2.4.1 by xtable 1.4-3 package
% Tue Apr 17 09:39:09 2007
\begin{table}[ht]
\begin{center}
\begin{tabular}{rll}
\hline
1 2 \\
\hline
1 \$$\backslash$sigma\_1\verb|^|2 0 \\
2 0 \$$\backslash$sigma\_2\verb|^|2\$ \\
\hline
\end{tabular}
\end{center}
\end{table}
sessionInfo ()
R version 2.4.1 (2006-12-18)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods
base
other attached packages:
xtable
1.4-3
version
_
platform i386-pc-mingw32
arch i386
os mingw32
system i386, mingw32
status
major 2
minor 4.1
year 2006
month 12
day18
svn rev40228
language R
version.string R version 2.4.1 (2006-12-18)
Thank You,
Andris Jankevics
On Pirmdiena, 16. Aprīlis 2007 22:38, Charles C. Berry wrote:
On Mon, 16 Apr 2007, Andris Jankevics wrote:
Dear R-helpers,
I am using xtable package to prepare a Latex code of some R tables.
Is this possible to have a greek symbols in xtable cells?
How can I get for example a string of : $\Delta$
$\Delta$
[1] $Delta$
And string: $\\Delta$
[1] $\\Delta$
Gives a latex aoutput like: \$$\backslash$Delta\$
The posting guide says
Sometimes it helps to provide a small example that someone can
actually run.
which you did not do.
This makes it hard for us to pinpoint the error.
For me, this works
mat - diag(c($\\sigma_1^2$,$\\sigma_2^2$))
xtable(mat)
% latex table generated in R 2.4.1 by xtable 1.4-2 package
% Mon Apr 16 12:27:54 2007
\begin{table}[ht]
\begin{center}
\begin{tabular}{rll}
\hline
1 2 \\
\hline
1 $\sigma_1^2$ 0 \\
2 0 $\sigma_2^2$ \\
\hline
\end{tabular}
\end{center}
\end{table}
and placing the text output in between
\documentclass{article}
\begin{document}
and
\end{document}
saving the file and running latex (or pdflatex) produces nicely set
output with the Greek letter 'sigma' properly rendered.
Thank You in advance
Andris Jankevics
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html and provide commented,
minimal, self-contained, reproducible code.
Charles C. Berry(858) 534-2098
Dept of Family/Preventive
Medicine E mailto:[EMAIL PROTECTED]UC San Diego
http://biostat.ucsd.edu/~cberry/ La Jolla, San Diego 92093-0901
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] It is possible to Expand RAM with HD space?
Jorge Cornejo-Donoso wrote: I'm using R in a winXp machine. I have problems with the size of my database (xyz), anyone know if is possible to use some kind of swap memory, in order to expand my real RAM and use HD space for the processes? This is done automatically by your operating system. You do not need to care about it. Uwe Ligges I'm trying to make a variogram, and I don't want to sacrifice the precision of the data, e.g. using mean or median by defined pixels or something like that (now I have to do that, or use other software like SAS). Thanks in advance! __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Extracting approximate Wald test (Chisq) from coxph(..frailty)
Dear List, How do I extract the approximate Wald test for the frailty (in the following example 17.89 value)? What about the P-values, other Chisq, DF, se(coef) and se2? How can they be extracted? ## kfitm1 Call: coxph(formula = Surv(time, status) ~ age + sex + disease + frailty(id, dist = gauss), data = kidney) coef se(coef) age0.00489 0.0150 sex -1.69703 0.4609 diseaseGN 0.17980 0.5447 diseaseAN 0.39283 0.5447 diseasePKD-1.13630 0.8250 frailty(id, dist = gauss se2Chisq DF age 0.0106 0.11 1.0 sex 0.3617 13.56 1.0 diseaseGN 0.3927 0.11 1.0 diseaseAN 0.3982 0.52 1.0 diseasePKD0.6173 1.90 1.0 frailty(id, dist = gauss17.89 12.1 p age 0.74000 sex 0.00023 diseaseGN 0.74000 diseaseAN 0.47000 diseasePKD0.17000 frailty(id, dist = gauss 0.12000 Iterations: 6 outer, 30 Newton-Raphson Variance of random effect= 0.493 Degrees of freedom for terms= 0.5 0.6 1.7 12.1 Likelihood ratio test=47.5 on 14.9 df, p=2.82e-05 n= 76 ## Thank you for your time. Thanks in advance. Mohammad Ehsanul Karim wildscop at yahoo dot com Institute of Statistical Research and Training University of Dhaka __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] colSum() in Matrix objects
Hi Jose,
Jose I'd like to simply add column-wise using Matrix objects (Csparse).
Jose It looks like one can apply mosty any base function to these objects
Jose (i.e., apply, colSums), but there is a nasty conversion to
traditional
Jose matrix objects if one does that.
not in this case, see below.
Jose Is there any workaround? I can see colSum listed in the help for
Class
colSums (final 's'!)
Jose 'CsparseMatrix' , but I wonder whether I'm using the default
colSums() or
Jose the one specific to CsparseMatrix...
#example
(z = Matrix(c(0,1,0,0), 10,10))
zr = rowSums(z)
class(zr) # numeric; I'd like it to be a CSparseMatrix object
selectMethod(colSums, class(z))
## or
showMethods(colSums)
both show you that you are using the class specific one.
However, why do you assume that colSums() should not return a
numeric vector? From the idea that colSums() and rowSums()
should be fast versions of apply(., marg, sum),
it must return a numeric vector, as it also does for
traditional matrices.
Are your objects so huge that even a 1-row {or 1-column} sparse
matrix would save a lot?
Regards,
Martin Maechler, ETH Zurich
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] background color
hi, I want to add different colors on the background of a classical plot. Each color is associated to an interval of the x axis. example: the background is red on the interval [1,10], blue on [11,20]. I try the rect function but it isn't appropriate for the background. Can any one can help me please? best regards. - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] It is possible to Expand RAM with HD space?
On Tue, 17 Apr 2007, Uwe Ligges wrote: Jorge Cornejo-Donoso wrote: I'm using R in a winXp machine. I have problems with the size of my database (xyz), anyone know if is possible to use some kind of swap memory, in order to expand my real RAM and use HD space for the processes? This is done automatically by your operating system. You do not need to care about it. Except perhaps that on his OS it is limited to 2GB, which can be increased to 3GB in some circumstances (see the rw-FAQ Q2.9). And anyone using less than 2GB RAM at current prices (except perhaps on a laptop) and is running out of RAM needs to buy some more: *far* cheaper than buying SAS. Uwe Ligges I'm trying to make a variogram, and I don't want to sacrifice the precision of the data, e.g. using mean or median by defined pixels or something like that (now I have to do that, or use other software like SAS). Thanks in advance! __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] help comparing two median with R
Dear R users, I am new to R and I would like to ask your help with the following topic. I have three sets of numeral data, 2 sets are paired and a third is independent of the other two. For each of these sets I have obtained their basic statistics (mean, median, stdv, range ...). Now I want to compare if these sets differ. I could compare the mean doing a basic T test . However, I was looking for a test to compare the medians using R. If that is possible I would love to hear the specifics. Thanks in advance for any help, Regards, pedro reche __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] background color
What do you mean by background? Maybe this is enough: plot(seq(-3,3,.01), dnorm(seq(-3,3,.01)), type=n, xlab=x, ylab=f(x), main=Normal density) polygon(x=c(-4,0,0,-4), y=c(-1,-1,.5,.5), col=red) polygon(x=c(4,0,0,4), y=c(-1,-1,.5,.5), col=blue) lines(seq(-3,3,.01), dnorm(seq(-3,3,.01)), type=l, lwd=2) Play a little bit with the polygon margins to get what you need. You can even generate them automatically based on your data. Petr yannig goude napsal(a): hi, I want to add different colors on the background of a classical plot. Each color is associated to an interval of the x axis. example: the background is red on the interval [1,10], blue on [11,20]. I try the rect function but it isn't appropriate for the background. Can any one can help me please? best regards. - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Petr Klasterecky Dept. of Probability and Statistics Charles University in Prague Czech Republic __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Use of argument '...'
Dear R list, I've read the function writing sections on both An introduction to R and R language Definition manuals but still don't understand why the following gives an error message: fun - function(x, ...) x + y fun(1, y=2) I get: Error in fun(1, y = 2) : object y not found I'd appreciate any help in understanding this. R version 2.4.1 (2006-12-18) i386-pc-mingw32 ... Juan Pablo Lewinger Department of Preventive Medicine Keck School of Medicine University of Southern California 1540 Alcazar Street, CHP-220 Los Angeles, CA 90089-9011, USA __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Use of argument '...'
Dear Juan,
You have to define your function like this:
fun - function (x,y=0, ...) {
x+y
}
The ... means that u can use the other parameters in par!
The variables, which you want to use, these you have to define in your
function!
In this function i put y=0 on a standard value zero,
if u try now fun(2), the result will be: 2
if u try fun(2,4), then the result will be: 6
Kind Regards,
Alin Soare
2007/4/17, Juan Lewinger [EMAIL PROTECTED]:
Dear R list,
I've read the function writing sections on both An introduction to R and
R language Definition manuals but still don't understand why the following
gives an error message:
fun - function(x, ...) x + y
fun(1, y=2)
I get:
Error in fun(1, y = 2) : object y not found
I'd appreciate any help in understanding this.
R version 2.4.1 (2006-12-18)
i386-pc-mingw32
...
Juan Pablo Lewinger
Department of Preventive Medicine
Keck School of Medicine
University of Southern California
1540 Alcazar Street, CHP-220
Los Angeles, CA 90089-9011, USA
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Use of argument '...'
Juan Lewinger wrote: Dear R list, I've read the function writing sections on both An introduction to R and R language Definition manuals but still don't understand why the following gives an error message: fun - function(x, ...) x + y fun(1, y=2) I get: Error in fun(1, y = 2) : object y not found I'd appreciate any help in understanding this. Putting it a bit sharply: The triple-dot mechanism is not a carte blanche to introduce arbitrary variables into a function body. If it did, then you could (intentionally or not) modify the behaviour of functions in weird and mysterious ways. As y is not among the formal argument names for fun, it will not be found inside ..., to find it in there, you need to be more explicit, as in x+list(...)$y. R version 2.4.1 (2006-12-18) i386-pc-mingw32 ... Juan Pablo Lewinger Department of Preventive Medicine Keck School of Medicine University of Southern California 1540 Alcazar Street, CHP-220 Los Angeles, CA 90089-9011, USA __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] help comparing two median with R
hi
+++
Now I want to compare if these sets differ. I could compare the mean
doing a basic T test . However, I was looking for a test to compare
the medians using R.
+++
@@@i hope but i am not sure that wilcoxon signed rank test may help :-)
cheers
Pedro A Reche [EMAIL PROTECTED]
Sent by: [EMAIL PROTECTED]
17-04-07 02:22 PM
To
r-help@stat.math.ethz.ch
cc
Subject
[R] help comparing two median with R
Dear R users,
I am new to R and I would like to ask your help with the following
topic. I have three sets of numeral data, 2 sets are paired and a
third is independent of the other two. For each of these sets I have
obtained their basic statistics (mean, median, stdv, range ...).
Now I want to compare if these sets differ. I could compare the mean
doing a basic T test . However, I was looking for a test to compare
the medians using R. If that is possible I would love to hear the
specifics.
Thanks in advance for any help,
Regards,
pedro reche
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
DISCLAIMER AND CONFIDENTIALITY CAUTION:\ \ This message and ...{{dropped}}
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] СНИЖЕНИЕ ВЕСА
ÇÎËÎÒÀß ÑÅÐÜÃÀ - ÈÇÁÀÂËÅÍÈÅ ÎÒ ËÈØÍÈÕ ÊÈËÎÃÐÀÌÌΠ(4 Êà  ÌÅÑ.)!! ÎÃÐÎÌÍÎÅ ÊÎËÈ×ÅÑÒÂÎ ËÞÄÅÉ ÂÎ ÂÑÅÌ ÌÈÐÅ ÏÎÕÓÄÅËÈ ÁËÀÃÎÄÀÐß ÇÎËÎÒÎÉ ÑÅÐÜÃÅ ÁÅÇ ÈÑÒßÇÀÞÙÈÕ ÄÈÅÒ, ÔÈÇÈ×ÅÑÊÈÕ ÓÏÐÀÆÍÅÍÈÉ, ÒÀÁËÅÒÎÊ! ÑÒÎÈÌÎÑÒÜ ÂÑÅÃÎ ËÈØÜ ÇÎÎÎ ÐÓÁËÅÉ!!! Ïðèíöèï ìåòîäà «Çîëîòàÿ Ñåðüãà» ñîçäàí íà îñíîâå ñåêðåòîâ òðàäèöèîííîé Êèòàéñêîé ìåäèöèíû. «Çîëîòàÿ ñåðüãà» óñòàíàâëèâàåòñÿ íà äëèòåëüíîå âðåìÿ è âîçäåéñòâóåò íà áèîëîãè÷åñêè àêòèâíûå òî÷êè óøíîé ðàêîâèíû.  ðåçóëüòàòå ýòîãî ðåçêî ñíèæàåòñÿ è ïîäàâëÿåòñÿ àïïåòèò. Ìîçã ïîëó÷àåò ñèãíàë, ÷òî ÷åëîâåê ñûò. ÏÐÅÈÌÓÙÅÑÒÂÀ ÄÀÍÍÎÃÎ ÌÅÒÎÄÀ: Âàì íå ïðèäåòñÿ çàíèìàòüñÿ ñïîðòîì è ïðèíèìàòü òàáëåòêè. Ïîõóäåòü ëåãêî Âû ïîñåùàåòå ñïåöèàëèñòà âñåãî 1 ðàç. Ïîäàâëÿåòñÿ ÷óâñòâî ãîëîäà è àïïåòèò. Âû èçáàâèòåñü îò ïðîáëåìíûõ çîí. Ñàìûé äîñòóïíûé è íåäîðîãîé ìåòîä ïîõóäåíèÿ. Ñíèæåíèå âåñà îò 3 äî 10 êã çà îäèí ìåñÿö. Ïîñëå îêîí÷àíèÿ êóðñà âåñ íå íàáèðàåòñÿ. Ýôôåêòèâíîñòü ìåòîäà çàâèñèò îò îñîáåííîñòåé îðãàíèçìà, íàñëåäñòâåííîé ñêëîííîñòè, ñòðîåíèÿ æèðîâîé òêàíè, è â êàæäîì ñëó÷àå îïðåäåëÿåòñÿ èíäèâèäóàëüíî. Çâîíèòå íàì, ÷òîáû óçíàòü, êàêèõ ðåçóëüòàòîâ ìîæíî äîáèòüñÿ â êîíêðåòíî Âàøåì ñëó÷àå!! Íàø òåëåôîí (495) 739 34~39 êðóãëîñóòî÷íî ñò. ì. Òóøèíñêàÿ, óë. Ìèòèíñêàÿ, ä.43, ì-í Ìèòèíî ì. «Ëþáëèíî» (100 ì îò ìåòðî), óë. Íîâîðîññèéñêàÿ, ä.28 [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] help comparing two median with R
Pedro A Reche wrote: Dear R users, I am new to R and I would like to ask your help with the following topic. I have three sets of numeral data, 2 sets are paired and a third is independent of the other two. For each of these sets I have obtained their basic statistics (mean, median, stdv, range ...). Now I want to compare if these sets differ. I could compare the mean doing a basic T test . However, I was looking for a test to compare the medians using R. If that is possible I would love to hear the specifics. Hi Pedro, You can use the Mann-Whitney test (wilcox with two samples), but you would have to check that the second and third moments of the variable distributions were the same, I think. Jim __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] About PLR
Hello Sir/Mam I am Nitish Kumar Mishra from IMTECH, Chandigarh, India. I want to calculate the Principal Component Analysis(PCA), PLS of the input file and find top 30 PCAs of this file using PLS in R. I am asking regarding Linux(Red Hat 9) operating system. I have downloaded PLS from Crains site and try to installed it, but don't installed it. How I can download and installed PLS(for PCA and PLS) in R using Linux. Please me how I can installed PLR in R(Linux). Thanking you. -- Nitish Kumar Mishra Junior Research Fellow BIC, IMTECH, Chandigarh, India E-Mail Address: [EMAIL PROTECTED] [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Find zeros of analytic functions
Robin Hankin wrote: If iterative methods are appropriate, it's perhaps worth pointing out that Newton-Rapheson works nicely for complex functions. Hmmm... I think there are many cases where Newton-Raphson diverges for complex functions, like those that generate beautiful fractals. Alberto Monteiro __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] help comparing two median with R
Thomas Lumley wrote:
On Tue, 17 Apr 2007, Robert McFadden wrote:
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Jim Lemon
Sent: Tuesday, April 17, 2007 12:37 PM
To: Pedro A Reche
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] help comparing two median with R
Pedro A Reche wrote:
Dear R users,
I am new to R and I would like to ask your help with the following
topic. I have three sets of numeral data, 2 sets are paired and a
third is independent of the other two. For each of these sets I have
obtained their basic statistics (mean, median, stdv, range ...).
Now I want to compare if these sets differ. I could compare
the mean
doing a basic T test . However, I was looking for a test to compare
the medians using R. If that is possible I would love to
hear the
specifics.
Hi Pedro,
You can use the Mann-Whitney test (wilcox with two
samples), but you would have to check that the second and
third moments of the variable distributions were the same, I think.
Jim
Use Mann-Whitney U test, but remember about 2 assumption:
1. samples come from continuous distribution (there are no tied
obserwations)
2. distributions are identical in shape. It's very similar to t-test but
Mann-Whitney U test is not as affected by violation of the homogeneity of
variance assumption as t-test is.
This turns out not to be quite correct.
If the two distributions differ only by a location shift then the
hypothesis that the shift is zero is equivalent to the medians being the
same (or the means, or the 3.14159th percentile), and the Mann-Whitney U
test will test this hypothesis. Otherwise the Mann-Whitney U test does not
test for equal medians.
The assumption that the distributions are continuous is for convenience --
it makes the distribution of the test statistic easier to calculate and
otherwise R uses a approximation. The assumption of a location shift is
critical -- otherwise it is easy to construct three data sets x,y,z so
that the Mann-Whitney U test thinks x is larger than y, y is larger than z
and z is larger than x (Google for Efron Dice). That is, the Mann-Whitney
U test cannot be a test for any location statistic.
There actually is an exact test for the median that does not assume a
location shift: dichotomize your data at the pooled median to get a 2x2
table of above/below median by group, and do Fisher's exact test on the
table. This is almost never useful (because it doesn't come with an
interval estimate), but is interesting because it (and the generalizations
to other quantiles) is the only exactly distribution-free location test
that does not have the 'non-transitivity' problem of the Mann-Whitney U
test. I believe this median test is attributed to Mood, but I have not
seen the primary source.
-thomas
The Mood test is so inefficient that its use is no longer recommended:
@Article{fri00sho,
author = {Freidlin, Boris and Gastwirth, Joseph L.},
title ={Should the median test be retired from
general use?},
journal = American Statistician,
year = 2000,
volume = 54,
number = 3,
pages ={161-164},
annote = {inefficiency of Mood median test}
}
The points that Thomas and Brian have made are certainly correct, if one
is truly interested in testing for differences in medians or means. But
the Wilcoxon test provides a valid test of x y more generally. The
test is consonant with the Hodges-Lehmann estimator: the median of all
possible differences between an X and a Y.
Frank
--
Frank E Harrell Jr Professor and Chair School of Medicine
Department of Biostatistics Vanderbilt University
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] extracting intercept from ppr fit
Hi, Is there a way, documented or not, to extract the intercept term (the alpha_0 the MASS book) from a ppr() (Projection Persuit Regression) fit? Thanks, Vadim ## Example: n - 1000 data - data.frame(x=rnorm(n), y=rnorm(n)) a - 10 data$z - evalq(a + atan(x + y) + rnorm(n), data) data.ppr - ppr(z ~ x + y, data=data, nterms=1) ## how to extract a = 10 from data.ppr? [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Modelling Heteroscedastic Multilevel Models
I think there are many who can help, but this question is quite vague.
This assumes we have access to the book you note and can make sense of
your question w/o sample data.
If you cannot find a sample data set please create a sample data file.
However, there are so many sample data sets in the mlmRev package and in
other places I doubt you will need to do this. For example, see the
egsingle or star data files that are education specific. But, if you for
some reason cannot do either at least give a good substantive
description of your data and the problem you are trying to solve.
In the code you have below, you have a random intercept for each school,
but you remove the intercept in the fixed portion of the call. Also,
does it make sense to model Sex as random? This is a repeatable factor
(I hope), how can it be treated as a random draw from a population?
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Rense
Nieuwenhuis
Sent: Monday, April 16, 2007 4:37 PM
To: r-help@stat.math.ethz.ch
Subject: [R] Modelling Heteroscedastic Multilevel Models
Dear ListeRs,
I am trying to fit a heteroscedastic multilevel model using
lmer{lme4- package). Take, for instance, the (fictive) model below.
lmer(test.result ~ homework + Sex -1 + (1 | School))
Suppose that I suspect the error terms in the predicted
values to differ between men and women (so, on the first
level). In order to model this, I want the 'Sex'-variable to
be random on the first level, as described in Snijders
Bosker, page 110.
Does anybody know if this is possible and how this can be
done using R?
Many thanks in advance.
Rense Nieuwenhuis
PS. Please excuse me for not providing a self-contained
example. I couldn't find a data-set in the lme4-package that
fitted my question.
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] extracting intercept from ppr fit
Hi, Is there a way, documented or not, to extract the intercept term (the alpha_0 the MASS book) from a ppr() (Projection Persuit Regression) fit? Thanks, Vadim [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] About PLR
From within R, you can give the command install.packages(pls) and R will download and install it for you (as long as you have access to the Internet). To install an already downloaded package, you can use R CMD INSTALL pls_2.0-0.tar.gz in a terminal window. -- Bjørn-Helge Mevik __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] help comparing two median with R
On Tue, 17 Apr 2007, Robert McFadden wrote: -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Jim Lemon Sent: Tuesday, April 17, 2007 12:37 PM To: Pedro A Reche Cc: r-help@stat.math.ethz.ch Subject: Re: [R] help comparing two median with R Pedro A Reche wrote: Dear R users, I am new to R and I would like to ask your help with the following topic. I have three sets of numeral data, 2 sets are paired and a third is independent of the other two. For each of these sets I have obtained their basic statistics (mean, median, stdv, range ...). Now I want to compare if these sets differ. I could compare the mean doing a basic T test . However, I was looking for a test to compare the medians using R. If that is possible I would love to hear the specifics. Hi Pedro, You can use the Mann-Whitney test (wilcox with two samples), but you would have to check that the second and third moments of the variable distributions were the same, I think. Jim Use Mann-Whitney U test, but remember about 2 assumption: 1. samples come from continuous distribution (there are no tied obserwations) 2. distributions are identical in shape. It's very similar to t-test but Mann-Whitney U test is not as affected by violation of the homogeneity of variance assumption as t-test is. Sorry, the other way round, as in R the (Welch)unequal-variance t test is the default, so it makes no 'homogeneity of variance assumption'. The Wilcoxon/Mann-Whitney test is sensitive to differences in shape, in rather complex ways, more so than t.test. (I've posted simulated comparisons here in the past.) I don't think anyone is answering the question here. If the distributions are identical in shape, the Wilcoxon/Mann-Whitney test is measuring differences in means, medians ... (they are all the same), and if they are not, it is not measuring differences in medians. To test specifically for a difference in medians, you can do almost the same thing as the Welch two-sample t-test, but use estimates of the asympototic standard errors of the medians (for which you will need density estimates at the medians). But I can only see this as being worth doing where the median is of substantive interest (as it could just be for lifetime distributions), and statistical significance is of importance (rather than useful improvement). -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] background color
On 4/17/07, yannig goude [EMAIL PROTECTED] wrote: I want to add different colors on the background of a classical plot. Each color is associated to an interval of the x axis. example: the background is red on the interval [1,10], blue on [11,20]. I try the rect function but it isn't appropriate for the background.. You can use rect together with par(usr) like this: usr - par(usr) plot(1:20, type = n) rect(1, usr[3], 10, usr[4], col = red) rect(11, usr[3], 20, usr[4], col = blue) points(1:20) There is also an example of this using polygon in place of rect here: http://www.mayin.org/ajayshah/KB/R/html/g5.html and an example of doing it with lattice graphics using xyplot.zoo (the same idea would work with xyplot) mid way through the examples in: library(zoo) ?xyplot.zoo __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] sample size 2-sample proportion test
Hi R-users, I want to calculate the sample size needed to carry out a 2-sample proprotion test (with alfa=0.05, beta=0.8) 1.- probability of success in subpopulation A: 0.8 2.- probability of success in subpopulation B: 0.05 3.- percentage of population in subpopulation A = 5%, 4.- percentage of population in subpopulation B= 95% (the rest, A and B complementary) (Note that it is dificult to find people from subpopulation A) Using library(Hmisc) (or others if exist?), I would like to know which one of these three options (if any) should I use: # FIRST: bsamsize(p1=0.8, p2=0.05, fraction=0.05, alpha=.05, power=.8) #n1n2 # 1.46 27.82 # SECOND: samplesize.bin(alpha=0.05, beta=0.8, pit=0.8, pic=0.05, rho=0.05) # n= 42 # THIRD: To calculate sample size assuming fraction=0.5 and then multiplying by 0.05: bsamsize(p1=0.8, p2=0.05, fraction=0.5, alpha=.05, power=.8) # n1 n2 # 5.47 5.47 # so that n=6*100/5=120 Taking into account the difficulty of obtaining the samples, which method would be recomended? Or in other words, how many people should I find, around 30-40 or 120? Thanks a lot in advance, Berta __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Special Limited Time Discount Offer from Chapman Hall/CRC Press
EXTRA SPECIAL 50% OFF DISCOUNT ON SELECT TITLES including: Robust Statistical Methods with R by Jana Jureckova and Jan Pickek List Price: $84.95 SALE PRICE: $42.48 An Introduction to the Bootstrap by Bradley Efron and R.J. Tibshirani List Price: $99.95 SALE PRICE: $49.98 HURRY to take advantage of this limited time offer. Click on this link to begin saving on these and other select titles: http://www.crcpress.com/util/showdiscountpromobyname.asp?promoname=Web+S pecials Nadja Nadja English Senior Marketing Manager CRC Press - Taylor Francis Group an informa business 6000 Broken Sound Parkway NW, Suite 300 Boca Raton, FL 33487 Phone (561) 998-2525 FAX(561) 998-0876 [EMAIL PROTECTED] www.crcpress.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 and provide commented, minimal, self-contained, reproducible code.
[R] Inverse of one function
Dear Partecipants to the list, I am searching a R function, which can calculate the inverse of one real value function. Does it exists a R code in order to make it? Many thanks for any kind of help and for Your availability. Enrico Foscolo -- Passa a Infostrada. ADSL e Telefono senza limiti e senza canone Telecom http://click.libero.it/infostrada __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] My First Function: cryptic error message
Weiwei Shi writes:
I cannot reproduce the error w/o the data called probands
I realized this when I posted; sharing the data with anyone is not
possible (I hope you understand). I'm hoping a more experienced useR
can spot a syntax or procedural error that I made.
I have made some progress. It appears that moving the function lower
down in the script stops the error from coming up and I'm now on to
writing other functions with other errors ;)
Any ideas why the placement in the script might make a difference?
Thanks,
Joel
--
Joel J. Adamson
Biostatistician
Pediatric Psychopharmacology Research Unit
Massachusetts General Hospital
Boston, MA 02114
(617) 643-1432
(303) 880-3109
The information transmitted in this electronic communication...{{dropped}}
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] A problem about all possible sequences
Dear All Suppose a sequence of length 10 generated by the following rule: the first element is 0 with 50% of probability or 1 with the same probability; the second element likewise; and so on. Is there some R command to obtain all possible different sequences formed by the above rule? I am aware that one could write a small program to do that, but I am speculating about whether a command is already existent. Thanks in advance, Paul __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] A problem about all possible sequences
On 4/17/07, Robin Hankin [EMAIL PROTECTED] wrote:
f - function(n){expand.grid(rep(list(0:1),n))}
f(10)
[snip]
Thanks, Robin, that is it!
Paul
On 17 Apr 2007, at 15:26, Paul Smith wrote:
Dear All
Suppose a sequence of length 10 generated by the following rule: the
first element is 0 with 50% of probability or 1 with the same
probability; the second element likewise; and so on.
Is there some R command to obtain all possible different sequences
formed by the above rule? I am aware that one could write a small
program to do that, but I am speculating about whether a command is
already existent.
Thanks in advance,
Paul
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-
guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] my.symbols and zoomplot functions now available in TeachingDemos package
In February (2/26) Jonathan Lees asked a question about plotting symbols whose shape/aspect ratio is independent of of the user coordinate system. A couple of other posts have also been made asking about plotting symbols of a given fixed size/aspect ratio. The latest version of the TeachingDemos package (1.5) now has a function called my.symbols that is similar to the symbols function except that it allows you to specify your own symbols. Your symbol can be a matrix or list with x and y coordinates (on the range -1, 1) that specify the shape of the symbol (to be connected with lines), or it can be a function (additional arguments are passed to this function) that plots your symbol. The my.symbols function will then place your symbols at the designated x,y coordinates at the specified size with the aspect ratio independent of the user coordinates. Some example symbol functions are also provided and the help page shows several examples of potential uses. Hopefully others will find this useful. There have also been several questions recently about changing the range of a plot after the plot has been created. TeachingDemos now also contains a function zoomplot that will do this (either zooming into a section of the existing plot, or zooming out to show parts that were previously clipped). This function is a quick and dirty hack and should not replace properly defining the ranges for serious graphs, but may be helpful for quick exploring. On a similar note, there is also a new function (plot2script) for creating a script file to recreate the current plot that can then be editted and rerun to make minor changes. These functions have only been tested on windows (R versions 2.4.1 and 2.5.0beta), so I would appreciate if someone could try them on other platforms and let me know if they work or not. Hope these are useful to others out there, -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare [EMAIL PROTECTED] (801) 408-8111 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] A problem about all possible sequences
Paul
f - function(n){expand.grid(rep(list(0:1),n))}
f(10)
[snip]
HTH
Robin
On 17 Apr 2007, at 15:26, Paul Smith wrote:
Dear All
Suppose a sequence of length 10 generated by the following rule: the
first element is 0 with 50% of probability or 1 with the same
probability; the second element likewise; and so on.
Is there some R command to obtain all possible different sequences
formed by the above rule? I am aware that one could write a small
program to do that, but I am speculating about whether a command is
already existent.
Thanks in advance,
Paul
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-
guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] PROC DISCRIM vs. lda( ) in MASS
On Tue, 17 Apr 2007, Greg Tarpinian wrote: Hello, I am using WinXP, R version 2.3.1, and SAS for PC version 8.1. I have mostly used SAS over the last 4 years and would like to compare the output of PROC DISCRIM to that of lda( ) with respect to a very specific aspect. My data have k=3 populations and there are 3 variates in the feature space. When using using the code PROC DISCRIM DATA = FOO OUT = FOO_OUT OUTSTAT = FOOSTAT METHOD = NORMAL LIST POOL = YES PCOV MANOVA; CLASS STRATA; PRIORS EQUAL; VAR X1 X2 X3; RUN; I am able to easily obtain the linear discriminant functions for the strata which allow computation of the three discriminant scores for a given observation. This information is contained in WORK.FOOTSTAT and may be extracted by subsetting: DATA LDFUNC; SET FOOSTAT(WHERE = (_TYPE_ = LINEAR)); RUN; To actually implement the linear discriminant functions takes a bit more formatting, but there it is. My question: Where is this information stored in R? I completely Nowhere (R does not store things like SAS does). However, if you read the help page you would see that the 'scaling' component of the result comes close. There are lots of definitional questions about what is meant by 'the individual lin. discrim. functions': lda() is support software for a book, and that explains its definitions. understand that predict( ) or predict.lda( ) are the preferable ways to obtain a classification prediction for new observations. Or even the linear discriminants applied to new or old observations. I still want to see the individual lin. discrim. functions and work with them myself. I have been using x.lda - lda(Strata ~ X1 + X2+ X3, data = foo.frame) to construct the analysis. Much thanks, Greg __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] help comparing two median with R
On Tue, 17 Apr 2007, Robert McFadden wrote: -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Jim Lemon Sent: Tuesday, April 17, 2007 12:37 PM To: Pedro A Reche Cc: r-help@stat.math.ethz.ch Subject: Re: [R] help comparing two median with R Pedro A Reche wrote: Dear R users, I am new to R and I would like to ask your help with the following topic. I have three sets of numeral data, 2 sets are paired and a third is independent of the other two. For each of these sets I have obtained their basic statistics (mean, median, stdv, range ...). Now I want to compare if these sets differ. I could compare the mean doing a basic T test . However, I was looking for a test to compare the medians using R. If that is possible I would love to hear the specifics. Hi Pedro, You can use the Mann-Whitney test (wilcox with two samples), but you would have to check that the second and third moments of the variable distributions were the same, I think. Jim Use Mann-Whitney U test, but remember about 2 assumption: 1. samples come from continuous distribution (there are no tied obserwations) 2. distributions are identical in shape. It's very similar to t-test but Mann-Whitney U test is not as affected by violation of the homogeneity of variance assumption as t-test is. This turns out not to be quite correct. If the two distributions differ only by a location shift then the hypothesis that the shift is zero is equivalent to the medians being the same (or the means, or the 3.14159th percentile), and the Mann-Whitney U test will test this hypothesis. Otherwise the Mann-Whitney U test does not test for equal medians. The assumption that the distributions are continuous is for convenience -- it makes the distribution of the test statistic easier to calculate and otherwise R uses a approximation. The assumption of a location shift is critical -- otherwise it is easy to construct three data sets x,y,z so that the Mann-Whitney U test thinks x is larger than y, y is larger than z and z is larger than x (Google for Efron Dice). That is, the Mann-Whitney U test cannot be a test for any location statistic. There actually is an exact test for the median that does not assume a location shift: dichotomize your data at the pooled median to get a 2x2 table of above/below median by group, and do Fisher's exact test on the table. This is almost never useful (because it doesn't come with an interval estimate), but is interesting because it (and the generalizations to other quantiles) is the only exactly distribution-free location test that does not have the 'non-transitivity' problem of the Mann-Whitney U test. I believe this median test is attributed to Mood, but I have not seen the primary source. -thomas __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Export multiple data files from R
You didn't describe the exact format of the .gpr files.
There are 32 heading lines in each file, which are now hidden from R
community.
You 'skip' 31 of them in read.table, and one more plays header ('header=T').
Since you are using read.table, your files are usual ascii files.
You shoud use 'cat' and/or 'write.table' R functions to accomplish your
task. There are some more functions exist.
Using them, you should explicitly write those 32 lines, and then write the
rest contents of the resulting data frame.
Jenny persson wrote:
Dear R-users,
I have 10 data files in gpr format (dat1.gpr,...,dat10.gpr). I want to
read in these files one by one in R and then add one extra column (called
log) to each data file as below
data.file=sort(dir(path ='C:/Documents and Settings/ Mina
dokument/data1, pattern = .gpr$,full.names = TRUE))
num.data.files- length(data.file)
num.data.files
i=1
### read one data file
data-read.table(file =
data.file[i],skip=31,header=T,sep='\t',na.strings=NA)
### Define the log ratio using values in column 2 resp 8
log=as.matrix(log((data[,2])/(data[,8])))
### append column called log to data frame data, for the reading data
file
data=cbind(data,log)
### Read remaining data files
for (i in 2:num.data.files) {
data-read.table(file=data.file[i],header=T,skip=31,sep='\t',na.strings=NA)
log=as.matrix(log((data[,2])/(data[,8])))
data=cbind(data,log)
}
Now I want to export these files (with an extra column in each) as
gpr-files in a folder called data2 but don't know exactly how to do it,
can you help me out ?
Thanks for your help,
Jenny
--
View this message in context:
http://www.nabble.com/Export-multiple-data-files-from-R-tf3587093.html#a10036489
Sent from the R help mailing list archive at Nabble.com.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] read.spss (package foreign) and SPSS 15.0 files
--- Charilaos Skiadas [EMAIL PROTECTED] wrote: On Apr 16, 2007, at 10:41 AM, John Kane wrote: --- Charilaos Skiadas [EMAIL PROTECTED] wrote: It is not an export option, it is a save as option. I don't have a 14 to check, but on a 15 I just go to File - Save As, and change the Save as type field to Comma Delimited (.*.csv). (I suppose tab delimited would be another option). Then there are two check- boxes below the window that allow a bit further customizing, one of them is about using value labels where defined instead of data values. I'm now back on a machine with SPSS 14. No csv option that I can see. Perhaps an enhancement to v15. I don't have a 14, but I did check a 13 today and you are correct, no csv option is there, which in my opinion is quite unacceptable for a statistical package that is on its 13/14th version. But there was an option for Excel 97 and ..., and that seemed to allow using value labels instead of the values (again you have to check the corresponding box). So perhaps that would be an option. Not my problem at the moment but a colleague pointed out that he had 750+ variables. Excel handles 256. Actually my problem is now a SAS one where I can get a clean csv export but lose the variable labels. My crude workaround was just to do a proc contents and cut-and-paste the results. A pain but it worked. I've got to figure out why I cannot get Hmisc to work! Thanks for the suggestion. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] A problem about all possible sequences
On 17-Apr-07 14:26:15, Paul Smith wrote:
Dear All
Suppose a sequence of length 10 generated by the following rule:
the first element is 0 with 50% of probability or 1 with the
same probability; the second element likewise; and so on.
You don't say whether the elements of the sequence are independent,
but plausibility suggests that this may be what you intend.
In which case:
Is there some R command to obtain all possible different sequences
formed by the above rule?
while(TRUE){print(sample(c(0,1),10,replace=TRUE))}
and just wait!
(Expected time to wait: about 7700 iterations, I think).
I am aware that one could write a small
program to do that, but I am speculating about whether a command
is already existent.
Taking my tongue out of my cheek, however, it's not clear what you
are really asking for.
If really you want to generate those sequences randomly according
to your probabilistic description, retaining as you go only those
which have not been sampled before, until you have all 2^10
possible sequences, then my suggestion above is not the way
to do it! And as far as I know there is not an R function which
does this by proceeding in exactly that way.
Better to recognise that your random scheme means that each
possible sequence is equally likely with all the others, and
so you can do the equivalent by sampling 1024 from (1:1024)
without replacement, i.e. putting (1:1024) in random order.
Then the binary representation of each element is such a
sequence.
So
S-sample((1:1024),1024)
is an existing R function which does the heart of the job.
(It remains to convert each integer K in S to binary form,
but as far as I know there is not an R function to convert
an integer K directly into a vector of binary 0/1 with a
given number of digits, i.e. not the equivalent of
to.binary(13,10) -- c(0,0,0,0,0,0,1,1,0,1)
except maybe in some special package, so I think you'll end
up writing your own for this bit anyway).
It gets more interesting if your example is just an illustraton,
and what you really want is more general.
E.g. if the different 0/1 outcomes in the 10 positions do
not have the same probabilities, but are still independent,
then you have to do more spadework (and again I'm pretty
sure there is no simple function in R to do it).
In that case it's definitely a programming job.
Even more so if the successive 0/1 outcomes are not independent,
whether P[0] = P[1] = 0.5 in each position or not. So again
a prgramming job.
Since you seem to be quite willling to do the programming
if necessary, I won't try to spoil your fun on that front!
Best wishes,
Ted.
E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 094 0861
Date: 17-Apr-07 Time: 16:40: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
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Greek symbols in xtable rows
Roger Bivand wrote:
Andris Jankevics andza at osi.lv writes:
Thank You for answer. I Tried code that You provided in reply to my question
on both PC's with Linux and Widows OS.
On linux box output is:
...
System information:
sessionInfo ()
R version 2.4.0 (2006-10-03)
i686-redhat-linux-gnu
locale:
LC_CTYPE=lv_LV.UTF-8;LC_NUMERIC=C;LC_TIME=lv_LV.UTF-8;LC_COLLATE=lv_LV.UTF-
8;LC_MONETARY=lv_LV.UTF-8;LC_MESSAGES=lv_LV.UTF-8;LC_PAPER=lv_LV.UTF-
8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=lv_LV.UTF-
8;LC_IDENTIFICATION=C
attached base packages:
[1] tcltk methods stats graphics grDevices utils
[7] datasets base
other attached packages:
xtable pls ellipse
1.4-3 1.2-1 0.3-2
And there is an output form my windows PC:
mat - diag (c($\\sigma_1^2,$\\sigma_2^2$))
xtable (mat)
... (replicates for me on same system)
sessionInfo ()
R version 2.4.1 (2006-12-18)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods
base
other attached packages:
xtable
1.4-3
Yes, I have also been bitten by the upgrade to xtable, which between 1.4-2
and 1.4-3 added code to try to be cleverer, but has broken output of LaTeX
markup directly. The offending code starts around line 177 in print.xtable.R
(#based on contribution from Jonathan Swinton jonathan at swintons.net in e-
mail dated Wednesday, January 17, 2007). I did try to write a
sanitize.text.function= solution, but failed, and backed off to an earlier
version. Could the maintainer David Dahl, please address this, and include a
relevant test? (I would CC him, but am travelling and posting via gmane to
keep
thr thread together - I apologize for pruning, but gmane won't post
otherwise).
So a solution is to install an earlier version of xtable from the package
archives, a harder but feasible task for Windows.
The solution is to use a workable sanitize.text.function.
mat - diag(c($\\sigma_1^2$,$\\sigma_2^2$))
tbl - xtable(mat)
print.xtable(tbl, sanitize.text.function = function(x) x)
% latex table generated in R 2.5.0 by xtable 1.4-3 package
% Tue Apr 17 11:00:05 2007
\begin{table}[ht]
\begin{center}
\begin{tabular}{rll}
\hline
1 2 \\
\hline
1 $\sigma_1^2$ 0 \\
2 0 $\sigma_2^2$ \\
\hline
\end{tabular}
\end{center}
\end{table}
Best,
Jim
Roger
...
Thank You,
Andris Jankevics
On Pirmdiena, 16. Aprīlis 2007 22:38, Charles C. Berry wrote:
On Mon, 16 Apr 2007, Andris Jankevics wrote:
mat - diag(c($\\sigma_1^2$,$\\sigma_2^2$))
xtable(mat)
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
James W. MacDonald, M.S.
Biostatistician
Affymetrix and cDNA Microarray Core
University of Michigan Cancer Center
1500 E. Medical Center Drive
7410 CCGC
Ann Arbor MI 48109
734-647-5623
**
Electronic Mail is not secure, may not be read every day, and should not be
used for urgent or sensitive issues.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Greek symbols in xtable rows
On Tue, 17 Apr 2007, James W. MacDonald wrote:
Roger Bivand wrote:
Andris Jankevics andza at osi.lv writes:
Thank You for answer. I Tried code that You provided in reply to my question
on both PC's with Linux and Widows OS.
On linux box output is:
...
System information:
sessionInfo ()
R version 2.4.0 (2006-10-03)
i686-redhat-linux-gnu
locale:
LC_CTYPE=lv_LV.UTF-8;LC_NUMERIC=C;LC_TIME=lv_LV.UTF-8;LC_COLLATE=lv_LV.UTF-
8;LC_MONETARY=lv_LV.UTF-8;LC_MESSAGES=lv_LV.UTF-8;LC_PAPER=lv_LV.UTF-
8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=lv_LV.UTF-
8;LC_IDENTIFICATION=C
attached base packages:
[1] tcltk methods stats graphics grDevices utils
[7] datasets base
other attached packages:
xtable pls ellipse
1.4-3 1.2-1 0.3-2
And there is an output form my windows PC:
mat - diag (c($\\sigma_1^2,$\\sigma_2^2$))
xtable (mat)
... (replicates for me on same system)
sessionInfo ()
R version 2.4.1 (2006-12-18)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods
base
other attached packages:
xtable
1.4-3
Yes, I have also been bitten by the upgrade to xtable, which between 1.4-2
and 1.4-3 added code to try to be cleverer, but has broken output of LaTeX
markup directly. The offending code starts around line 177 in print.xtable.R
(#based on contribution from Jonathan Swinton jonathan at swintons.net in e-
mail dated Wednesday, January 17, 2007). I did try to write a
sanitize.text.function= solution, but failed, and backed off to an earlier
version. Could the maintainer David Dahl, please address this, and include a
relevant test? (I would CC him, but am travelling and posting via gmane to
keep
thr thread together - I apologize for pruning, but gmane won't post
otherwise).
So a solution is to install an earlier version of xtable from the package
archives, a harder but feasible task for Windows.
The solution is to use a workable sanitize.text.function.
mat - diag(c($\\sigma_1^2$,$\\sigma_2^2$))
tbl - xtable(mat)
print.xtable(tbl, sanitize.text.function = function(x) x)
% latex table generated in R 2.5.0 by xtable 1.4-3 package
% Tue Apr 17 11:00:05 2007
\begin{table}[ht]
\begin{center}
\begin{tabular}{rll}
\hline
1 2 \\
\hline
1 $\sigma_1^2$ 0 \\
2 0 $\sigma_2^2$ \\
\hline
\end{tabular}
\end{center}
\end{table}
OK, thanks, I was trying too hard to be clever in my attempts. Now back to
many Rnw files to add the sanatize... functions ...
Roger
Best,
Jim
Roger
...
Thank You,
Andris Jankevics
On Pirmdiena, 16. Aprīlis 2007 22:38, Charles C. Berry wrote:
On Mon, 16 Apr 2007, Andris Jankevics wrote:
mat - diag(c($\\sigma_1^2$,$\\sigma_2^2$))
xtable(mat)
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
James W. MacDonald, M.S.
Biostatistician
Affymetrix and cDNA Microarray Core
University of Michigan Cancer Center
1500 E. Medical Center Drive
7410 CCGC
Ann Arbor MI 48109
734-647-5623
**
Electronic Mail is not secure, may not be read every day, and should not be
used for urgent or sensitive issues.
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: [EMAIL PROTECTED]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] read.spss (package foreign) and SPSS 15.0 files
John Kane wrote: --- Charilaos Skiadas [EMAIL PROTECTED] wrote: On Apr 16, 2007, at 10:41 AM, John Kane wrote: --- Charilaos Skiadas [EMAIL PROTECTED] wrote: It is not an export option, it is a save as option. I don't have a 14 to check, but on a 15 I just go to File - Save As, and change the Save as type field to Comma Delimited (.*.csv). (I suppose tab delimited would be another option). Then there are two check- boxes below the window that allow a bit further customizing, one of them is about using value labels where defined instead of data values. I'm now back on a machine with SPSS 14. No csv option that I can see. Perhaps an enhancement to v15. I don't have a 14, but I did check a 13 today and you are correct, no csv option is there, which in my opinion is quite unacceptable for a statistical package that is on its 13/14th version. But there was an option for Excel 97 and ..., and that seemed to allow using value labels instead of the values (again you have to check the corresponding box). So perhaps that would be an option. Not my problem at the moment but a colleague pointed out that he had 750+ variables. Excel handles 256. Actually my problem is now a SAS one where I can get a clean csv export but lose the variable labels. My crude workaround was just to do a proc contents and cut-and-paste the results. A pain but it worked. I've got to figure out why I cannot get Hmisc to work! And note that Hmisc has other SAS import options besides sas.get that keep labels. Some of them require you to run PROC CONTENTS CNTLOUT= to save metadata in a dataset. Frank Thanks for the suggestion. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Frank E Harrell Jr Professor and Chair School of Medicine Department of Biostatistics Vanderbilt University __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] GREP - Choosing values between two borders
Hello, I import datas from an file with: readLines But I need only a part of all measurments of this file. These are between two borders START and END. Can you tell me the syntax of grep(), to choose values between two borders? My R Code was not succesful, and I can't finde anything in the help. Thank's a lot. Felix # R-CODE ### file- file-content Measure - grep([START-END],file) #Measure - grep([START|END],file) FILE-CONTENT ## EXAM NUM:2 - EXAM #1 ASTIG:-2.4D AXIS:4.8 START OF HEIGHT DATA 0 0.0 0. 0 0.1 0.00055643 9 4.9 1.67278117 9 5.0 1.74873257 10 0.0 0. 10 0.1 0.00075557 99 5.3 1.94719490 END OF HEIGHT DATA X POS:-0.299mm Y POS:0.442mm Z POS:-0.290mm - EXAM #2 ASTIG:-2.4D AXIS:4.8 START OF HEIGHT DATA 0 0.0 0. 0 0.1 0.00055643 9 4.9 1.67278117 9 5.0 1.74873257 10 0.0 0. 10 0.1 0.00075557 99 5.3 1.94719490 END OF HEIGHT DATA X POS:-0.299mm Y POS:0.442mm Z POS:-0.290mm __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] installing R on Ubuntu
I am not sure what the question is but maybe this provides a starting point? http://cran.r-project.org/bin/linux/ubuntu/README HTH, Roland On 4/16/07, Erin Hodgess [EMAIL PROTECTED] wrote: __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Extracting approximate Wald test (Chisq) from coxph(..frailty)
Assign the output of coxph to some object, and use the $ extractor function to obtain what you need. ie: rtfm - coxph(formula = Surv(time, status) ~ age + sex + disease + frailty(id, dist = gauss), data = kidney) Age - coef(rtfm)[age] OR Sex - rtfm$coef[sex] Hope this helps. Paul Mohammad Ehsanul Karim wrote: Dear List, How do I extract the approximate Wald test for the frailty (in the following example 17.89 value)? What about the P-values, other Chisq, DF, se(coef) and se2? How can they be extracted? ## kfitm1 Call: coxph(formula = Surv(time, status) ~ age + sex + disease + frailty(id, dist = gauss), data = kidney) coef se(coef) age0.00489 0.0150 sex -1.69703 0.4609 diseaseGN 0.17980 0.5447 diseaseAN 0.39283 0.5447 diseasePKD-1.13630 0.8250 frailty(id, dist = gauss se2Chisq DF age 0.0106 0.11 1.0 sex 0.3617 13.56 1.0 diseaseGN 0.3927 0.11 1.0 diseaseAN 0.3982 0.52 1.0 diseasePKD0.6173 1.90 1.0 frailty(id, dist = gauss17.89 12.1 p age 0.74000 sex 0.00023 diseaseGN 0.74000 diseaseAN 0.47000 diseasePKD0.17000 frailty(id, dist = gauss 0.12000 Iterations: 6 outer, 30 Newton-Raphson Variance of random effect= 0.493 Degrees of freedom for terms= 0.5 0.6 1.7 12.1 Likelihood ratio test=47.5 on 14.9 df, p=2.82e-05 n= 76 ## Thank you for your time. Thanks in advance. Mohammad Ehsanul Karim wildscop at yahoo dot com Institute of Statistical Research and Training University of Dhaka __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- View this message in context: http://www.nabble.com/Extracting-approximate-Wald-test-%28Chisq%29-from-coxph%28..frailty%29-tf3589257.html#a10038426 Sent from the R help mailing list archive at Nabble.com. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Extracting approximate Wald test (Chisq) from coxph(..frailty)
On Tue, 17 Apr 2007, Mohammad Ehsanul Karim wrote: Dear List, How do I extract the approximate Wald test for the frailty (in the following example 17.89 value)? The example you give silently invokes print.coxph() to produce that output. You _can_ use tmp - capture.output( print( your example ) ) and then further process tmp. A _better_ solution for most purposes is to look at the object produced by coxph() and figure out how to calculate the Wald statistic from that object. See ?coxph.object and ?str Another tactic is to look at how print.coxph() does its work and use the code in it to produce just the output you desire. Look at page( survival:::print.coxph, print ) What about the P-values, other Chisq, DF, se(coef) and se2? How can they be extracted? ## kfitm1 Call: coxph(formula = Surv(time, status) ~ age + sex + disease + frailty(id, dist = gauss), data = kidney) coef se(coef) age0.00489 0.0150 sex -1.69703 0.4609 diseaseGN 0.17980 0.5447 diseaseAN 0.39283 0.5447 diseasePKD-1.13630 0.8250 frailty(id, dist = gauss se2Chisq DF age 0.0106 0.11 1.0 sex 0.3617 13.56 1.0 diseaseGN 0.3927 0.11 1.0 diseaseAN 0.3982 0.52 1.0 diseasePKD0.6173 1.90 1.0 frailty(id, dist = gauss17.89 12.1 p age 0.74000 sex 0.00023 diseaseGN 0.74000 diseaseAN 0.47000 diseasePKD0.17000 frailty(id, dist = gauss 0.12000 Iterations: 6 outer, 30 Newton-Raphson Variance of random effect= 0.493 Degrees of freedom for terms= 0.5 0.6 1.7 12.1 Likelihood ratio test=47.5 on 14.9 df, p=2.82e-05 n= 76 ## Thank you for your time. Thanks in advance. Mohammad Ehsanul Karim wildscop at yahoo dot com Institute of Statistical Research and Training University of Dhaka __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Charles C. Berry(858) 534-2098 Dept of Family/Preventive Medicine E mailto:[EMAIL PROTECTED] UC San Diego http://biostat.ucsd.edu/~cberry/ La Jolla, San Diego 92093-0901 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] A problem about all possible sequences
Paul Smith wrote: Dear All Suppose a sequence of length 10 generated by the following rule: the first element is 0 with 50% of probability or 1 with the same probability; the second element likewise; and so on. Is there some R command to obtain all possible different sequences formed by the above rule? I am aware that one could write a small program to do that, but I am speculating about whether a command is already existent. levels(interaction(c(0,1), c(0,1), c(0,1), c(0,1), c(0,1), c(0,1), c(0,1), c(0,1), c(0,1), c(0,1), sep=)) Thanks in advance, Paul __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Chuck Cleland, Ph.D. NDRI, Inc. 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] PROC DISCRIM vs. lda( ) in MASS
Hello, I am using WinXP, R version 2.3.1, and SAS for PC version 8.1. I have mostly used SAS over the last 4 years and would like to compare the output of PROC DISCRIM to that of lda( ) with respect to a very specific aspect. My data have k=3 populations and there are 3 variates in the feature space. When using using the code PROC DISCRIM DATA = FOO OUT = FOO_OUT OUTSTAT = FOOSTAT METHOD = NORMAL LIST POOL = YES PCOV MANOVA; CLASS STRATA; PRIORS EQUAL; VAR X1 X2 X3; RUN; I am able to easily obtain the linear discriminant functions for the strata which allow computation of the three discriminant scores for a given observation. This information is contained in WORK.FOOTSTAT and may be extracted by subsetting: DATA LDFUNC; SET FOOSTAT(WHERE = (_TYPE_ = LINEAR)); RUN; To actually implement the linear discriminant functions takes a bit more formatting, but there it is. My question: Where is this information stored in R? I completely understand that predict( ) or predict.lda( ) are the preferable ways to obtain a classification prediction for new observations. I still want to see the individual lin. discrim. functions and work with them myself. I have been using x.lda - lda(Strata ~ X1 + X2+ X3, data = foo.frame) to construct the analysis. Much thanks, Greg __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Greek symbols in xtable rows
Andris Jankevics andza at osi.lv writes: Thank You for answer. I Tried code that You provided in reply to my question on both PC's with Linux and Widows OS. On linux box output is: ... System information: sessionInfo () R version 2.4.0 (2006-10-03) i686-redhat-linux-gnu locale: LC_CTYPE=lv_LV.UTF-8;LC_NUMERIC=C;LC_TIME=lv_LV.UTF-8;LC_COLLATE=lv_LV.UTF- 8;LC_MONETARY=lv_LV.UTF-8;LC_MESSAGES=lv_LV.UTF-8;LC_PAPER=lv_LV.UTF- 8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=lv_LV.UTF- 8;LC_IDENTIFICATION=C attached base packages: [1] tcltk methods stats graphics grDevices utils [7] datasets base other attached packages: xtable pls ellipse 1.4-3 1.2-1 0.3-2 And there is an output form my windows PC: mat - diag (c($\\sigma_1^2,$\\sigma_2^2$)) xtable (mat) ... (replicates for me on same system) sessionInfo () R version 2.4.1 (2006-12-18) i386-pc-mingw32 locale: LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: xtable 1.4-3 Yes, I have also been bitten by the upgrade to xtable, which between 1.4-2 and 1.4-3 added code to try to be cleverer, but has broken output of LaTeX markup directly. The offending code starts around line 177 in print.xtable.R (#based on contribution from Jonathan Swinton jonathan at swintons.net in e- mail dated Wednesday, January 17, 2007). I did try to write a sanitize.text.function= solution, but failed, and backed off to an earlier version. Could the maintainer David Dahl, please address this, and include a relevant test? (I would CC him, but am travelling and posting via gmane to keep thr thread together - I apologize for pruning, but gmane won't post otherwise). So a solution is to install an earlier version of xtable from the package archives, a harder but feasible task for Windows. Roger ... Thank You, Andris Jankevics On Pirmdiena, 16. Aprīlis 2007 22:38, Charles C. Berry wrote: On Mon, 16 Apr 2007, Andris Jankevics wrote: mat - diag(c($\\sigma_1^2$,$\\sigma_2^2$)) xtable(mat) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] help comparing two median with R
-Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Jim Lemon Sent: Tuesday, April 17, 2007 12:37 PM To: Pedro A Reche Cc: r-help@stat.math.ethz.ch Subject: Re: [R] help comparing two median with R Pedro A Reche wrote: Dear R users, I am new to R and I would like to ask your help with the following topic. I have three sets of numeral data, 2 sets are paired and a third is independent of the other two. For each of these sets I have obtained their basic statistics (mean, median, stdv, range ...). Now I want to compare if these sets differ. I could compare the mean doing a basic T test . However, I was looking for a test to compare the medians using R. If that is possible I would love to hear the specifics. Hi Pedro, You can use the Mann-Whitney test (wilcox with two samples), but you would have to check that the second and third moments of the variable distributions were the same, I think. Jim Use Mann-Whitney U test, but remember about 2 assumption: 1. samples come from continuous distribution (there are no tied obserwations) 2. distributions are identical in shape. It's very similar to t-test but Mann-Whitney U test is not as affected by violation of the homogeneity of variance assumption as t-test is. Rob __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Tcltk
I have problems with Tcl/Tk in R 2.4.1, when running it on Windows Vista
(see error message below).
Regards, Sofia
library(tcltk)
Loading Tcl/Tk interface ... Error in fun(...) : Can't find a usable
init.tcl in the following directories:
{C:\Program\R\R-2.4.1/Tcl/lib/tcl8.4}
{C:\Program\R\R-2.4.1/Tcl/lib/tcl8.4} C:/Program/R/R-2.4.1/Tcl/lib/tcl8.4
C:/Program/R/R-2.4.1/Tcl/lib/tcl8.4
This probably means that Tcl wasn't installed properly.
Error: .onLoad failed in 'loadNamespace' for 'tcltk'
Error: package/namespace load failed for 'tcltk'
_
Sofia Wikström, PhD
AquaBiota Water Research
Svante Arrhenius väg 21A, SE-104 05 Stockholm, Sweden
Phone: (+46) 8 16 10 07
[EMAIL PROTECTED]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Modelling Heteroscedastic Multilevel Models
May I suggest the Exam data from the mlmRev package as an example.
If you wish to have a random effect for Sex by school you could write
the model as
lmer(test.result ~ homework + Sex + (Sex|school))
which gives correlated random effects for the overall achievement in
schools and the differential effect of Sex by school. Alternatively
you could write
lmer(test.result ~ homework + Sex + (1|school) + (1|Sex:school))
which gives uncorrelated effects for the intercept by school and for
the effect of sex by school.
On 4/17/07, Doran, Harold [EMAIL PROTECTED] wrote:
I think there are many who can help, but this question is quite vague.
This assumes we have access to the book you note and can make sense of
your question w/o sample data.
If you cannot find a sample data set please create a sample data file.
However, there are so many sample data sets in the mlmRev package and in
other places I doubt you will need to do this. For example, see the
egsingle or star data files that are education specific. But, if you for
some reason cannot do either at least give a good substantive
description of your data and the problem you are trying to solve.
In the code you have below, you have a random intercept for each school,
but you remove the intercept in the fixed portion of the call. Also,
does it make sense to model Sex as random? This is a repeatable factor
(I hope), how can it be treated as a random draw from a population?
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Rense
Nieuwenhuis
Sent: Monday, April 16, 2007 4:37 PM
To: r-help@stat.math.ethz.ch
Subject: [R] Modelling Heteroscedastic Multilevel Models
Dear ListeRs,
I am trying to fit a heteroscedastic multilevel model using
lmer{lme4- package). Take, for instance, the (fictive) model below.
lmer(test.result ~ homework + Sex -1 + (1 | School))
Suppose that I suspect the error terms in the predicted
values to differ between men and women (so, on the first
level). In order to model this, I want the 'Sex'-variable to
be random on the first level, as described in Snijders
Bosker, page 110.
Does anybody know if this is possible and how this can be
done using R?
Many thanks in advance.
Rense Nieuwenhuis
PS. Please excuse me for not providing a self-contained
example. I couldn't find a data-set in the lme4-package that
fitted my question.
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] GREP - Choosing values between two borders
You can adapt this to your situation: http://finzi.psych.upenn.edu/R/Rhelp02a/archive/22195.html On 4/17/07, Felix Wave [EMAIL PROTECTED] wrote: Hello, I import datas from an file with: readLines But I need only a part of all measurments of this file. These are between two borders START and END. Can you tell me the syntax of grep(), to choose values between two borders? My R Code was not succesful, and I can't finde anything in the help. Thank's a lot. Felix # R-CODE ### file- file-content Measure - grep([START-END],file) #Measure - grep([START|END],file) FILE-CONTENT ## EXAM NUM:2 - EXAM #1 ASTIG:-2.4D AXIS:4.8 START OF HEIGHT DATA 0 0.0 0. 0 0.1 0.00055643 9 4.9 1.67278117 9 5.0 1.74873257 10 0.0 0. 10 0.1 0.00075557 99 5.3 1.94719490 END OF HEIGHT DATA X POS:-0.299mm Y POS:0.442mm Z POS:-0.290mm - EXAM #2 ASTIG:-2.4D AXIS:4.8 START OF HEIGHT DATA 0 0.0 0. 0 0.1 0.00055643 9 4.9 1.67278117 9 5.0 1.74873257 10 0.0 0. 10 0.1 0.00075557 99 5.3 1.94719490 END OF HEIGHT DATA X POS:-0.299mm Y POS:0.442mm Z POS:-0.290mm __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] langage R
2007/4/16, elyakhlifi mustapha [EMAIL PROTECTED]: bonjour, je vous écris pour savoir si vous aviez trouver réponse pour l'écriture du test de Newman Keuls sous R car j'en ai aussi besoin merci encore cordialement. - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Hello, Maybe this post can help you, http://tolstoy.newcastle.edu.au/R/help/04/09/4292.html Rod. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Greatest common divisor of two numbers
Dear Sir/Madam: Could you please let me know which function shall I use to get the Greatest common divisor of two numbers. Thank you so much for your attention to this matter, and i look forward to hear from you soon. Regards; Abou == AbouEl-Makarim Aboueissa, Ph.D. Assistant Professor of Statistics Department of Mathematics Statistics University of Southern Maine 96 Falmouth Street P.O. Box 9300 Portland, ME 04104-9300 Tel: (207) 228-8389 Email: [EMAIL PROTECTED] [EMAIL PROTECTED] Office: 301C Payson Smith __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Greatest common divisor of two numbers
S Poetry has a function for that. Patrick Burns [EMAIL PROTECTED] +44 (0)20 8525 0696 http://www.burns-stat.com (home of S Poetry and A Guide for the Unwilling S User) AbouEl-Makarim Aboueissa wrote: Dear Sir/Madam: Could you please let me know which function shall I use to get the Greatest common divisor of two numbers. Thank you so much for your attention to this matter, and i look forward to hear from you soon. Regards; Abou == AbouEl-Makarim Aboueissa, Ph.D. Assistant Professor of Statistics Department of Mathematics Statistics University of Southern Maine 96 Falmouth Street P.O. Box 9300 Portland, ME 04104-9300 Tel: (207) 228-8389 Email: [EMAIL PROTECTED] [EMAIL PROTECTED] Office: 301C Payson Smith __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] how to estimate dose from respond given drc package result
Dear all, I can use the very nice drc package (multdrc()) to model and plot a dataframe containing dose and response values. I can also use predict.drc() to yield response values given a dose. I need to do the opposite, estimate a dose given the response. The general predict documentation seems to say that this is possible, but it does not appear that predict.drc has that capability. This makes sense because it is easy to compute y values from x values, given the fitted equation, but it is hard to do the opposite. It seems that I need to do some kind of successive approximation. Is there a package that takes a multdrc generated model and estimates dose given response? Thank you all for your time, John __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] about contrasts
Hi all R-users Please, I would like to learn about 'contrasts'. I do not know much about importance and several types (of contrasts) over regression and data analysis. I would like references for my study! Thanks in advanced... klebyn ___ Experimente já e veja as novidades. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Greatest common divisor of two numbers
Hi, wouldn't it be enough for simple purposes (as stated in the original message: Greatest common divisor of two numbers) to use: gcd - function(a,b) ifelse (b==0, a, gcd(b, a %% b)) gcd(12,4) [1] 4 gcd(4,12) [1] 4 gcd(123456789,987654321) [1] 9 gcd(987654321,123456789) [1] 9 Best, Roland On 4/17/07, Patrick Burns [EMAIL PROTECTED] wrote: S Poetry has a function for that. Patrick Burns [EMAIL PROTECTED] +44 (0)20 8525 0696 http://www.burns-stat.com (home of S Poetry and A Guide for the Unwilling S User) AbouEl-Makarim Aboueissa wrote: Dear Sir/Madam: Could you please let me know which function shall I use to get the Greatest common divisor of two numbers. Thank you so much for your attention to this matter, and i look forward to hear from you soon. Regards; Abou == AbouEl-Makarim Aboueissa, Ph.D. Assistant Professor of Statistics Department of Mathematics Statistics University of Southern Maine 96 Falmouth Street P.O. Box 9300 Portland, ME 04104-9300 Tel: (207) 228-8389 Email: [EMAIL PROTECTED] [EMAIL PROTECTED] Office: 301C Payson Smith __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] value of complexity parameter in ridge regression
Hi, What is the optimum range to look for a value of lambda while doing ridge regression. Can/ should lambda be greater than 1 ? I have conflicting (or what appears conflicting to me) sources that use lambda = 0, without any upper limit, but that makes the search space infinite.. right ?? So, perhaps my question is: is there an upper limit to lambda. Does the value of lambda convey something about my data ? Thanks a lot, Sikander - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] ICA for less common data
Hi everyone, I am not sure this is the appropriate list I should put this question to, but I hope you will re-direct me to the most appropriate one if necessary. I am doing an independent component analysis on a dataset that represents different metrics for patchreefs such as depth, area, volume, relative relief, shape index and rugosity. Doing different multivariate analyses and correlations among variables we discovered that the average reef depth data is heavily bi-modal, and rugosity behaves differently in each of the depth populations. Other reef metrics behaves differently with depth as well. Since not all metrics are really independent we did a PCA analysis followed by cluster analysis and we tried to compare the results with the results from the depth analysis. It was not too conclusive I am afraid, and trying to understand the results I came across independent component analysis (ICA). So .. Ive run it on a combination of principal components and it seems that certainly we have 2 independent components that keep popping up (if I can use this expression) when we run the analysis with 2, 3, or 4 components. So I guess these 2 components are the strongest ones . If I can say so. My next question is . How can I relate these 2 components to the initial data??? Ive plotted each component and if I add a loess line to each, visually it seems that one independent component is an unknown function of rugosity while the other component is an unknown function of reef geometry and depth. But, of course, I would like something more than a visual similarity. Also the 2 independent components seem to split the data in 3 classes, rather than 2, as the analysis of the depth data suggested. Looking back at the depth histogram it is obvious that there are some data that actually are not quite modeled by the 2 mixing functions I came up with. These data correspond to the deepest patch reefs, a category clearly singled out by ICA classification. The bottom line is that I am trying to understand what each independent component tells me about the patch reefs and how I can relate that to the patch reef morphometrics, biology, other factors that impact some reefs but not others, etc. If you have any clarifying thoughts or if you know about any other literature about the subject that can help (except articles that deal with ICA and image analysis or wave form data) I will really appreciate. Thank you very much for your consideration, Monica _ The average US Credit Score is 675. The cost to see yours: $0 by Experian. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Is this a bug?
I have found a strange ifelse behaviour (I think) This works: ifelse(T,1+1,1+2) [1] 2 ifelse(F,1+1,1+2) [1] 3 Maybe I missed something about R internals, but why ifelse(T,print(hello),print(goodbye)) [1] hello [1] hello ifelse(F,print(hello),print(goodbye)) [1] goodbye [1] goodbye values are returned two times? I'm not sure: if it's a bug I'll post it immediately Thank You Luca Version: platform = i486-pc-linux-gnu arch = i486 os = linux-gnu system = i486, linux-gnu status = major = 2 minor = 4.1 year = 2006 month = 12 day = 18 svn rev = 40228 language = R version.string = R version 2.4.1 (2006-12-18) Locale: LC_CTYPE=it_IT.UTF-8;LC_NUMERIC=C;LC_TIME=it_IT.UTF-8;LC_COLLATE=it_IT.UTF-8;LC_MONETARY=it_IT.UTF-8;LC_MESSAGES=it_IT.UTF-8;LC_PAPER=it_IT.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=it_IT.UTF-8;LC_IDENTIFICATION=C Search Path: .GlobalEnv, package:MASS, package:utils, package:stats, package:graphics, package:grDevices, package:methods, Autoloads, package:base __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Is this a bug?
one is returned value, the other one is the result from print t0 - ifelse(T, print(h), print(e)) [1] h t0 [1] h HTH, weiwei On 4/17/07, Luca Braglia [EMAIL PROTECTED] wrote: I have found a strange ifelse behaviour (I think) This works: ifelse(T,1+1,1+2) [1] 2 ifelse(F,1+1,1+2) [1] 3 Maybe I missed something about R internals, but why ifelse(T,print(hello),print(goodbye)) [1] hello [1] hello ifelse(F,print(hello),print(goodbye)) [1] goodbye [1] goodbye values are returned two times? I'm not sure: if it's a bug I'll post it immediately Thank You Luca Version: platform = i486-pc-linux-gnu arch = i486 os = linux-gnu system = i486, linux-gnu status = major = 2 minor = 4.1 year = 2006 month = 12 day = 18 svn rev = 40228 language = R version.string = R version 2.4.1 (2006-12-18) Locale: LC_CTYPE=it_IT.UTF-8;LC_NUMERIC=C;LC_TIME=it_IT.UTF-8;LC_COLLATE=it_IT.UTF-8;LC_MONETARY=it_IT.UTF-8;LC_MESSAGES=it_IT.UTF-8;LC_PAPER=it_IT.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=it_IT.UTF-8;LC_IDENTIFICATION=C Search Path: .GlobalEnv, package:MASS, package:utils, package:stats, package:graphics, package:grDevices, package:methods, Autoloads, package:base __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Weiwei Shi, Ph.D Research Scientist GeneGO, Inc. Did you always know? No, I did not. But I believed... ---Matrix III __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Is this a bug?
On 4/17/07, Luca Braglia [EMAIL PROTECTED] wrote: I have found a strange ifelse behaviour (I think) Don't you think it is rather consistent behavior? ifelse(T,1+1,1+2) [1] 2 ifelse(F,1+1,1+2) [1] 3 ifelse(T,hello,goodbye) [1] hello ifelse(F,hello,goodbye) [1] goodbye ifelse(T,print(hello),print(goodbye)) [1] hello [1] hello ifelse(F,print(hello),print(goodbye)) [1] goodbye [1] goodbye ifelse(T,print(1+1),print(1+2)) [1] 2 [1] 2 ifelse(F,print(1+1),print(1+2)) [1] 3 [1] 3 This works: ifelse(T,1+1,1+2) [1] 2 ifelse(F,1+1,1+2) [1] 3 Maybe I missed something about R internals, but why ifelse(T,print(hello),print(goodbye)) [1] hello [1] hello ifelse(F,print(hello),print(goodbye)) [1] goodbye [1] goodbye values are returned two times? I'm not sure: if it's a bug I'll post it immediately Thank You Luca Version: platform = i486-pc-linux-gnu arch = i486 os = linux-gnu system = i486, linux-gnu status = major = 2 minor = 4.1 year = 2006 month = 12 day = 18 svn rev = 40228 language = R version.string = R version 2.4.1 (2006-12-18) Locale: LC_CTYPE=it_IT.UTF-8;LC_NUMERIC=C;LC_TIME=it_IT.UTF-8;LC_COLLATE=it_IT.UTF-8;LC_MONETARY=it_IT.UTF-8;LC_MESSAGES=it_IT.UTF-8;LC_PAPER=it_IT.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=it_IT.UTF-8;LC_IDENTIFICATION=C Search Path: .GlobalEnv, package:MASS, package:utils, package:stats, package:graphics, package:grDevices, package:methods, Autoloads, package:base __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Greatest common divisor of two numbers
Thank you so much All. It works. Abou == AbouEl-Makarim Aboueissa, Ph.D. Assistant Professor of Statistics Department of Mathematics Statistics University of Southern Maine 96 Falmouth Street P.O. Box 9300 Portland, ME 04104-9300 Tel: (207) 228-8389 Email: [EMAIL PROTECTED] [EMAIL PROTECTED] Office: 301C Payson Smith __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] help comparing two median with R
Thomas Lumley wrote: There actually is an exact test for the median that does not assume a location shift: dichotomize your data at the pooled median to get a 2x2 table of above/below median by group, and do Fisher's exact test on the table. Fascinating. But can one be sure that the fisher test actually has the correct distribution in this setup? I feel somewhat unconvinced. It is fairly clear that the test is not independent of the joint sample median for instance (think two normals, one with a very small variance). This is almost never useful (because it doesn't come with an interval estimate), but is interesting because it (and the generalizations to other quantiles) is the only exactly distribution-free location test that does not have the 'non-transitivity' problem of the Mann-Whitney U test. I believe this median test is attributed to Mood, but I have not seen the primary source. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Is this a bug?
On 17/04/07 - 14:59, Roland Rau wrote: On 4/17/07, Luca Braglia [EMAIL PROTECTED] wrote: ifelse(T,1+1,1+2) [1] 2 ifelse(F,1+1,1+2) [1] 3 ifelse(T,hello,goodbye) [1] hello ifelse(F,hello,goodbye) [1] goodbye ifelse(T,print(hello),print(goodbye)) [1] hello [1] hello ifelse(F,print(hello),print(goodbye)) [1] goodbye [1] goodbye ifelse(T,print(1+1),print(1+2)) [1] 2 [1] 2 ifelse(F,print(1+1),print(1+2)) [1] 3 [1] 3 Thank you , Weiwei and Roland, all right now: I was thinking wrong! bye Luca __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Tcltk
I suspect tcl's own version of 'access', but can you please confirm that
this still happens under 'Run as Administrator', assuming 'C:\Program' is
a system area in Swedish Windows Vista?
I will be able to take a closer look, but not before 2.5.0 (which is in
code freeze and I have limited acccess to a Vista machine).
On Tue, 17 Apr 2007, Sofia Wikström wrote:
I have problems with Tcl/Tk in R 2.4.1, when running it on Windows Vista
(see error message below).
Regards, Sofia
library(tcltk)
Loading Tcl/Tk interface ... Error in fun(...) : Can't find a usable
init.tcl in the following directories:
{C:\Program\R\R-2.4.1/Tcl/lib/tcl8.4}
{C:\Program\R\R-2.4.1/Tcl/lib/tcl8.4} C:/Program/R/R-2.4.1/Tcl/lib/tcl8.4
C:/Program/R/R-2.4.1/Tcl/lib/tcl8.4
This probably means that Tcl wasn't installed properly.
Error: .onLoad failed in 'loadNamespace' for 'tcltk'
Error: package/namespace load failed for 'tcltk'
_
Sofia Wikström, PhD
AquaBiota Water Research
Svante Arrhenius väg 21A, SE-104 05 Stockholm, Sweden
Phone: (+46) 8 16 10 07
[EMAIL PROTECTED]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Brian D. Ripley, [EMAIL PROTECTED]
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax: +44 1865 272595__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Question of the impact of the pilot experiment on the overal statistic interpretation of the subsequent work
Hi, I have a question regarding the impact of the pilot experiment on the overall statistic interpretation of the subsequent work. The context is as following: In a lab there are one Professor and three graduate students A, B, and C. They are working on analysis of some disease to discover genes differentiate + and - categories like disease or non-disease. Number of samples of (+) is about m while that of the samples of (-) is n. Both m, n are sufficiently large, e.g. bigger than 100. Pilot experiment: In order to save efforts and resources, the Professor decided to pool the samples in each category with equal amount such that he got only two pooled samples of (+) and (-). His argument is that if there is no difference in the pooled samples then he would decide to abandon the project. Graduate student A did a microarray analysis of the pooled (+) and (-) and found gene X, Y, Z have fold of change bigger than 100. Professor thought this was interesting and encouraging based upon his biological insight of gene X, Y and Z and the potential disease link of these genes. (1) He asked graduate student B to do a protein analysis, using a different technique (western blot), of all the original samples (m, n) and found gene Y is truly differential. Based upon the protein analysis data, graduate student A calculated P value using t test to describe the statistic significance of gene Y differentiating (+) and (-) categories. (2) Simultaneously, he also asked graduate student C to do a full scale microarray experiment using all m, n samples individually. It is a very laborious work but graduate student C finished everything and using some off the shelf microarray statistical packages, he calculated and found gene Y, Z and another un-identified gene W to be statistically significant. He calculated the false discovery rate and P value of these genes differentiating (+) and (-) categories. The professor presented his students A, B, C's work including the calculated statistics in a conference. In the audience, statistician D commented that professor has made a mistake here: because he is using the SAME samples, whether pooled or individual, in both the pilot and subsequent experiments, statistically the professor is cheating and his students' calculated statistics are no longer valid. Can statisticians in this mailing list comment on this story? One thing I want to emphasize here is that nobody disputes that it is highly critical to use a different set of samples to validate the discoveries. The question here is that contingent on the pilot experiment of the pooled samples, whether the subsequent full scale experiments using the SAME samples can yield meaningful statistics to describe the differences of the discovered features. Thanks. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Dealing with data frame column names beginning with a numeric
On Mon, 16-Apr-2007 at 10:12PM +0100, Prof Brian Ripley wrote:
| On Mon, 16 Apr 2007, Duncan Murdoch wrote:
| The name change happens in the conversion to a dataframe, so why not
| change the name afterwards? That is:
|
| df - data.frame(mat)
| names(df)[3] - 5T
| boxplot(df, main=blah blah blah)
|
| Or use
|
| boxplot(as.data.frame(mat))
|
| which seems more natural than data.frame(mat, check.names=FALSE) (which
| also does the job) or even data.frame(mat).
boxplot(data.frame(mat)) is how it's described in the help, but it
produces the problem that began my discussion (i.e. prepends an X to
the name for a dataframe column name).
However, boxplot(as.data.frame(mat)) as suggested by Brian works fine.
So I needed only three keystrokes to get round the problem.
Thank you all.
--
~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.
___Patrick Connolly
{~._.~} Great minds discuss ideas
_( Y )_Middle minds discuss events
(:_~*~_:)Small minds discuss people
(_)-(_) . Anon
~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] extracting intercept from ppr fit
Sorry for triple-posting : I seem to have a problem w/ my mail client. Hi, Is there a way, documented or not, to extract the intercept term (the alpha_0 the MASS book) from a ppr() (Projection Persuit Regression) fit? Thanks, Vadim ## Example: n - 1000 data - data.frame(x= rnorm (n), y= rnorm (n)) a - 10 data$z - evalq(a + atan (x + y) + rnorm (n), data) data.ppr - ppr(z ~ x + y, data=data, nterms =1) ## how to extract a = 10 from data.ppr? [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Runing R in a bash script
Hello!
I am having issues trying to plot to a ong (or jpg) when the R-code in a
bash script is executed from cron.
I can generate a pdf file, but when I try to write to a png, the file is
created, but nothing is written. If I execute the bash script from my
console, everything works file. Any ideas?
In my cron I have SHELL=/bin/bash - otherwise /bin/shell is used and the
folowing enery, so example is executed every minute
* * * * * [path]/example.sh
I am running
R version 2.4.1 (2006-12-18)
Here's a minimal example - two files one R-script ('example.r') and one
bash-script ('example.sh')
example.r
# Example R-script
x - c(1:10)
y - x^2
png(file=example2.png)
#pdf(file=example2.pdf)
plot(x,y)
graphics.off()
example.sh
#/bin/bash
#
# Hello world is written to exhotext every time cron executes this script
echo Hello world echotext
# This works, but not when executed from cron
n=`R --save example.r`
# using exec as in `exec R --save example.r` dosent work with cron either
# This also works, but nothing is written to the png when executed from cron
R --save RSCRIPT
x - c(1:10)
y - x^2
png(file=example2.png)
#pdf(file=example2.pdf)
plot(x,y)
graphics.off()
#dev.off() dosent work at all when executed from cron
RSCRIPT
Thanks
Ulrik
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Question of the impact of the pilot experiment on the overal statistic interpretation of the subsequent work
Bruce, Far below you ask for comment on an issue of statistical interpretation of a collection of biological experiments. Your university (judging from your email handle) has one of the best statistics departments in the world and some of the best biostatisticians in the galaxy. You would do far better (than posting your query here) to take your issue up with one of them. In a face-to-face meeting with one of them you will get a much better analysis and discussion of the issues than you could hope for from a list like this, notwithstanding that some statisticians sometimes provide thoughtful answers to posts asking for statistical help. Chuck On Tue, 17 Apr 2007, Bruce Ling wrote: Hi, I have a question regarding the impact of the pilot experiment on the overall statistic interpretation of the subsequent work. The context is as following: In a lab there are one Professor and three graduate students A, B, and C. They are working on analysis of some disease to discover genes differentiate + and - categories like disease or non-disease. Number of samples of (+) is about m while that of the samples of (-) is n. Both m, n are sufficiently large, e.g. bigger than 100. Pilot experiment: In order to save efforts and resources, the Professor decided to pool the samples in each category with equal amount such that he got only two pooled samples of (+) and (-). His argument is that if there is no difference in the pooled samples then he would decide to abandon the project. Graduate student A did a microarray analysis of the pooled (+) and (-) and found gene X, Y, Z have fold of change bigger than 100. Professor thought this was interesting and encouraging based upon his biological insight of gene X, Y and Z and the potential disease link of these genes. (1) He asked graduate student B to do a protein analysis, using a different technique (western blot), of all the original samples (m, n) and found gene Y is truly differential. Based upon the protein analysis data, graduate student A calculated P value using t test to describe the statistic significance of gene Y differentiating (+) and (-) categories. (2) Simultaneously, he also asked graduate student C to do a full scale microarray experiment using all m, n samples individually. It is a very laborious work but graduate student C finished everything and using some off the shelf microarray statistical packages, he calculated and found gene Y, Z and another un-identified gene W to be statistically significant. He calculated the false discovery rate and P value of these genes differentiating (+) and (-) categories. The professor presented his students A, B, C's work including the calculated statistics in a conference. In the audience, statistician D commented that professor has made a mistake here: because he is using the SAME samples, whether pooled or individual, in both the pilot and subsequent experiments, statistically the professor is cheating and his students' calculated statistics are no longer valid. Can statisticians in this mailing list comment on this story? One thing I want to emphasize here is that nobody disputes that it is highly critical to use a different set of samples to validate the discoveries. The question here is that contingent on the pilot experiment of the pooled samples, whether the subsequent full scale experiments using the SAME samples can yield meaningful statistics to describe the differences of the discovered features. Thanks. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Charles C. Berry(858) 534-2098 Dept of Family/Preventive Medicine E mailto:[EMAIL PROTECTED] UC San Diego http://biostat.ucsd.edu/~cberry/ La Jolla, San Diego 92093-0901 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Runing R in a bash script
Ulrik Stervbo wrote:
Hello!
I am having issues trying to plot to a ong (or jpg) when the R-code in a
bash script is executed from cron.
I can generate a pdf file, but when I try to write to a png, the file is
created, but nothing is written. If I execute the bash script from my
console, everything works file. Any ideas?
In my cron I have SHELL=/bin/bash - otherwise /bin/shell is used and the
folowing enery, so example is executed every minute
* * * * * [path]/example.sh
I am running
R version 2.4.1 (2006-12-18)
Here's a minimal example - two files one R-script ('example.r') and one
bash-script ('example.sh')
example.r
# Example R-script
x - c(1:10)
y - x^2
png(file=example2.png)
#pdf(file=example2.pdf)
plot(x,y)
graphics.off()
example.sh
#/bin/bash
#
# Hello world is written to exhotext every time cron executes this script
echo Hello world echotext
# This works, but not when executed from cron
n=`R --save example.r`
# using exec as in `exec R --save example.r` dosent work with cron
either
# This also works, but nothing is written to the png when executed
from cron
R --save RSCRIPT
x - c(1:10)
y - x^2
png(file=example2.png)
#pdf(file=example2.pdf)
plot(x,y)
graphics.off()
#dev.off() dosent work at all when executed from cron
RSCRIPT
The png() device requires an X server for the image rendering. You might
be able to get away with exporting the DISPLAY environment variable
export DISPLAY=:0.0 # try and connect to X server on display 0.0
within your script, but it will only work if the script is executed by
the same user as is running the X server, *and* the X server is running
at the time the script is executed.
There are a handful of packages that will create a png without the
presence of an X server, and I'm partial to Cairo (since I've done some
work on it). You can install the latest version like this:
install.packages(Cairo,,'http://rforge.net/',type='source')
Cairo can also outputs nice pdf's with embedded fonts... useful if you
want to embed high-quality OpenType or TrueType fonts.
Best,
Jeff
--
http://biostat.mc.vanderbilt.edu/JeffreyHorner
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Extracting approximate Wald test (Chisq) from coxph(..frailty)
Dear list, I need to extract the approximate Wald test (Chisq) so that I can put it in a loop. str seemed like a great idea, but I cannot seem to find the approximate Wald test for frailty (in the example data below: 17.89 and its p-value 0.12000) there. I cannot seem to find it in capture.output either as numeric form. Do I need to modify some given values? If yes, please give me a clue for the example: library(survival) kfitm1-coxph(formula = Surv(time, status) ~ age + sex +disease + frailty(id, dist = gauss), data = kidney) str(kfitm1) capture.output( print(kfitm1) ) Mohammad Ehsanul Karim (R - 2.3.1 on windows) wildscop at yahoo dot com Institute of Statistical Research and Training University of Dhaka On Tue, 17 Apr 2007, Mohammad Ehsanul Karim wrote: You _can_ use tmp - capture.output( print( your example ) ) and then further process tmp. A _better_ solution for most purposes is to look at the object produced by coxph() and figure out how to calculate the Wald statistic from that object. See ?coxph.object and ?str Another tactic is to look at how print.coxph() does its work and use the code in it to produce just the output you desire. Look at page( survival:::print.coxph, print ) Assign the output of coxph to some object, and use the $ extractor function to obtain what you need. ie: rtfm - coxph(formula = Surv(time, status) ~ age + sex + disease + frailty(id, dist = gauss), data = kidney) Age - coef(rtfm)[age] OR Sex - rtfm$coef[sex] Mohammad Ehsanul Karim wrote: Dear List, How do I extract the approximate Wald test for the frailty (in the following example 17.89 value)? What about the P-values, other Chisq, DF, se(coef) and se2? How can they be extracted? ## kfitm1 Call: coxph(formula = Surv(time, status) ~ age + sex + disease + frailty(id, dist = gauss), data = kidney) coef se(coef) age0.00489 0.0150 sex -1.69703 0.4609 diseaseGN 0.17980 0.5447 diseaseAN 0.39283 0.5447 diseasePKD-1.13630 0.8250 frailty(id, dist = gauss se2Chisq DF age 0.0106 0.11 1.0 sex 0.3617 13.56 1.0 diseaseGN 0.3927 0.11 1.0 diseaseAN 0.3982 0.52 1.0 diseasePKD0.6173 1.90 1.0 frailty(id, dist = gauss17.89 12.1 p age 0.74000 sex 0.00023 diseaseGN 0.74000 diseaseAN 0.47000 diseasePKD0.17000 frailty(id, dist = gauss 0.12000 Iterations: 6 outer, 30 Newton-Raphson Variance of random effect= 0.493 Degrees of freedom for terms= 0.5 0.6 1.7 12.1 Likelihood ratio test=47.5 on 14.9 df, p=2.82e-05 n= 76 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Extracting approximate Wald test (Chisq) fromcoxph(..frailty)
Hi Mohammed,
Here's one way to get the information you need.
First I examined the output of your coxph() call:
library(survival)
kfitm1-coxph(formula = Surv(time, status) ~ age +
+ sex +disease + frailty(id, dist = gauss),
+ data = kidney)
class(kfitm1)
[1] coxph.penal coxph
attributes(kfitm1)
$names
[1] coefficients var var2
[4] loglikiter linear.predictors
[7] residuals means method
[10] frail fvar df
[13] df2 penalty pterms
[16] assign2 history coxlist1
[19] printfun n terms
[22] assignwald.test y
[25] formula call
$class
[1] coxph.penal coxph
So the returned object is actually of class coxph.penal,
not coxph. Thus you'll want to look into the functions
survival:::summary.coxph.penal
survival:::print.coxph.penal
I checked the wald piece of the returned object,
that's not what you needed.
kfitm1$wald.test
[1] 14.96798
I see that what you want is in the temp
matrix composed in the
survival:::summary.coxph.penal
function.
Make your own copy of the function
my.summary.coxph.penal - survival:::summary.coxph.penal
Edit the function and return the temp matrix
(see modified my.summary.coxph.penal function below,
I added a return list including the temp matrix.)
kfitm1.my.summary - my.summary.coxph.penal(kfitm1)
Call:
coxph(formula = Surv(time, status) ~ age + sex + disease + frailty(id,
dist = gauss), data = kidney)
n= 76
coef se(coef) se2Chisq DF p
age0.00489 0.0150 0.0106 0.11 1.0 0.74000
sex -1.69703 0.4609 0.3617 13.56 1.0 0.00023
diseaseGN 0.17980 0.5447 0.3927 0.11 1.0 0.74000
diseaseAN 0.39283 0.5447 0.3982 0.52 1.0 0.47000
diseasePKD-1.13630 0.8250 0.6173 1.90 1.0 0.17000
frailty(id, dist = gauss 17.89 12.1 0.12000
exp(coef) exp(-coef) lower .95 upper .95
age1.005 0.9950.9759 1.035
sex0.183 5.4580.0742 0.452
diseaseGN 1.197 0.8350.4116 3.481
diseaseAN 1.481 0.6750.5093 4.308
diseasePKD 0.321 3.1150.0637 1.617
Iterations: 6 outer, 30 Newton-Raphson
Variance of random effect= 0.493
Degrees of freedom for terms= 0.5 0.6 1.7 12.1
Rsquare= 0.465 (max possible= 0.993 )
Likelihood ratio test= 47.5 on 14.9 df, p=2.82e-05
Wald test= 15.0 on 14.9 df, p=0.446
attributes(kfitm1.my.summary)
$names
[1] temp tmp
kfitm1.my.summary$temp
coef se(coef) se2 Chisq DF p
age0.00489 0.0150 0.0106 0.11 1.0 0.74000
sex -1.69703 0.4609 0.3617 13.56 1.0 0.00023
diseaseGN 0.17980 0.5447 0.3927 0.11 1.0 0.74000
diseaseAN 0.39283 0.5447 0.3982 0.52 1.0 0.47000
diseasePKD-1.13630 0.8250 0.6173 1.90 1.0 0.17000
frailty(id, dist = gauss17.89 12.1 0.12000
class(kfitm1.my.summary$temp)
[1] matrix
kfitm1.my.summary$temp[grep(frailty,
dimnames(kfitm1.my.summary$temp)[[1]]), Chisq]
[1] 17.89
So you can get the information you need from the returned
temp matrix as above.
There are many other ways to do this, but the above
ideas can get you going.
Hope this helps
my.summary.coxph.penal -
function (object, conf.int = 0.95, scale = 1, terms = FALSE,
maxlabel = 25, digits = max(options()$digits - 4, 3), ...)
{
if (!is.null(object$call)) {
cat(Call:\n)
dput(object$call)
cat(\n)
}
if (!is.null(object$fail)) {
cat( Coxreg failed., object$fail, \n)
return()
}
savedig - options(digits = digits)
on.exit(options(savedig))
omit - object$na.action
if (length(omit))
cat( n=, object$n, (, naprint(omit), )\n, sep = )
else cat( n=, object$n, \n)
coef - object$coef
if (length(coef) == 0 length(object$frail) == 0)
stop(Penalized summary function can't be used for a null model)
if (length(coef) 0) {
nacoef - !(is.na(coef))
coef2 - coef[nacoef]
if (is.null(coef) | is.null(object$var))
stop(Input is not valid)
se - sqrt(diag(object$var))
}
pterms - object$pterms
nterms - length(pterms)
npenal - sum(pterms 0)
print.map - rep(0, nterms)
if (!is.null(object$printfun)) {
temp - unlist(lapply(object$printfun, is.null))
print.map[pterms 0] - (1:npenal) * (!temp)
}
print1 - NULL
pname1 - NULL
if (is.null(object$assign2))
alist - object$assign[-1]
else alist -
Re: [R] help comparing two median with R
On Tue, 17 Apr 2007, Frank E Harrell Jr wrote: The points that Thomas and Brian have made are certainly correct, if one is truly interested in testing for differences in medians or means. But the Wilcoxon test provides a valid test of x y more generally. The test is consonant with the Hodges-Lehmann estimator: the median of all possible differences between an X and a Y. Yes, but there is no ordering of distributions (taken one at a time) that agrees with the Wilcoxon two-sample test, only orderings of pairs of distributions. The Wilcoxon test provides a test of xy if it is known a priori that the two distributions are stochastically ordered, but not under weaker assumptions. Otherwise you can get xyzx. This is in contrast to the t-test, which orders distributions (by their mean) whether or not they are stochastically ordered. Now, it is not unreasonable to say that the problems are unlikely to occur very often and aren't worth worrying too much about. It does imply that it cannot possibly be true that there is any summary of a single distribution that the Wilcoxon test tests for (and the same is true for other two-sample rank tests, eg the logrank test). I know Frank knows this, because I gave a talk on it at Vanderbilt, but most people don't know it. (I thought for a long time that the Wilcoxon rank-sum test was a test for the median pairwise mean, which is actually the R-estimator corresponding to the *one*-sample Wilcoxon test). -thomas __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Runing R in a bash script
Or see png2() in R.utils, which imitates png() but uses bitmap(),
which in turn uses postscript-to-png via ghostscript. BTW, personally
I think PNGs generated via bitmap() look way better than the ones
generated via png().
/Henrik
On 4/17/07, Jeffrey Horner [EMAIL PROTECTED] wrote:
Ulrik Stervbo wrote:
Hello!
I am having issues trying to plot to a ong (or jpg) when the R-code in a
bash script is executed from cron.
I can generate a pdf file, but when I try to write to a png, the file is
created, but nothing is written. If I execute the bash script from my
console, everything works file. Any ideas?
In my cron I have SHELL=/bin/bash - otherwise /bin/shell is used and the
folowing enery, so example is executed every minute
* * * * * [path]/example.sh
I am running
R version 2.4.1 (2006-12-18)
Here's a minimal example - two files one R-script ('example.r') and one
bash-script ('example.sh')
example.r
# Example R-script
x - c(1:10)
y - x^2
png(file=example2.png)
#pdf(file=example2.pdf)
plot(x,y)
graphics.off()
example.sh
#/bin/bash
#
# Hello world is written to exhotext every time cron executes this script
echo Hello world echotext
# This works, but not when executed from cron
n=`R --save example.r`
# using exec as in `exec R --save example.r` dosent work with cron
either
# This also works, but nothing is written to the png when executed
from cron
R --save RSCRIPT
x - c(1:10)
y - x^2
png(file=example2.png)
#pdf(file=example2.pdf)
plot(x,y)
graphics.off()
#dev.off() dosent work at all when executed from cron
RSCRIPT
The png() device requires an X server for the image rendering. You might
be able to get away with exporting the DISPLAY environment variable
export DISPLAY=:0.0 # try and connect to X server on display 0.0
within your script, but it will only work if the script is executed by
the same user as is running the X server, *and* the X server is running
at the time the script is executed.
There are a handful of packages that will create a png without the
presence of an X server, and I'm partial to Cairo (since I've done some
work on it). You can install the latest version like this:
install.packages(Cairo,,'http://rforge.net/',type='source')
Cairo can also outputs nice pdf's with embedded fonts... useful if you
want to embed high-quality OpenType or TrueType fonts.
Best,
Jeff
--
http://biostat.mc.vanderbilt.edu/JeffreyHorner
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Data Manipulation using R
Dear Friends, I have data set with around 220,000 rows and 17 columns. One of the columns is an id variable which is grouped from 1000 through 9000. I need to perform the following operations. 1) Remove all the observations with id's between 6000 and 6999 I tried using this method. remdat1 - subset(data, ID6000) remdat2 - subset(data, ID=7000) donedat - rbind(remdat1, remdat2) I check the last and first entry and found that it did not have ID values 6000. Therefore I think that this might be correct, but is this the most efficient way of doing this? 2) I need to remove observations within columns 3, 4, 6 and 8 when they are negative. For instance if the number in column 3 is -4, then I need to delete the entire observation. Can somebody help me with this too. Thank and Regards Anup - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Einladung - Invito
(Versione italiana più sotto) Sehr geehrte Damen und Herren Wir haben bemerkt, dass viele Firmenemailadressen keine HTML-Emails empfangen dürfen. Aus diesem Grund senden wir Ihnen unsere Einladung auch in Textformat zu, damit Sie diese rechtzeitig erhalten, denn die Anmeldefrist läuft bald ab. Hiermit laden wir Sie herzlich ein unsere exklusive Veranstaltung zu besuchen: Der bulgarische Finanzmarkt - Gelegenheiten und Perspektiven in Zürich am 9. Mai 2007, 9:00-14:00 Uhr, Hotel Baur au Lac Die Information ist sowohl für institutionelle als auch für private Anleger wichtig, da sie dadurch einen Insider-Einblick in den boomenden Finanzmarkt Bulgariens erhalten. Die Veranstaltung eröffnet neue Perspektiven und bietet bis jetzt nicht wahrgenommenen Chancen für den Profi-Investor. Für alle Details öffnen Sie bitte diesen Link in Ihrem Browser: http://www.swisscham-bulgaria.org/index.php?id=215tx_rlmpeventdb_pi1[showUid]=3no_cache=1 Sagen Sie uns Ihre Meinung über diese Veranstaltung in unserer anonymen Umfrage hier: http://www.swisscham-bulgaria.org/index.php?id=290 Wenn Sie keine weiteren Einladungen von der WKSB erhalten möchten, folgen Sie bitte diesem Link: http://www.swisscham-bulgaria.org/cgi-bin/acmelist/l.cgi?subid=18uss=1[EMAIL PROTECTED] Mit freundlichen Grüssen, Biliana Monova, Sekretariat Versione italiana Gentili signore, egregi signori, Abbiamo notato, che tante ditte non ricevano le e-mail in HTML. Per questo motivo inviamo il nostro invito in formato testo, per farlo diventare ricevibile e per arrivarvi in tempo utile, visto che i tempi per la registrazione si stanno esaurendo. Con il presente vinvitiamo a partecipare al nostro evento esclusivo: Il mercato monetario bulgaro - occasioni e prospettive 10 maggio 2007, 9:00-14:00 nello Splendid Royal Hotel di Lugano Linformazione è estremamente utile sia per gli investitori istituzionali, sia per quelli di carattere privato. In questo modo potete avere uno sguardo dallinterno nel prospero mercato monetario della Bulgaria. Levento apre nuove prospettive e propone agli investitori professionali opportunità fin ora sconosciute o sottovalutate. Per ottenere i dettagli riguardo allevento, vi preghiamo di aprire questo link nel Vostro Browser: http://www.swisscham-bulgaria.org/index.php?id=215L=5tx_rlmpeventdb_pi1[showUid]=4no_cache=1 Comunicateci il Vostro prezioso parere riguardo allevento attraverso la qui presente inchiesta anonima: http://www.swisscham-bulgaria.org/index.php?id=290L=5 Se volete ottenere un maggior numero dinviti dalla Camera di Economia SvizzeraBulgaria, seguite questo link: http://www.swisscham-bulgaria.org/cgi-bin/acmelist/l.cgi?subid=18uss=1[EMAIL PROTECTED] Cordiali saluti, Biliana Monova, segretario __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to estimate dose from respond given drc package result
Try this: library(drc) # get dose corresponding to response (y = 2) model1 - multdrc(rootl ~ conc, data = ryegrass) f - function(x, y) y - predict(model1, data.frame(conc = x))[1] uniroot(f, c(0, 10), y = 2) On 4/17/07, John McNeil [EMAIL PROTECTED] wrote: Dear all, I can use the very nice drc package (multdrc()) to model and plot a dataframe containing dose and response values. I can also use predict.drc() to yield response values given a dose. I need to do the opposite, estimate a dose given the response. The general predict documentation seems to say that this is possible, but it does not appear that predict.drc has that capability. This makes sense because it is easy to compute y values from x values, given the fitted equation, but it is hard to do the opposite. It seems that I need to do some kind of successive approximation. Is there a package that takes a multdrc generated model and estimates dose given response? Thank you all for your time, John __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Data Manipulation using R
On Apr 17, 2007, at 8:03 PM, Anup Nandialath wrote: Dear Friends, I have data set with around 220,000 rows and 17 columns. One of the columns is an id variable which is grouped from 1000 through 9000. I need to perform the following operations. 1) Remove all the observations with id's between 6000 and 6999 I tried using this method. remdat1 - subset(data, ID6000) remdat2 - subset(data, ID=7000) donedat - rbind(remdat1, remdat2) I check the last and first entry and found that it did not have ID values 6000. Therefore I think that this might be correct, but is this the most efficient way of doing this? The rbind is a bit unnecessary probably. I think all you are missing for both questions is the or operator, |. ( ?| ) Simply: donedat - subset(data, ID 6000 | ID =7000) would do for this. Not sure about efficiency, but if the code is fast as it stands I wouldn't worry too much about it. 2) I need to remove observations within columns 3, 4, 6 and 8 when they are negative. For instance if the number in column 3 is -4, then I need to delete the entire observation. Can somebody help me with this too. The following should do it (untested, not sure if it would handle NA's): toremove - data[,3] 0 | data[,4] 0 | data[,6] 0 | data[,8] 0 data[!toremove,] If you want more columns than those 4, then we could perhaps look for a better line than the first line above. Thank and Regards Anup Haris Skiadas Department of Mathematics and Computer Science Hanover College __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to estimate dose from respond given drc package result
Wow, thanks, that does it perfectly John On Apr 17, 2007, at 5:51 PM, Gabor Grothendieck wrote: Try this: library(drc) # get dose corresponding to response (y = 2) model1 - multdrc(rootl ~ conc, data = ryegrass) f - function(x, y) y - predict(model1, data.frame(conc = x))[1] uniroot(f, c(0, 10), y = 2) On 4/17/07, John McNeil [EMAIL PROTECTED] wrote: Dear all, I can use the very nice drc package (multdrc()) to model and plot a dataframe containing dose and response values. I can also use predict.drc() to yield response values given a dose. I need to do the opposite, estimate a dose given the response. The general predict documentation seems to say that this is possible, but it does not appear that predict.drc has that capability. This makes sense because it is easy to compute y values from x values, given the fitted equation, but it is hard to do the opposite. It seems that I need to do some kind of successive approximation. Is there a package that takes a multdrc generated model and estimates dose given response? Thank you all for your time, John __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting- guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] R-OSX error: 'memory not mapped'
Hi, I often get a following error with R *** caught segfault *** address 0x78807e00, cause 'memory not mapped' Possible actions: 1: abort (with core dump) 2: normal R exit 3: exit R without saving workspace 4: exit R saving workspace Selection: The system is OSX 4.9 and R-version 2.4.1. Is there something to d0? Atte Tenkanen __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Runing R in a bash script
As I har problems installing the Cairo package, I went for Henriks solution
- and it works almost perfect. I would like to have been able to generate
transparent png.
Thanks for the help
Ulrik
On 18/04/07, Henrik Bengtsson [EMAIL PROTECTED] wrote:
Or see png2() in R.utils, which imitates png() but uses bitmap(),
which in turn uses postscript-to-png via ghostscript. BTW, personally
I think PNGs generated via bitmap() look way better than the ones
generated via png().
/Henrik
On 4/17/07, Jeffrey Horner [EMAIL PROTECTED] wrote:
Ulrik Stervbo wrote:
Hello!
I am having issues trying to plot to a ong (or jpg) when the R-code
in a
bash script is executed from cron.
I can generate a pdf file, but when I try to write to a png, the file
is
created, but nothing is written. If I execute the bash script from my
console, everything works file. Any ideas?
In my cron I have SHELL=/bin/bash - otherwise /bin/shell is used and
the
folowing enery, so example is executed every minute
* * * * * [path]/example.sh
I am running
R version 2.4.1 (2006-12-18)
Here's a minimal example - two files one R-script ('example.r') and
one
bash-script ('example.sh')
example.r
# Example R-script
x - c(1:10)
y - x^2
png(file=example2.png)
#pdf(file=example2.pdf)
plot(x,y)
graphics.off()
example.sh
#/bin/bash
#
# Hello world is written to exhotext every time cron executes this
script
echo Hello world echotext
# This works, but not when executed from cron
n=`R --save example.r`
# using exec as in `exec R --save example.r` dosent work with cron
either
# This also works, but nothing is written to the png when executed
from cron
R --save RSCRIPT
x - c(1:10)
y - x^2
png(file=example2.png)
#pdf(file=example2.pdf)
plot(x,y)
graphics.off()
#dev.off() dosent work at all when executed from cron
RSCRIPT
The png() device requires an X server for the image rendering. You might
be able to get away with exporting the DISPLAY environment variable
export DISPLAY=:0.0 # try and connect to X server on display 0.0
within your script, but it will only work if the script is executed by
the same user as is running the X server, *and* the X server is running
at the time the script is executed.
There are a handful of packages that will create a png without the
presence of an X server, and I'm partial to Cairo (since I've done some
work on it). You can install the latest version like this:
install.packages(Cairo,,'http://rforge.net/',type='source')
Cairo can also outputs nice pdf's with embedded fonts... useful if you
want to embed high-quality OpenType or TrueType fonts.
Best,
Jeff
--
http://biostat.mc.vanderbilt.edu/JeffreyHorner
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
