Re: [R] Help!!! from R beginner
Hello Jim,
Thank you so much!!! It is a magic!!! Now I understand the use of
ifelse. Thank you again!!!
Tae-Jin
On Feb 5, 2011, at 12:49 AM, jim holtman wrote:
You should be able to use 'ifelse'
Os.chr4.gene.new$color -
ifelse(Os.chr4.gene.new$if_TE_related == TE_related, black,
orange)
On Fri, Feb 4, 2011 at 7:09 PM, Tae-Jin Lee tj...@ncsu.edu wrote:
Hello,
I'm trying to add a column to the following data frame. The new
column
will contain black when the 5th column(if_TE_related) is
TE_related, or orange when the 4th column is (space).
chromoMSU_locus end5 end3 if_TE_related
chr04 LOC_Os04g0100610322679TE_related
chr04 LOC_Os04g0100876363951TE_related
chr04 LOC_Os04g01010952110296 TE_related
chr04 LOC_Os04g0102017165 17437
chr04 LOC_Os04g0103029372 18440 TE_related
chr04 LOC_Os04g0104030637 37300 TE_related
...
So, after a data manipulation, it should look like the following...
chromoMSU_locus end5 end3 if_TE_related
color
chr04 LOC_Os04g0100610322679TE_related
black
chr04 LOC_Os04g0100876363951TE_related
black
chr04 LOC_Os04g01010952110296 TE_related
black
chr04 LOC_Os04g0102017165 17437 orange
chr04 LOC_Os04g0103029372 18440 TE_related
black
chr04 LOC_Os04g0104030637 37300 TE_related
black
...
I have worked on the following code to do this job using function and
loop, but it is not working. If someone help me, I would really
appreciate!!!
The original data frame is Os.chr4.gene.new.
Gene - Os.chr4.gene.new[, c(if_TE_related)]
Genecolor - function(Gene) {
lg -length(Gene)
for(i in 1:lg) {
if (Gene == TE_related) {D1 - (Gene == black)}
if (Gene == ) {D1 - (Gene == orange)}
}
Gene.color - cbind(Gene, D1)
write.table(Gene.color, file=Gene_color1.txt, sep=\t,
row.names=F)
}
Genecolor(Gene)
Tae-Jin
Researcher in NC State University
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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
Data Munger Guru
What is the problem that you are trying to solve?
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Create factorial design
I have got data with one column indicating the area where the data
was recorded:
R: n - 43
R: df - data.frame(area=sample(1:7, n, repl=T), dat=rnorm(n))
In each of the 7 different areas I want to implement one of 7
specific strategies. The assignment should be random. Therefore, I
pair 7 areas with 7 strategies randomly by
R: ass - as.data.frame(cbind(area=sample(1:7, 7),
strategy=sample(1:7, 7)))
Now I want to create a new variable indicating, which case in the
original data should be assigned to which strategy. I thought
about
R: x - numeric(n)
R: for(i in 1:7){
x[df[, area]==i] - ass[ ass[, area]==i , strategy]
}
and then binding the new variable to the data frame
R: str(df2 - as.data.frame(cbind(df, strategy=x)))
which works fine. My question is whether there is a more elegant
way?
Thanks, *S*
--
Sascha Vieweg, saschav...@gmail.com
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Failing to install {rggobi} on win-7 R 2.12.0
Dear Prof Brian Ripley and others, After (finally) checking again, I found that ggobi doesn't work with the newer GTK (probably needed for R 2.12.0). Here is the results of my experimentations: *What works:* Installing ggobi and the GTK (version 2.12.9) provided on: http://www.ggobi.org/downloads/ Will *work* with R 2.11.1 + rggobi (version 2.1.14) (notice that the newest version 2.1.16, is not available for windows users. And the newer version 2.1.16-3 is not available on CRAN: http://cran.r-project.org/web/packages/rggobi/index.html) *What doesn't work:* Trying this with R 2.12.0 with rggobi (version 2.1.16-3), will *crash* with the error: the procedure entry point g_malloc_n could not be located in the dynamic link library libglic-2.0-0.dll When then trying to install gtk+ from the dialog box offered, then it downloads GTK (version 2.22.0 instead of 2.12.9) from here: http://downloads.sourceforge.net/gtk-win/gtk2-runtime-2.22.0-2010-10-21-ash.exe?download ' After doing this, ggobi (not rggobi) itself, won't run. It will crash with the error: the procedure entry point g_assertion_message_error could not be located in the dynamic link library libglib-2.0-0.dll Trying then to insert the dll's from the packages Prof Brian Ripley suggested (libxml2.dll and iconv.dll) won't fix the problem. Also using the libglib-2.0-0.dll from the old GTK package won't help. Also, reinstalling ggobi won't work. Any suggestions on how to make ggobi work with the newer version of GTK ? Thanks, Tal Contact Details:--- Contact me: tal.gal...@gmail.com | 972-52-7275845 Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) | www.r-statistics.com (English) -- On Wed, Jan 26, 2011 at 4:50 PM, Tal Galili tal.gal...@gmail.com wrote: I checked it using: Sys.getenv(PATH) And the output includes the PATH to the GTK2 installation (it's the last item in the following list): C:\\Windows\\system32;C:\\Windows;C:\\Windows\\System32\\Wbem;C:\\Windows\\System32\\WindowsPowerShell\\v1.0\\;C:\\Program Files (x86)\\Common Files\\Ulead Systems\\MPEG;C:\\Program Files\\TortoiseGit\\bin;C:\\Program Files (x86)\\QuickTime\\QTSystem\\;C:\\Program Files (x86)\\ggobi;C:\\Program Files (x86)\\GTK2-Runtime\\bin What else might I try? (Thanks) Contact Details:--- Contact me: tal.gal...@gmail.com | 972-52-7275845 Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) | www.r-statistics.com (English) -- On Wed, Jan 26, 2011 at 4:23 PM, Prof Brian Ripley rip...@stats.ox.ac.ukwrote: Your GTK+ installation is not being found: check your PATH. On Wed, 26 Jan 2011, Tal Galili wrote: Hello Prof Brian Ripley, Yihui and Tom, Thank you for your suggestions. It seemed to have made some differences in the error massages - but rggobi still fails to load. Steps taken: 1) I removed the old GTK (through the uninstall interface) 2) I ran library(RGtk2) which downloaded the new GTK-runtime version 2.22.0-2010-10-21 (instead of the one I got from ggobi, which was 2.12.9-2). 3) I downloaded both ftp://ftp.zlatkovic.com/libxml/libxml2-2.7.7.win32.zip and ftp://ftp.zlatkovic.com/libxml/iconv-1.9.2.win32.zip Unzipped them, and moved their dll's (from their bin directory), into - C:\Program Files (x86)\GTK2-Runtime\bin 4) I then tried starting rggobi: library(rggobi) and got the following error massages: Error 1: the program can't start because libgdk-win32-2.0-0.dll is missing from your computer. Try reinstalling the program to fix this problem. It then tried to reinstall GTK, and after I refused to, it sent the second Error massage: the program can't start because libfreetype-6.dll is missing from your computer. Try reinstalling the program to fix this problem. Any suggestions what else I should try? Many thanks for helping, Tal Contact Details:--- Contact me: tal.gal...@gmail.com | 972-52-7275845 Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) | www.r-statistics.com (English) --- --- On Wed, Jan 26, 2011 at 9:17 AM, Prof Brian Ripley rip...@stats.ox.ac.uk wrote: On Tue, 25 Jan 2011, Tom La Bone wrote: I recall that my problem on Windows was related to having a number of stray versions of GTK+ installed. I went back and deleted all versions and reinstalled the latest GTK+ and that seemed to fix things. However, when I went
Re: [R] Create factorial design
On Sat, Feb 05, 2011 at 11:01:33AM +0100, Sascha Vieweg wrote:
I have got data with one column indicating the area where the data
was recorded:
R: n - 43
R: df - data.frame(area=sample(1:7, n, repl=T), dat=rnorm(n))
In each of the 7 different areas I want to implement one of 7
specific strategies. The assignment should be random. Therefore, I
pair 7 areas with 7 strategies randomly by
R: ass - as.data.frame(cbind(area=sample(1:7, 7),
strategy=sample(1:7, 7)))
Now I want to create a new variable indicating, which case in the
original data should be assigned to which strategy. I thought
about
R: x - numeric(n)
R: for(i in 1:7){
x[df[, area]==i] - ass[ ass[, area]==i , strategy]
}
and then binding the new variable to the data frame
R: str(df2 - as.data.frame(cbind(df, strategy=x)))
which works fine. My question is whether there is a more elegant
way?
Hello.
If the table ass is sorted according to area, then its second
column may be used as a function mapping area to strategy. This
leads to the following
ass2 - ass[order(ass[, area]), strategy]
y - ass2[df[, area]]
identical(x, y + 0)
[1] TRUE
This also suggests that the same distribution on the random assignments
is obtained, if area is created already sorted and only the second
column of ass is random
ass - as.data.frame(cbind(area=1:7, strategy=sample(1:7, 7)))
Whether creating only this table is sufficient, depends on the application.
Hope this helps.
Petr Savicky.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Failing to install {rggobi} on win-7 R 2.12.0
Hi Tal, Thanks for working through this. GGobi needs to be rebuilt for the new version of GTK+. I'm probably the person to do that, but my time is short these days. I'll try to get to it soon. The new binary will just include the necessary DLLs, so that this GTK+ installation step is unnecessary. Michael On Sat, Feb 5, 2011 at 2:57 AM, Tal Galili tal.gal...@gmail.com wrote: Dear Prof Brian Ripley and others, After (finally) checking again, I found that ggobi doesn't work with the newer GTK (probably needed for R 2.12.0). Here is the results of my experimentations: *What works:* Installing ggobi and the GTK (version 2.12.9) provided on: http://www.ggobi.org/downloads/ Will *work* with R 2.11.1 + rggobi (version 2.1.14) (notice that the newest version 2.1.16, is not available for windows users. And the newer version 2.1.16-3 is not available on CRAN: http://cran.r-project.org/web/packages/rggobi/index.html) *What doesn't work:* Trying this with R 2.12.0 with rggobi (version 2.1.16-3), will *crash* with the error: the procedure entry point g_malloc_n could not be located in the dynamic link library libglic-2.0-0.dll When then trying to install gtk+ from the dialog box offered, then it downloads GTK (version 2.22.0 instead of 2.12.9) from here: http://downloads.sourceforge.net/gtk-win/gtk2-runtime-2.22.0-2010-10-21-ash.exe?download ' After doing this, ggobi (not rggobi) itself, won't run. It will crash with the error: the procedure entry point g_assertion_message_error could not be located in the dynamic link library libglib-2.0-0.dll Trying then to insert the dll's from the packages Prof Brian Ripley suggested (libxml2.dll and iconv.dll) won't fix the problem. Also using the libglib-2.0-0.dll from the old GTK package won't help. Also, reinstalling ggobi won't work. Any suggestions on how to make ggobi work with the newer version of GTK ? Thanks, Tal Contact Details:--- Contact me: tal.gal...@gmail.com | 972-52-7275845 Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) | www.r-statistics.com (English) -- On Wed, Jan 26, 2011 at 4:50 PM, Tal Galili tal.gal...@gmail.com wrote: I checked it using: Sys.getenv(PATH) And the output includes the PATH to the GTK2 installation (it's the last item in the following list): C:\\Windows\\system32;C:\\Windows;C:\\Windows\\System32\\Wbem;C:\\Windows\\System32\\WindowsPowerShell\\v1.0\\;C:\\Program Files (x86)\\Common Files\\Ulead Systems\\MPEG;C:\\Program Files\\TortoiseGit\\bin;C:\\Program Files (x86)\\QuickTime\\QTSystem\\;C:\\Program Files (x86)\\ggobi;C:\\Program Files (x86)\\GTK2-Runtime\\bin What else might I try? (Thanks) Contact Details:--- Contact me: tal.gal...@gmail.com | 972-52-7275845 Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) | www.r-statistics.com (English) -- On Wed, Jan 26, 2011 at 4:23 PM, Prof Brian Ripley rip...@stats.ox.ac.ukwrote: Your GTK+ installation is not being found: check your PATH. On Wed, 26 Jan 2011, Tal Galili wrote: Hello Prof Brian Ripley, Yihui and Tom, Thank you for your suggestions. It seemed to have made some differences in the error massages - but rggobi still fails to load. Steps taken: 1) I removed the old GTK (through the uninstall interface) 2) I ran library(RGtk2) which downloaded the new GTK-runtime version 2.22.0-2010-10-21 (instead of the one I got from ggobi, which was 2.12.9-2). 3) I downloaded both ftp://ftp.zlatkovic.com/libxml/libxml2-2.7.7.win32.zip and ftp://ftp.zlatkovic.com/libxml/iconv-1.9.2.win32.zip Unzipped them, and moved their dll's (from their bin directory), into - C:\Program Files (x86)\GTK2-Runtime\bin 4) I then tried starting rggobi: library(rggobi) and got the following error massages: Error 1: the program can't start because libgdk-win32-2.0-0.dll is missing from your computer. Try reinstalling the program to fix this problem. It then tried to reinstall GTK, and after I refused to, it sent the second Error massage: the program can't start because libfreetype-6.dll is missing from your computer. Try reinstalling the program to fix this problem. Any suggestions what else I should try? Many thanks for helping, Tal Contact Details:--- Contact me: tal.gal...@gmail.com | 972-52-7275845 Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) | www.r-statistics.com (English)
Re: [R] Quadratic regression: estimating the maximizing value
From: greg.s...@imail.org
To: ghe...@blm.gov; r-help@r-project.org
Date: Fri, 4 Feb 2011 17:49:51 -0700
Subject: Re: [R] Quadratic regression: estimating the maximizing value
No, your approach is not correct. For one you have not taken the covariance
between B and C into account, another thing is that the ratio of 2 normal
random variables is not necessarily normal, in some cases it can even follow
a Cauchy distribution. Also note that with only 1 degree of freedom the
Central Limit Theorem is not justified for using normal theory for non normal
distributions. So the normal based confidence intervals on the coefficients
are only reasonable if you are certain that the true error distribution is
normal or nearly normal.
To make this df issue more clear, take some code like this,
for ( i in 1:4)
+ {
+ Ldat - data.frame(X=c(3,7,14,24), Y=c(1,5,8,0))
+ Ldat=Ldat[which( 1:4 != i),]
+ LM-lm(formula = Y ~ X + I(X^2), data = Ldat)
+ print ( DZ(LM$coefficients[2],LM$coefficients[3]))
+ }
X
13.46512
X
13.1519
X
13.11364
X
14.625
and try to compute confidence intervals with 3 data points fitted to 3
parameters.
What exactly does this mean? The result is exact if you only have 3 points or
that there is really no way to tell until you have a lot more data than
parameters?
With your original data, look at residuals,
plot(LM)
and note the residuals are all zero with any 3 data points.
There is only so much you can do with limited data esp as number of points
approaches number of fitting parameters ( this is a common joke, always have
more parameters than data points ). Whatever reasons you have for expecting
the form to fit and things like measurement noise may be important to include
in your analysis to get anything meaningful out of it.
Some possibilities that may work for you:
Use the nls function with the curve parameterized to use the vertex point as
one of the parameters (still need normality).
Do bootstrapping (with only 4 points you are not going to get much with
non-parametric bootstrap, but parametric with some assumptions may work).
Use a Bayesian approach (this may be best if you can find some meaningful
prior information).
--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
project.org] On Behalf Of ghe...@blm.gov
Sent: Friday, February 04, 2011 2:34 PM
To: r-help@r-project.org
Subject: [R] Quadratic regression: estimating the maximizing value
A bioligist colleague sent me the following data.
x Y
3 1
7 5
14 8
24 0
(Yes, only four data points.) I don't know much about the
application, but apparently there are good empirical
reasons to use a quadratic model.
The goal is to find the X value which maximizes the
response Y, and to find a confidence interval for this X
value.
Finding the maximizing X value is pretty straightforward:
Ldat - data.frame(X=c(3,7,14,24), Y=c(1,5,8,0))
(LM-lm(formula = Y ~ X + I(X^2), data = Ldat))
Call:
lm(formula = Y ~ X + I(X^2), data = Ldat)
Coefficients:
(Intercept) X I(X^2)
-3.86978 1.77762 -0.06729
DZ-function(B,C) { (-B)/(2*C) } # Solve d/dx(A + Bx + Cx^2) = 0
DZ(LM$coefficients[2],LM$coefficients[3])
X
13.20961
To find a confidence interval, I used confint().
Default confidence level of 95% was not useful; used 80% instead,
and then computed DZ with the extreme X and I(X^2) coefficients:
(CI80-confint(LM,level=0.8))
10 % 90 %
(Intercept) -5.6147948 -2.12476138
X 1.4476460 2.10759306
I(X^2) -0.0790312 -0.05553898
DZ(CI80[2,1],CI80[3,1])
[1] 9.1587
DZ(CI80[2,2],CI80[3,2])
[1] 18.97400
Conclusion: the 80% confidence level for the maximizing X value is
included in the range (9.158, 18.974)
#
Questions:
1) Is this the right procedure, or totally off base?
2) The coefficient of the Intercept is irrelevant to calculating
the maximizing value of X. Is there a way to reduce the size of
the confidence interval by doing a computation that leaves out this
parameter?
3) I believe that confint() indicates the axes of an ellipsoid,
rather than the corners of a box, in parameter space;
so that the above procedure is (slightly) too conservative.
4) Where are the calculations for confint() documented ?
Thanks,
George Heine
(Embedded image moved to file: pic23206.gif)Please consider the
environment
before printing this page
=
=
George Heine, PhD
Mathematical Analyst
National Operations
Center
Bureau of Land Management
voice (303) 236-0099
fax (303) 236-1974
cell (303) 905-5296
=
=
__
R-help@r-project.org mailing list
Re: [R] stochastic models for population growth
Hello, It's not easy to express clearly what I have in mind. I've been working with a couple of titles -Stochastic Models in Biology, for one, but both date back a decade or more. I'll try to illustrate the model a little more. A reasonably comparable situation is Bolker's analysis of how many flowers are created on a given plant among a patch of plants, then comparing patches of plants at different locations. Then, I believe it's possible to determine the degree to which number of flowers appears to be dependent on other variables, like plant height, or root depth. Finally, I would like to add to these analyses, one which, having determined a stochastically(if I'm using this word correctly) distributed parameter for the average number of flowers for each patch (patch 1(N(mu, sigma)), patch 2(N(mu, sigma)),...) would constitute a growth model for each patch over time, I think maybe something like: df(patch 1) =r1*f0+ exp^(r1*t) dt df(patch 2) =r2*f0+ exp^(r2*t) dt It would be nice to be able to add the other effects here: plant height, and others, as further parameters in each equation. I think it's likely these models have been done in biology, though I don't seem to find exactly how to do it in Bolker's book, and as a non-biologist, I 'm not sure where to look to find it. Maybe one can see that I'd like to be able to transfer this kind of model to use with other kinds of data from other fields! -regards, s --- On Sun, 1/23/11, Michael D mike...@gmail.com wrote: From: Michael D mike...@gmail.com Subject: Re: [R] stochastic models for population growth To: shv...@yahoo.com, r-help@r-project.org Received: Sunday, January 23, 2011, 3:41 AM It's not too clear to me what you plan to do. You want to model population growth using a normal distribution? You should consider the classic differential equation of population growth and look at variants with species interaction. For modeling a single species you want to have dP/dt = r*P (r is rate of increase accounting for birth and death) which is P(t)=c1 e^(rt) (c1 you solve from your initial conditions) Then you can look to add a N(0,s) random component for If you want to do some competition for resources/predator-prey modeling there's already lots done on those. As far as an R library you can download and just enter a few parameters... I'm sure there's something in: http://cran.r-project.org/web/views/Environmetrics.html Michael D On Thu, Jan 20, 2011 at 23:57:18 -0800 (PST), Vassily Shvets shv736_at_yahoo.comwrote: Hello, Having measured two populations' characteristics at one particular time[with great precision] with R, I would like to extend this to measuring the same populations starting at t1, and then again at t2, and try to develop a growth model (something like dpop1/dt=r*pop^(...),dpop2/dt=r*pop^(...)). I think the idea is to create a model that will predict the growth of a population(N(mu, sigma)) within a margin of error. This kind of modeling isn't well known or publicized in terms of R, am I right? regards, s [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] spline interpolation
Hello R-help I have the following data for a standard curve concentration(nM),fluorescence 0,48.34 2,58.69 5,70.83 10,94.73 20,190.8 50,436.0 100, 957.9 (1)Is there function in R to plot a spline. (2)How can I interpolation,say 1000 point from 0nM-100nM and store this as a data frame of concentration,fluorescence (3)How can I modify the code below so that instead of retrieving a concentration with the exact value of fluorescence, it gives me concentration for the value that is closest to that fluorescence. subset(df,fluorescence==200.3456,select=concentration) [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] clustering fuzzy
After ordering the table of membership degrees , i must get the difference
between the first and second coloumns , between the first and second largest
membership degree of object i. This for K=2,K=3,to K.max=6.
This difference is multiplyed by the Crisp silhouette index vector (si). Too
it dependending on K=2,...,K.max=6; the result divided by the sum of these
differences
I need a final vector composed of the indexes for each clustering
(K=2,...,K.max=6).
There is a method, i think that is classe.memb, but i can't to solve problem
because trasformation of the membership degrees matrix( (ris$membership) and
of list object (ris$silinfo), does not permitme to use classe.memb
propertyes.
.
Σí(uί1-uí2)sí/Σí(uí1-uí2)
head(t(A.sort)) membership degrees table ordering by max to min value
[,1] [,2] [,3] [,4]
1 0.66 0.30 0.04 0.01
2 0.89 0.09 0.02 0.00
3 0.92 0.06 0.01 0.01
4 0.71 0.21 0.07 0.01
5 0.85 0.10 0.04 0.01
6 0.91 0.04 0.02 0.02
head(t(A.sort))
[,1] [,2] [,3] [,4]
1 0.66 0.30 0.04 0.01
2 0.89 0.09 0.02 0.00
3 0.92 0.06 0.01 0.01
4 0.71 0.21 0.07 0.01
5 0.85 0.10 0.04 0.01
6 0.91 0.04 0.02 0.02
H.Asort=head(t(A.sort))
H.Asort[,1]-H.Asort[,2]
123456
0.36 0.80 0.86 0.50 0.75 0.87
H.Asort=t(H.Asort[,1]-H.Asort[,2])
This is the differences vector by multiplying trasformed table ris$silinfo.
ris$silinfo
$widths
cluster neighbor sil_width
72 13 0.43820207
54 13 0.43427773
29 16 0.41729079
62 16 0.40550562
64 16 0.32686757
32 13 0.30544722
45 13 0.30428723
79 13 0.30192624
12 13 0.30034472
60 16 0.29642495
41 13 0.29282778
113 0.28000788
85 13 0.24709237
74 13 0.239
P=ris$silinfo
P=P[1]
P=as.data.frame(P)
V4=rownames(P)
mode(V4)=numeric
P[,4]=V4
P[order(P$V4),]
widths.cluster widths.neighbor widths.sil_width V4
1 1 3 0.28000788 1
2 2 4 0.07614849 2
3 2 3 -0.11676440 3
4 2 4 0.15436648 4
5 2 3 0.14693927 5
6 3 1 0.57083836 6
7 4 5 0.36391826 7
8 5 4 0.63491118 8
9 4 2 0.54458733 9
10 5 4 0.51059626 10
11 2 5 0.03908952 11
12 1 3 0.30034472 12
13 1 3 -0.04928562 13
14 4 3 0.20337180 14
15 3 4 0.46164324 15
18 5 4 0.52066782 18
20 4 3 0.45517287 20
21 3 4 0.39405507 21
22 4 5 0.05574547 22
23 6 1 -0.06750403 23
P= P[order(P$V4),]
P=P[,3]
This is trasformed vector ris$silinfo =P.
I can't to use this vector object in the classe.memb.
K=2
K.max=6
while (K=K.max)
{
ris=fanny(frj,K,memb.exp=m,metric=SqEuclidean,stand=TRUE,maxit=1000,tol=1e-6)
ris$centroid=matrix(0,nrow=K,ncol=J)
for (k in 1:K)
{
ris$centroid[k,]=(t(ris$membership[,k]^m)%*%as.matrix(frj))/sum(ris$membership[,k]^m)
}
rownames(ris$centroid)=1:K
colnames(ris$centroid)=colnames(frj)
print(K)
print(round(ris$centroid,2))
print(classe.memb(ris$membership)$table.U)
print(ris$silinfo$avg.width)
K=K+1
}
this should be scheme clearly are determined centroid based on classe.memb.
classe.memb=function(U)
{
info.U=cbind(max.col(U),apply(U,1,max))
i=1
while (i = nrow(U))
{
if (apply(U,1,max)[i]0.5) info.U[i,1]=0
i=i+1
}
K=ncol(U)
table.U=matrix(0,nrow=K,ncol=4)
cl=1
while (cl = K)
{
table.U[cl,1] = length(which(info.U[info.U[,1]==cl,2]=.90))
table.U[cl,2] = length(which(info.U[info.U[,1]==cl,2]=.70)) -
table.U[cl,1]
table.U[cl,3] = length(which(info.U[info.U[,1]==cl,2]=.50)) -
table.U[cl,1] - table.U[cl,2]
table.U[cl,4] = sum(table.U[cl,])
cl = cl+1
}
rownames(table.U) = c(1:K)
colnames(table.U) = c(Alto, Medio, Basso, Totale)
out=list()
out$info.U=round(info.U,2)
out$table.U=table.U
return(out)
}
--
View this message in context:
http://r.789695.n4.nabble.com/clustering-fuzzy-tp3229853p3261837.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] [OT] BLUE vs. BLUP
On Fri, Feb 4, 2011 at 8:44 PM, Laura Smith smithlaura...@gmail.com wrote: Hi! Does anyone have a numeric example for calculating BLUE and BLUP, please? You will need to be more specific. BLUE is an acronym for Best Linear Unbiased Estimator and BLUP for Best Linear Unbiased Predictor. So the phrase calculating BLUE is incomplete. You need to say calculating the BLUE of __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] different results in MASS's mca and SAS's corresp
On Feb 4, 2011, at 7:06 PM, Gong-Yi Liao wrote: Dear list: I have tried MASS's mca function and SAS's PROC corresp on the farms data (included in MASS, also used as mca's example), the results are different: R: mca(farms)$rs: 1 2 1 0.059296637 0.0455871427 2 0.043077902 -0.0354728795 3 0.059834286 0.0730485572 4 0.059834286 0.0730485572 5 0.012900181 -0.0503121890 6 0.038846577 -0.0340961617 7 0.005886752 -0.0438516465 8 -0.015108789 -0.0247221783 9 0.007505626 -0.0646608108 10 0.006631230 -0.0362117073 11 0.013309217 -0.0680733730 12 0.056549933 0.0010773359 13 0.015681958 0.0334320046 14 -0.065598990 0.0151619769 15 -0.046868229 0.0357782553 16 -0.003048803 0.0128157261 17 -0.051281437 0.0278941743 18 -0.051819085 0.0004327598 19 -0.072814626 0.0195622280 20 -0.072814626 0.0195622280 And in SAS's corresp output: Row Coordinates Dim1 Dim2 1 1.0607-0.8155 2 0.7706 0.6346 3 1.0703-1.3067 4 1.0703-1.3067 5 0.2308 0.9000 6 0.6949 0.6099 7 0.1053 0.7844 8 -0.2703 0.4422 9 0.1343 1.1567 10 0.1186 0.6478 11 0.2381 1.2177 12 1.0116-0.0193 13 0.2805-0.5980 14 -1.1735-0.2712 15 -0.8384-0.6400 16 -0.0545-0.2293 17 -0.9174-0.4990 18 -0.9270-0.0077 19 -1.3025-0.3499 20 -1.3025-0.3499 Is MASS's mca developed with different definition to SAS's corresp ? No, it's just that the values can only be defined up to a scaling factor (the same situation as with eigenvector decompostion). Take a look at the two dimensions, when each is put on the same scale: cbind(scale(rmca$D1),scale(smca$Dim1) ) [,1][,2] [1,] 1.2824421 1.28242560 [2,] 0.9316703 0.93168561 [3,] 1.2940701 1.29403231 [4,] 1.2940701 1.29403231 [5,] 0.2789996 0.27905048 [6,] 0.8401570 0.84016193 [7,] 0.1273161 0.12731705 [8,] -0.3267664 -0.32679513 [9,] 0.1623284 0.16237896 [10,] 0.1434174 0.14339716 [11,] 0.2878460 0.28787641 [12,] 1.2230376 1.22306216 [13,] 0.3391626 0.33913934 [14,] -1.4187467 -1.41879225 [15,] -1.0136458 -1.01364584 [16,] -0.0659382 -0.06588616 [17,] -1.1090928 -1.10915932 [18,] -1.1207208 -1.12076602 [19,] -1.5748033 -1.57475730 [20,] -1.5748033 -1.57475730 cbind(scale(rmca$D2),scale(smca$Dim2) ) [,1][,2] [1,] 1.06673426 -1.06677626 [2,] -0.83006158 0.83012474 [3,] 1.70932841 -1.70932351 [4,] 1.70932841 -1.70932351 [5,] -1.17729983 1.17729909 [6,] -0.79784653 0.79781424 [7,] -1.02612383 1.02608072 [8,] -0.57849632 0.57844296 [9,] -1.51305605 1.51309282 [10,] -0.84735007 0.84739189 [11,] -1.59290964 1.59288798 [12,] 0.02520954 -0.02525321 [13,] 0.78230533 -0.78226073 [14,] 0.35478864 -0.35476797 [15,] 0.83720734 -0.83720166 [16,] 0.29988662 -0.29995785 [17,] 0.65272069 -0.65275711 [18,] 0.01012653 -0.01007904 [19,] 0.45775404 -0.45771681 [20,] 0.45775404 -0.45771681 -- David. Thank you for any comments! -- Gong-Yi Liao Department of Statistics University of Connecticut -- David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] spline interpolation
On Feb 5, 2011, at 9:29 AM, Asan Ramzan wrote: Hello R-help I have the following data for a standard curve concentration(nM),fluorescence 0,48.34 2,58.69 5,70.83 10,94.73 20,190.8 50,436.0 100, 957.9 (1)Is there function in R to plot a spline. ?spline (2)How can I interpolation,say 1000 point from 0nM-100nM and store this as a data frame of concentration,fluorescence ?approxfun (3)How can I modify the code below so that instead of retrieving a concentration with the exact value of fluorescence, it gives me concentration for the value that is closest to that fluorescence. ?which.min #(although I would be using [ rather than subset for that purpose.) subset(df,fluorescence==200.3456,select=concentration) [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] plot two pearson correlation matrix on one heatmap
Hi all, I am trying to plot two pearson correlation coefficient matrices on one heatmap using heatmap.2. I combined two matrices as following A = [0.8 0.9; 0.75 0.95] B = [0.82 0.87 0.94 0.74] combined = [0.8 0.9 0 0; 0.75 0.95 0 0 0 0 0.82 0.87 0 0 0.94 0.74] heatmap.2(combined,distfun=function(c) as.dist(1-c),hclustfun=function(c) hclust(c,method=average), race=none, col=greenred(75),trace=none,dendrogram =column) The resultant heatmap having a scale from 0 to 1. I want to ignore the zero elements in the combined matrix and display them in black. Also, I want the scale to be greater than 0.7. I looked at the heatmap.2 webpage, but did not have a clue. Could anybody help me? Thank you in advance. Wendy -- View this message in context: http://r.789695.n4.nabble.com/plot-two-pearson-correlation-matrix-on-one-heatmap-tp3261882p3261882.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] spline interpolation
Also ?splinefun (like approxfun but using splines instead of lines). -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare greg.s...@imail.org 801.408.8111 -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- project.org] On Behalf Of David Winsemius Sent: Saturday, February 05, 2011 8:24 AM To: Asan Ramzan Cc: r-help@r-project.org Subject: Re: [R] spline interpolation On Feb 5, 2011, at 9:29 AM, Asan Ramzan wrote: Hello R-help I have the following data for a standard curve concentration(nM),fluorescence 0,48.34 2,58.69 5,70.83 10,94.73 20,190.8 50,436.0 100, 957.9 (1)Is there function in R to plot a spline. ?spline (2)How can I interpolation,say 1000 point from 0nM-100nM and store this as a data frame of concentration,fluorescence ?approxfun (3)How can I modify the code below so that instead of retrieving a concentration with the exact value of fluorescence, it gives me concentration for the value that is closest to that fluorescence. ?which.min #(although I would be using [ rather than subset for that purpose.) subset(df,fluorescence==200.3456,select=concentration) [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting- guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] spline interpolation
On Feb 5, 2011, at 10:37 AM, Mike Marchywka wrote: From: dwinsem...@comcast.net To: asanram...@yahoo.com Date: Sat, 5 Feb 2011 10:24:04 -0500 CC: r-help@r-project.org Subject: Re: [R] spline interpolation On Feb 5, 2011, at 9:29 AM, Asan Ramzan wrote: Hello R-help I have the following data for a standard curve concentration(nM),fluorescence 0,48.34 2,58.69 5,70.83 10,94.73 20,190.8 50,436.0 100, 957.9 (1)Is there function in R to plot a spline. ?spline (2)How can I interpolation,say 1000 point from 0nM-100nM and store this as a data frame of concentration,fluorescence ?approxfun I was waiting for some other code to finish and thought I would give this a try. I'm still learning R, maybe you can comment on this approach # made data file f1-read.table(d1.txt,sep=,) ss-smooth.spline(f1$V1,f1$V2) # it is alwys good to look first. plot(ss) lines(ss) str(ss) ssi-smooth.spline(f1$V2,f1$V1) predict(ssi,200.3456) predict(ssi,1:10) predict(ssi,1:100)$y Works for me. smooth.spline() is better than spline() for noisy data. spline() hits every point. And there are further good examples on ? spline and ?smooth.spline. The splines package alos provides the bs() and ns() functions for B-splines and natural splines (=cubic splines) as terms in model formulae. (There's no protection against extrapolation outside the domain) f1 - data.frame(x=1:100, y=log(1:100)+rnorm(100)) plot(f1) ss-smooth.spline(f1$x,f1$y) plot(ss, type=l) points(f1) predict(ss, 200.3456) ## $x [1] 200.3456 $y [1] 8.185104 etc (3)How can I modify the code below so that instead of retrieving a concentration with the exact value of fluorescence, it gives me concentration for the value that is closest to that fluorescence. ?which.min #(although I would be using [ rather than subset for that purpose.) subset(df,fluorescence==200.3456,select=concentration) David Winsemius, MD West Hartford, CT David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] spline interpolation
From: dwinsem...@comcast.net To: asanram...@yahoo.com Date: Sat, 5 Feb 2011 10:24:04 -0500 CC: r-help@r-project.org Subject: Re: [R] spline interpolation On Feb 5, 2011, at 9:29 AM, Asan Ramzan wrote: Hello R-help I have the following data for a standard curve concentration(nM),fluorescence 0,48.34 2,58.69 5,70.83 10,94.73 20,190.8 50,436.0 100, 957.9 (1)Is there function in R to plot a spline. ?spline (2)How can I interpolation,say 1000 point from 0nM-100nM and store this as a data frame of concentration,fluorescence ?approxfun I was waiting for some other code to finish and thought I would give this a try. I'm still learning R, maybe you can comment on this approach # made data file f1-read.table(d1.txt,sep=,) ss-smooth.spline(f1$V1,f1$V2) # it is alwys good to look first. plot(ss) lines(ss) str(ss) ssi-smooth.spline(f1$V2,f1$V1) predict(ssi,200.3456) predict(ssi,1:10) predict(ssi,1:100)$y etc (3)How can I modify the code below so that instead of retrieving a concentration with the exact value of fluorescence, it gives me concentration for the value that is closest to that fluorescence. ?which.min #(although I would be using [ rather than subset for that purpose.) subset(df,fluorescence==200.3456,select=concentration) David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Where does R look for library packages - there is no package called 'BRugs'
On 05.02.2011 00:32, conor1725 wrote: I created a file called .Renviron and set R_LIBS_USER to the same path I had set R_LIBS to. I put this file in C:/Users/myusername/Documents/mysubdirectory. I also commented out the line I had put in the Renviron.site file. Now I get the same error as I originally got. When I put .Renvrion in C:/Users/myusername/Documents instead, it does work, but the original reason I posted my question was my desire to not have anything R (or BRugs) related directly in my C:/Users/myusername/Documents folder. If there's a way to do this the correct way without putting any files or folders in my C:/Users/myusername/Documents folder as opposed to my C:/Users/myusername/Documents/mysubdirectory folder or my C:/Program Files/R/R-2.12.0, I'd still be interested in knowing how. Then just define an environment variable in your OS settings. Uwe Ligges Prof Brian Ripley wrote: For the record, this is not a good answer to the question (and the real answer is 'explained fully' on the help page for .libPaths, and also in the rw-FAQ, Q4.2). If you move your personal library directory, you need to set the environment variable R_LIBS_USER . And personal (not site) environment variables should be set in ~/.Renviron, not in Renviron.site. If you want a site library, that is explained fully in the references in my first paragraph. On Thu, 3 Feb 2011, conor1725 wrote: David Winsemius wrote: On Jan 30, 2011, at 5:42 PM, conor1725 wrote: I am have installed R on Windows 7 machine. R got installed in the directory C:\Program Files\R\R-2.12.0. Then I installed the package BRugs using the install.packages command. I did not get an option during the installation as to what directory I wanted BRugs installed in. It ended up in C:\myusername\Documents as C:\myusername\Documents\R\win-library \2.12\BRugs. When I open and R and run the command library(BRugs) it works fine. However, I would like to move the BRugs stuff out of my Documents folder, either to a subdirectory under Documents of my choosing or to somewhere in C:\Program Files\R\R-2.12.0. I tried just moving the whole folder and subfolders R\win-library\2.12\... from Documents to a subfolder. I don't think this matters, but just to complete, there is also a folder called Coda with the BRugs folder in the 2.12 folder and of course there are subfolders and files under each of the folders Coda and BRugs. Then when I try library(BRugs) I get the error message Error in library(BRugs) : there is no package called 'BRugs'. So it seems R only knows to look for the BRugs package in it's original installation location. How can I get it to work having BRugs in one of my desired locations? Then you need to understand where R will be able to find your packages, and how one might change that: ?libPaths -- David Winsemius, MD West Hartford, CT I'm posting a specific answer to my own question for the record since The help for libPaths (actually should be ?.libPaths) doesn't explain this fully. I created the file Renviron.site in the folder C:\Program Files\R\R-2.12.0\etc and put the text R_LIBS = C:/Users/myusername/Documents/.../R/win-library/2.12 in that file. The ... represents the path for the subdirectory that I moved R\win-library\2.12\... to. -- View this message in context: http://r.789695.n4.nabble.com/Where-does-R-look-for-library-packages-there-is-no-package-called-BRugs-tp3247748p3259520.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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, rip...@stats.ox.ac.uk 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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] spline interpolation
Hi, Just pure curiosity: may I ask why you want to do spline interplation on fluorescence intensity as function of concentration? Particularly as it looks quite typical for an unknown problem's calibration plot? Claudia On 02/05/2011 03:29 PM, Asan Ramzan wrote: Hello R-help I have the following data for a standard curve concentration(nM),fluorescence 0,48.34 2,58.69 5,70.83 10,94.73 20,190.8 50,436.0 100,�957.9 � (1)Is there function in R�to plot a spline. (2)How can I interpolation,say 1000 point from 0nM-100nM and store this as a data frame of concentration,fluorescence (3)How can I modify�the code below so that instead of retrieving a concentration with the exact value of fluorescence, it gives me concentration for the value that is closest to that fluorescence. � subset(df,fluorescence==200.3456,select=concentration) [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Claudia Beleites Dipartimento dei Materiali e delle Risorse Naturali Università degli Studi di Trieste Via Alfonso Valerio 6/a I-34127 Trieste phone: +39 0 40 5 58-37 68 email: cbelei...@units.it __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 maximum value of the legend on an image.plot and animation
On 01.02.2011 16:10, steve_fried...@nps.gov wrote:
Hello,
I'm doing the following:
library(ncdf)
library(fields)
library(animation)
saline- open.ncdf(salinity_1990.nc)
salt = get.var.ncdf(nc=saline, varid=Salinity)
# create an animation of the complete temporal domain in the ncdf file.
saveHTML({
for (i in 1:364) {
image.plot(salt[, , i])
}
}, img.name = salinity.img, imgdir = salinity_dir, htmlfile =
salinity.html, zmax = c(0, 80),
autobrowse = FALSE, title = TIME SALINITY PREDICTIONS,
description = c(This should plot 1 years daily salinity
predictions in Florida Bay)
)
Almost all of the data sets I work with are multi-temporal spatial forms.
Being able to view the time sequence is very important to us, The animation
procedure is very good. I appreciate its development.
Here are my questions.
1) I would like to find a procedure to set the maximum value of the legend
of an image.plot. Each time step has a different maximum value, but for
the animation to be valuable, I need to standardize these to a maximum such
that the scale is equal in each time step.
Specify the breaks manually using the argument breaks that can be
passed to ... in image.plot().
2) Secondly, when setting imgdir it only goes to a temp directory. If I
define a directory, the new directory path is added to the default temp
directory path. Is there a way to fix this so that the output directory is
truly defined in the saveHTML statements?
I don't see a way, so just ask the maintainer of the animation package
to allow for absolute paths or so.
Uwe Ligges
I'm working on a windows XP machine using R 2.12.1.
Thanks
Steve
Steve Friedman Ph. D.
Ecologist / Spatial Statistical Analyst
Everglades and Dry Tortugas National Park
950 N Krome Ave (3rd Floor)
Homestead, Florida 33034
steve_fried...@nps.gov
Office (305) 224 - 4282
Fax (305) 224 - 4147
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] using SNOW and clusterApplyLB to run jobs parallel
On 01.02.2011 23:07, kparamas wrote:
I have this function and want to run it parallel with different sets of data.
Using SNOW and clusterApplyLB.
Nice, but what is your question?
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Uwe Ligges
system.time(out- mclapply(cData, plotGraph)) #each cData contains 100X6000
doubles
system.time(out- mclapply(cData2, plotGraph))
system.time(out- mclapply(cData3, plotGraph))
system.time(out- mclapply(cData4, plotGraph))
system.time(out- mclapply(cData5, plotGraph))
system.time(out- mclapply(cData6, plotGraph))
plotGraph()- function(cData)
{
cl = unname(cor(cData))
result = cbind(as.vector(row(cl)),as.vector(col(cl)),as.vector(cl))
result = result[result[,1] != result[,2],]
corm = result
corm =corm[abs(corm[,3])= CORRELATION, ]
# remove low cor pairs
library(network); library(sna)
net- network(corm, directed = F)
# the network
cd- component.dist(net)
# component analysis
delete.vertices(net, which(cd$csize[cd$membership] == 1))
# delete genes not connected with others
plot(net)
}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Performance of graph and igraph package
On 05.02.2011 06:21, William Tu wrote: Dear R users, I'm using graph library to create a mesh-like network topology and implement a load balance routing algorithm. The current implementation uses graph, RBGL, and Rgraphviz libraries. I have a few attributes on every edge to represent the network loading and capacity, and I frequently update these values. However, I found it works so slow and right now I'm consider rewriting it using C or Python. Then I found another library, which is igraph. I made a performance comparison between graph and igraph as below. For both graph, I create two attributes for each edge and set/get the value of all edges and measure the total elapsed time. A. using graph package: 108 nodes, 832 edges: 17 sec 140 nodes, 4800 edges: 576 sec B. using igraph package: 100 nodes, 4950 edges: 4 sec 200 nodes, 19900 edges: 111 sec igraph is much much faster than graph! Is this reasonable? or I did something wrong? Calculating graphs efficiently is not a trivial task and you can optimize for different szenarios. Unless you doubt about the accuracy of results, I would not be worried about a factor of 5 in different implementations. Uwe Ligges ps. I'm using 2 core 3GHz CPU with 2G ram. Thanks in advance, William __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] clustering fuzzy
After ordering the table of membership degrees , i must get the difference
between the first and second coloumns , between the first and second largest
membership degree of object i. This for K=2,K=3,to K.max=6.
This difference is multiplyed by the Crisp silhouette index vector (si). Too
it dependending on K=2,...,K.max=6; the result divided by the sum of these
differences
I need a final vector composed of the indexes for each clustering
(K=2,...,K.max=6).
There is a method, i think that is classe.memb, but i can't to solve problem
because trasformation of the membership degrees matrix( (ris$membership) and
of list object (ris$silinfo), does not permitme to use classe.memb
propertyes.
.
Σí(uί1-uí2)sí/Σí(uí1-uí2)
head(t(A.sort)) membership degrees table ordering by max to min value
[,1] [,2] [,3] [,4]
1 0.66 0.30 0.04 0.01
2 0.89 0.09 0.02 0.00
3 0.92 0.06 0.01 0.01
4 0.71 0.21 0.07 0.01
5 0.85 0.10 0.04 0.01
6 0.91 0.04 0.02 0.02
head(t(A.sort))
[,1] [,2] [,3] [,4]
1 0.66 0.30 0.04 0.01
2 0.89 0.09 0.02 0.00
3 0.92 0.06 0.01 0.01
4 0.71 0.21 0.07 0.01
5 0.85 0.10 0.04 0.01
6 0.91 0.04 0.02 0.02
H.Asort=head(t(A.sort))
H.Asort[,1]-H.Asort[,2]
123456
0.36 0.80 0.86 0.50 0.75 0.87
H.Asort=t(H.Asort[,1]-H.Asort[,2])
This is the differences vector by multiplying trasformed table ris$silinfo.
ris$silinfo
$widths
cluster neighbor sil_width
72 13 0.43820207
54 13 0.43427773
29 16 0.41729079
62 16 0.40550562
64 16 0.32686757
32 13 0.30544722
45 13 0.30428723
79 13 0.30192624
12 13 0.30034472
60 16 0.29642495
41 13 0.29282778
113 0.28000788
85 13 0.24709237
74 13 0.239
P=ris$silinfo
P=P[1]
P=as.data.frame(P)
V4=rownames(P)
mode(V4)=numeric
P[,4]=V4
P[order(P$V4),]
widths.cluster widths.neighbor widths.sil_width V4
1 1 3 0.28000788 1
2 2 4 0.07614849 2
3 2 3 -0.11676440 3
4 2 4 0.15436648 4
5 2 3 0.14693927 5
6 3 1 0.57083836 6
7 4 5 0.36391826 7
8 5 4 0.63491118 8
9 4 2 0.54458733 9
10 5 4 0.51059626 10
11 2 5 0.03908952 11
12 1 3 0.30034472 12
13 1 3 -0.04928562 13
14 4 3 0.20337180 14
15 3 4 0.46164324 15
18 5 4 0.52066782 18
20 4 3 0.45517287 20
21 3 4 0.39405507 21
22 4 5 0.05574547 22
23 6 1 -0.06750403 23
P= P[order(P$V4),]
P=P[,3]
This is trasformed vector ris$silinfo =P.
I can't to use this vector object in the classe.memb.
K=2
K.max=6
while (K=K.max)
{
ris=fanny(frj,K,memb.exp=m,metric=SqEuclidean,stand=TRUE,maxit=1000,tol=1e-6)
ris$centroid=matrix(0,nrow=K,ncol=J)
for (k in 1:K)
{
ris$centroid[k,]=(t(ris$membership[,k]^m)%*%as.matrix(frj))/sum(ris$membership[,k]^m)
}
rownames(ris$centroid)=1:K
colnames(ris$centroid)=colnames(frj)
print(K)
print(round(ris$centroid,2))
print(classe.memb(ris$membership)$table.U)
print(ris$silinfo$avg.width)
K=K+1
}
this should be scheme clearly are determined centroid based on classe.memb.
classe.memb=function(U)
{
info.U=cbind(max.col(U),apply(U,1,max))
i=1
while (i = nrow(U))
{
if (apply(U,1,max)[i]0.5) info.U[i,1]=0
i=i+1
}
K=ncol(U)
table.U=matrix(0,nrow=K,ncol=4)
cl=1
while (cl = K)
{
table.U[cl,1] = length(which(info.U[info.U[,1]==cl,2]=.90))
table.U[cl,2] = length(which(info.U[info.U[,1]==cl,2]=.70)) -
table.U[cl,1]
table.U[cl,3] = length(which(info.U[info.U[,1]==cl,2]=.50)) -
table.U[cl,1] - table.U[cl,2]
table.U[cl,4] = sum(table.U[cl,])
cl = cl+1
}
rownames(table.U) = c(1:K)
colnames(table.U) = c(Alto, Medio, Basso, Totale)
out=list()
out$info.U=round(info.U,2)
out$table.U=table.U
return(out)
--
View this message in context:
http://r.789695.n4.nabble.com/clustering-fuzzy-tp3262027p3262027.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] (no subject)
-- View this message in context: http://r.789695.n4.nabble.com/no-subject-tp3262024p3262024.html Sent from the R help mailing list archive at Nabble.com. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] multivariate regression
On 04.02.2011 13:54, Deniz SIGIRLI wrote: How can I run multivariate linear regression in R (I have got 3 dependent variables and only 1 independent variable)? I tried lm function, but it gave different R2 and p values for every dependent variable. I need one R2 and p value for the model. Well, then you have to tell us how this one R^2 / p-value should be defined from your point of view. Uwe Ligges [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] spline interpolation
Date: Sat, 5 Feb 2011 18:08:43 +0100 From: cbelei...@units.it To: r-help@r-project.org Subject: Re: [R] spline interpolation Hi, Just pure curiosity: may I ask why you want to do spline interplation on fluorescence intensity as function of concentration? Particularly as it looks quite typical for an unknown problem's calibration plot? I would mention of course the suggestions given are not typically the end of the story. Instead of a bunch of splines, you may get a better result by writing a differential equation for the kinetics of the system you have, setting time derivatives to zero, and try to estimate things like lifetime-concentration combinations. Is that your point or you think this is a homework problem LOL? Claudia On 02/05/2011 03:29 PM, Asan Ramzan wrote: Hello R-help I have the following data for a standard curve concentration(nM),fluorescence 0,48.34 2,58.69 5,70.83 10,94.73 20,190.8 50,436.0 100,�957.9 � (1)Is there function in R�to plot a spline. (2)How can I interpolation,say 1000 point from 0nM-100nM and store this as a data frame of concentration,fluorescence (3)How can I modify�the code below so that instead of retrieving a concentration with the exact value of fluorescence, it gives me concentration for the value that is closest to that fluorescence. � subset(df,fluorescence==200.3456,select=concentration) [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Claudia Beleites Dipartimento dei Materiali e delle Risorse Naturali Università degli Studi di Trieste Via Alfonso Valerio 6/a I-34127 Trieste phone: +39 0 40 5 58-37 68 email: cbelei...@units.it __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] customizing names in an ordination biplot
On 04.02.2011 12:21, swertie wrote: Hello, I would like to change the names of the sites in an ordination biplot (resulting from the function rda in vegan). Can somebody give me some trick? Thank you very much See ?plot.cca how to put things together separately and control labels. Probably more easily (also for the rest of your analyses): Change the labels before applying the cca() function, then everything is displayed with the labels you specified before. Uwe Ligges __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Error in solve.default(inf, tol = tol.solve) :
On 04.02.2011 08:37, Marie-Line Glaesener wrote: Hello, I'm trying to run a lagsarlm (maximum likelihood estimation of a spatial lag model) in the spdep library ; but R gives following error message: Error in solve.default(inf, tol = tol.solve) : system is computationally singular: reciprocal condition number = 4.20137e-12 May we guess that you have almost perfectly correlated columns in your data - or factor levels with very few observations? Or many levels? But since we cannot know since we do not have your data. Read about how to invert matrices numerically and you will find that it is not easy if your matrix is almost singular. Uwe Ligges I get the same message when I try to run de lagsarlm with a bigger data set (4333 regions). The command I used: TestLag-lagsarlm( mean_price ~ transcount+ C1_5_1+ C1_5_2+ C1_5_3+ C1_5_4, data=DataB,transcecB.listw) summary(transcecB.listw) Characteristics of weights list object: Neighbour list object: Number of regions: 521 Number of nonzero links: 2904 Percentage nonzero weights: 1.069846 Average number of links: 5.573896 Link number distribution: 1 2 3 4 5 6 7 8 9 10 11 13 2 17 46 86 116 108 68 43 23 7 2 3 2 least connected regions: 12070513 12070504 with 1 link 3 most connected regions: 12060503 11040601 11020701 with 13 links Weights style: W Weights constants summary: n nn S0 S1 S2 W 521 271441 521 200.3434 2169.626 If any further information is needed, please let me know. I hope this is the right platform to get some help in solving this problem. Best regards, Marie-Line Glaesener Phd. Student Unité de Recherche IPSE (Identités. Politiques, Sociétés, Espaces) Laboratoire de Géographie et Aménagement du Territoire UNIVERSITÉ DU LUXEMBOURG CAMPUS WALFERDANGE Route de Diekirch / BP 2 L-7201 Walferdange Luxembourg www.geo.ipse.uni.luhttp://www.geo.ipse.uni.lu [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 summary
Perhaps you display an object but forgot to load the package that ships the summary method for the class of the object? I guess library(VGAM) summary(RES_TOBIT) should solve your problem. BTW: Please always report which package you are talking about (vghlm is not part of R base, obviously). Uwe Ligges On 03.02.2011 21:03, Hongwei Dong wrote: Hi, dear R users, I'm running a Tobit regression and using summary to display model results. It used to work. But this time, I kept getting this: summary(RES_TOBIT) Length Class Mode 1 vglm S4 Could anyone tell me how to solve this problem, and why I got this problem? Thanks. Gary [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 with Help
On 02.02.2011 21:15, John Filben wrote: I have recently been reading several books on data mining which contain a few data sets. The books offer some perspective on model choices, tuning decisions, result interpretation. Are there any good resources that can walk me through the thought process of an experienced data miner with the datasets that are included for packages such as rpart kmeans. kmeans is not a package - at least not on CRAN. I guess you are talking about a kmeans function in anotehr package. I think you should read two kinds of books: 1. books about the methodology and 2. book about using R. If you understood both, it is not too hard to put the two things together. There is a web page for R related book at http://www.r-project.org/ I currently use MS Access VBA for a bulk of my data manipulation and then use R Commander and Rattle for my statistical analysis. Thank you. Many of us will certainly think that neither MS Access VBA nor the R Commander (the latter is really nice for teaching certain kinds of non-statistics students) will be used by an experienced data miner. Best, Uwe Ligges Regards, John John Filben Cell Phone - 773.401.2822 Email - johnfil...@yahoo.com [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] drawing from one cell to another using layout() - possible?
On 02.02.2011 16:49, Mark Heckmann wrote: Is it possible to cross the cell boundaries set by layout using base graphics? I.e. I want to draw e.g. a line from one layout cell to another. Is there a way to do that? layout(matrix(c(1,2), byrow=TRUE, ncol=2)) plot.new() text(0,0,paste(rep(a, 200), collapse=), xpd=T) layout.show(2) I would like the a's to not end at the layout borders of the left cell. Use xpd=NA rather than xpd=TRUE. See ?par for the reason why. Uwe Ligges Thanks in advance, Mark PS. I need to use base not groid graphics, though it may be simpler... __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] problem with parLapply from snow
I just tried with SOCK rather than MPI on my machine and it worked. Just
ask R to
update.packages(checkBuilt=TRUE)
and try again. If you have the same problems, try SOCK rather than MPI.
If MPI is the culprit, please debug your MPI setup.
Uwe Ligges
On 03.02.2011 03:44, mark.pal...@csiro.au wrote:
Hi,
The following function use to work, but now it doesn't giving the error
CallSnow(, 100)
Using snow package, asking for 2 nodes
2 slaves are spawned successfully. 0 failed.
Error in checkForRemoteErrors(val) :
2 nodes produced errors; first error: no applicable method for 'lapply' applied to an
object of class list
.
Where this is the problem line of code yusum- parLapply(cl, yu, sum)
The function is
CallSnow- function (Nnodes = 2, Nsamples = Nnodes)
{
require(Rmpi)
require(snow)
cat(Using snow package, asking for , Nnodes, nodes \n)
cl- makeCluster(Nnodes, type=MPI)
on.exit(stopCluster(cl))
#print(do.call(rbind, clusterCall(cl, function(cl) Sys.info()[nodename])))
#
## uses RSPRNG if there
#
#clusterSetupSPRNG(cl)
clusterSetupRNGstream(cl, seed = rep(123456, 6))
yu- clusterCall(cl, runif, Nsamples)
yusum- parLapply(cl, yu, sum)
print(yusum)
yn- clusterCall(cl, rnorm, Nsamples)
print(yn)
return()
}
This is under R-2.12.1. on a windows Xp machine. This function still runs
satisfactorily on a machine running r-2.11.1 under Suse10.3/sles.
Thanks
Mark Palmer
Senior Statistician
CSIRO Mathematics, Informatics and Statistics
Phone: +61 8 9333 6293 | Fax: +61 8 9333 6121 | Mobile: 0427502353
mark.pal...@csiro.aumailto:mark.pal...@csiro.au | www.csiro.au |
www.csiro.au/cmishttp://www.csiro.au/cmis
Address: Private bag 5, Wembley, WA 6913, Australia
PLEASE NOTE
The information contained in this email may be confidential or privileged. Any
unauthorised use or disclosure is prohibited. If you have received this email
in error, please delete it immediately and notify the sender by return email.
Thank you. To the extent permitted by law, CSIRO does not represent, warrant
and/or guarantee that the integrity of this communication has been maintained
or that the communication is free of errors, virus, interception or
interference.
Please consider the environment before printing this email.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Getting variable names in function output
I realy think you are on a completely non-R way of implementing
something that could be done much better. But since I do not know what
you are going to do, here is the answer for your question:
Just add
names(end) - end
one line before your return() statement.
Uwe Ligges
On 03.02.2011 01:55, Sébastien Bihorel wrote:
Dear R-users,
I would like to have some advises about a problem illustrated by the
following snippet. Within myf, I need to evaluate a piece of R code that is
passed as a character argument and then return the objects that are created
by this code. The difficulty comes from the fact that the content of the
code is variable and unknown to me (obviously not in this illustration!).
With the following script, myf returns the content of the objects created by
the code, but I could not find a way to get the names of the objects in
addition to their contents.
Any suggestion will be welcome
Sebastien
rm(list=ls(all=T))
code1- paste('l- list(a=1,b=2,c=3)',
'm- matrix(1:15,nrow=3)', sep='\n')
code2- 'g- rnorm(100,10,2)'
myf- function(code=NULL){
start- ls(all=T)
eval(parse(text=code))
end- ls(all=T)
end- end[-which(end=='start')]
return(lapply(end[!end%in%
start],function(x) eval(parse(text=x
}
myf(code1)
myf(code2)
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Lightweight store plots as images.
Dear all I would like to store a few hundrends images to my hard disk so I need
your comment that is the most efficient way in R to do
for i in c(1:100){
A. create plot with title # I do not want to give output to screen\
B. store it as an image using plot's title as filename
C. Destroy that 'object'
}
I do not have any experience so can I have some advice from your side?
Best Regards
Alex.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Seeking help to define method for '+'
Dear all, I am trying to define + method for my newly defined s4 class which is as follows: setClass(Me, sealed=F,representation(x1 = numeric, x2 = character)) new1 - new(Me, x1=2, x2=comment1) new2 - new(Me, x1=3, x2=comment1) setMethod(+, Me, definition=function(x,y) cat(x@x1 + y@x1, \n)) However while am trying to set method, it fails: Error in conformMethod(signature, mnames, fnames, f, fdef, definition) : in method for '+' with signature 'e1=Me': formal arguments (e1 = Me, e2 = Me) omitted in the method definition cannot be in the signature I am just started learning s4 class system therefore please forgive me if I have some any bad mistakes. Anyway I would really appreciate if someone point me how should I proceed. Thanks and reagrds, [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] beanplot
Hello, sample1, sample2 and sample3 are vectors of numbers. The boxplot works properly but beanplot gives the following error for following run: boxplot(sample1, sample2, sample3, ylim = ylim, main = boxplot, names = 1:3) beanplot(sample1, sample2, sample3, ylim = ylim, main = beanplot, col = c(#CAB2D6, #33A02C, #B2DF8A), border = #CAB2D6, names=c(A,B,C)) log=y selected Warning message: In plot.window(xlim = c(0.5, 3.5), ylim = c(0, 100), log = y) : nonfinite axis limits [GScale(-inf,2,2, .); log=1] I cannot duplicate the error and warning with the following set of commands: par(mfrow = c(1, 2), mai = c(0.5, 0.5, 0.5, 0.1)) mu - 2 si - 0.6 c - 500 mu1 = 60 bimodal1 - c(rnorm(100, 60, 10), rnorm(900, 30, 10)) bimodal2 -c(rnorm(999, 60, 10), rnorm(1, 30, 10)) bimodal3 -c(rnorm(500, 60, 10), rnorm(500, 30, 10)) ylim - c(0, 100) boxplot(bimodal1, bimodal2, bimodal3, ylim = ylim, main = boxplot, names = 1:3) beanplot(bimodal1, bimodal2, bimodal3, ylim = ylim, main = beanplot, col = c(#CAB2D6, #33A02C, #B2DF8A), border = #CAB2D6, names=c(A,B,C)) Here, bimodal1, bimodal2 and bimodal3 are vectors with numbers with the same range of values as above for sample1,sample2 and sample3. Please let me know how to fix the error. Thank you very much, Paul __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] survfit - help avoiding re-inventing the wheel.
Greetings, R-ians: I would like to calculate the loglikelihood based on a Surv object and non-mle parameter estimates. I already have the mles as provided by survreg( Surv(time, time2, event, type) ). I can compute a loglikelihood based on a given density, censoring types, and parameter values, but find myself doing what must already be available, if I knew how to access it. The reason for wanting to compute a loglikelihood at non-mle values is to estimate the log likelihood ratios. Any help is much appreciated. Thanks. Charles Annis, P.E. charles.an...@statisticalengineering.com 561-352-9699 http://www.StatisticalEngineering.com [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Seeking help to define method for '+'
On 02/05/2011 02:19 PM, Bogaso Christofer wrote: Dear all, I am trying to define + method for my newly defined s4 class which is as follows: setClass(Me, sealed=F,representation(x1 = numeric, x2 = character)) new1 - new(Me, x1=2, x2=comment1) new2 - new(Me, x1=3, x2=comment1) setMethod(+, Me, definition=function(x,y) cat(x@x1 + y@x1, \n)) However while am trying to set method, it fails: Error in conformMethod(signature, mnames, fnames, f, fdef, definition) : in method for '+' with signature 'e1=Me': formal arguments (e1 = Me, e2 = Me) omitted in the method definition cannot be in the signature From getGeneric(+) standardGeneric for + defined from package base belonging to group(s): Arith function (e1, e2) standardGeneric(+, .Primitive(+)) environment: 0x19a9918 Methods may be defined for arguments: e1, e2 Use showMethods(+) for currently available ones. you can see that the generic is defined to dispatch on two arguments in a signature function(e1, e2). So setMethod(+, c(Me, Me), function(e1, e2) e1@x1 + e2@x1) but S4 also has group generics, and setMethod(Arith, c(Me, Me), function(e2, e2) callGeneric(e1@x1, e2@x1)) will define not just + but also - and other functions in the 'Arith' group. See ?GroupGenericFunctions Martin I am just started learning s4 class system therefore please forgive me if I have some any bad mistakes. Anyway I would really appreciate if someone point me how should I proceed. Thanks and reagrds, [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Computational Biology Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 Location: M1-B861 Telephone: 206 667-2793 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Seeking help to define method for '+'
Thanks Martin for your help. Although your suggestion is working for somewhat restrictive case, I would like to make thing more generic. For example: setMethod(Arith, c(Me, Me), function(e1, e2) callGeneric(e1@x1, e2@x1)) with in (new1 + new2) is working perfectly, however if I want to add arbitrary number of objects then it fails: new1 + new2 - new1 Error in new1 + new2 - new1 : non-numeric argument to binary operator What can I do to handle this? Thanks and regards, -Original Message- From: Martin Morgan [mailto:mtmor...@fhcrc.org] Sent: 06 February 2011 03:57 To: Bogaso Christofer Cc: r-help@r-project.org Subject: Re: [R] Seeking help to define method for '+' On 02/05/2011 02:19 PM, Bogaso Christofer wrote: Dear all, I am trying to define + method for my newly defined s4 class which is as follows: setClass(Me, sealed=F,representation(x1 = numeric, x2 = character)) new1 - new(Me, x1=2, x2=comment1) new2 - new(Me, x1=3, x2=comment1) setMethod(+, Me, definition=function(x,y) cat(x@x1 + y@x1, \n)) However while am trying to set method, it fails: Error in conformMethod(signature, mnames, fnames, f, fdef, definition) : in method for '+' with signature 'e1=Me': formal arguments (e1 = Me, e2 = Me) omitted in the method definition cannot be in the signature From getGeneric(+) standardGeneric for + defined from package base belonging to group(s): Arith function (e1, e2) standardGeneric(+, .Primitive(+)) environment: 0x19a9918 Methods may be defined for arguments: e1, e2 Use showMethods(+) for currently available ones. you can see that the generic is defined to dispatch on two arguments in a signature function(e1, e2). So setMethod(+, c(Me, Me), function(e1, e2) e1@x1 + e2@x1) but S4 also has group generics, and setMethod(Arith, c(Me, Me), function(e2, e2) callGeneric(e1@x1, e2@x1)) will define not just + but also - and other functions in the 'Arith' group. See ?GroupGenericFunctions Martin I am just started learning s4 class system therefore please forgive me if I have some any bad mistakes. Anyway I would really appreciate if someone point me how should I proceed. Thanks and reagrds, [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Computational Biology Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 Location: M1-B861 Telephone: 206 667-2793 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Seeking help to define method for '+'
On 05/02/2011 5:58 PM, Bogaso Christofer wrote: Thanks Martin for your help. Although your suggestion is working for somewhat restrictive case, I would like to make thing more generic. For example: setMethod(Arith, c(Me, Me), function(e1, e2) callGeneric(e1@x1, e2@x1)) with in (new1 + new2) is working perfectly, however if I want to add arbitrary number of objects then it fails: new1 + new2 - new1 Error in new1 + new2 - new1 : non-numeric argument to binary operator What can I do to handle this? Your function returns a number, not a Me object, so R ends up with a number minus a Me.You could define what to do in that case, or have your function construct a Me object and return that, in which case you'll only see the signature (Me, Me). In the latter case you won't be able to handle things like (new1 + 1). Duncan Murdoch Thanks and regards, -Original Message- From: Martin Morgan [mailto:mtmor...@fhcrc.org] Sent: 06 February 2011 03:57 To: Bogaso Christofer Cc: r-help@r-project.org Subject: Re: [R] Seeking help to define method for '+' On 02/05/2011 02:19 PM, Bogaso Christofer wrote: Dear all, I am trying to define + method for my newly defined s4 class which is as follows: setClass(Me, sealed=F,representation(x1 = numeric, x2 = character)) new1- new(Me, x1=2, x2=comment1) new2- new(Me, x1=3, x2=comment1) setMethod(+, Me, definition=function(x,y) cat(x@x1 + y@x1, \n)) However while am trying to set method, it fails: Error in conformMethod(signature, mnames, fnames, f, fdef, definition) : in method for '+' with signature 'e1=Me': formal arguments (e1 = Me, e2 = Me) omitted in the method definition cannot be in the signature From getGeneric(+) standardGeneric for + defined from package base belonging to group(s): Arith function (e1, e2) standardGeneric(+, .Primitive(+)) environment: 0x19a9918 Methods may be defined for arguments: e1, e2 Use showMethods(+) for currently available ones. you can see that the generic is defined to dispatch on two arguments in a signature function(e1, e2). So setMethod(+, c(Me, Me), function(e1, e2) e1@x1 + e2@x1) but S4 also has group generics, and setMethod(Arith, c(Me, Me), function(e2, e2) callGeneric(e1@x1, e2@x1)) will define not just + but also - and other functions in the 'Arith' group. See ?GroupGenericFunctions Martin I am just started learning s4 class system therefore please forgive me if I have some any bad mistakes. Anyway I would really appreciate if someone point me how should I proceed. Thanks and reagrds, [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Seeking help to define method for '+'
On 02/05/2011 02:58 PM, Bogaso Christofer wrote:
Thanks Martin for your help. Although your suggestion is working for
somewhat restrictive case, I would like to make thing more generic. For
example:
setMethod(Arith, c(Me, Me), function(e1, e2) callGeneric(e1@x1,
e2@x1))
with in (new1 + new2) is working perfectly, however if I want to add
arbitrary number of objects then it fails:
new1 + new2 - new1
Error in new1 + new2 - new1 : non-numeric argument to binary operator
as written, new1 + new2 returns a numeric, so (new1 + new2) - new1 would
require a method with signature c(numeric, Me). But perhaps what you
want is an implementation along the lines of
setMethod(Arith, c(Me, Me), function(e1, e2) {
new(Me, x1=callGeneric(e1@x1, e2@x2)
x2=c(e1@x2, e2@x2))
})
i.e., return an instance of Me.
Martin
What can I do to handle this?
Thanks and regards,
-Original Message-
From: Martin Morgan [mailto:mtmor...@fhcrc.org]
Sent: 06 February 2011 03:57
To: Bogaso Christofer
Cc: r-help@r-project.org
Subject: Re: [R] Seeking help to define method for '+'
On 02/05/2011 02:19 PM, Bogaso Christofer wrote:
Dear all, I am trying to define + method for my newly defined s4 class
which is as follows:
setClass(Me, sealed=F,representation(x1 = numeric, x2 = character))
new1 - new(Me, x1=2, x2=comment1)
new2 - new(Me, x1=3, x2=comment1)
setMethod(+, Me, definition=function(x,y) cat(x@x1 + y@x1, \n))
However while am trying to set method, it fails:
Error in conformMethod(signature, mnames, fnames, f, fdef, definition) :
in method for '+' with signature 'e1=Me': formal arguments (e1 = Me,
e2 = Me) omitted in the method definition cannot be in the signature
From
getGeneric(+)
standardGeneric for + defined from package base
belonging to group(s): Arith
function (e1, e2)
standardGeneric(+, .Primitive(+))
environment: 0x19a9918
Methods may be defined for arguments: e1, e2
Use showMethods(+) for currently available ones.
you can see that the generic is defined to dispatch on two arguments in
a signature function(e1, e2). So
setMethod(+, c(Me, Me), function(e1, e2) e1@x1 + e2@x1)
but S4 also has group generics, and
setMethod(Arith, c(Me, Me),
function(e2, e2) callGeneric(e1@x1, e2@x1))
will define not just + but also - and other functions in the 'Arith'
group. See
?GroupGenericFunctions
Martin
I am just started learning s4 class system therefore please forgive me if
I
have some any bad mistakes. Anyway I would really appreciate if someone
point me how should I proceed.
Thanks and reagrds,
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Computational Biology
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109
Location: M1-B861
Telephone: 206 667-2793
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 maximum value of the legend on an image.plot and animation
Steve,
I think that your second question may be answered by
either setting the appropriate value for 'outdir' in
ani.options() or by adding an outdir=yourdir to the
argument list in saveHTML.
Peter Ehlers
On 2011-02-05 09:22, Uwe Ligges wrote:
On 01.02.2011 16:10, steve_fried...@nps.gov wrote:
Hello,
I'm doing the following:
library(ncdf)
library(fields)
library(animation)
saline- open.ncdf(salinity_1990.nc)
salt = get.var.ncdf(nc=saline, varid=Salinity)
# create an animation of the complete temporal domain in the ncdf file.
saveHTML({
for (i in 1:364) {
image.plot(salt[, , i])
}
}, img.name = salinity.img, imgdir = salinity_dir, htmlfile =
salinity.html, zmax = c(0, 80),
autobrowse = FALSE, title = TIME SALINITY PREDICTIONS,
description = c(This should plot 1 years daily salinity
predictions in Florida Bay)
)
Almost all of the data sets I work with are multi-temporal spatial forms.
Being able to view the time sequence is very important to us, The animation
procedure is very good. I appreciate its development.
Here are my questions.
1) I would like to find a procedure to set the maximum value of the legend
of an image.plot. Each time step has a different maximum value, but for
the animation to be valuable, I need to standardize these to a maximum such
that the scale is equal in each time step.
Specify the breaks manually using the argument breaks that can be
passed to ... in image.plot().
2) Secondly, when setting imgdir it only goes to a temp directory. If I
define a directory, the new directory path is added to the default temp
directory path. Is there a way to fix this so that the output directory is
truly defined in the saveHTML statements?
I don't see a way, so just ask the maintainer of the animation package
to allow for absolute paths or so.
Uwe Ligges
I'm working on a windows XP machine using R 2.12.1.
Thanks
Steve
Steve Friedman Ph. D.
Ecologist / Spatial Statistical Analyst
Everglades and Dry Tortugas National Park
950 N Krome Ave (3rd Floor)
Homestead, Florida 33034
steve_fried...@nps.gov
Office (305) 224 - 4282
Fax (305) 224 - 4147
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] very basic HLM question
Hi everyone, I need to get a between-component variance (e.g. random effects Anova), but using lmer I don't get the same results (variance component) than using random effects Anova. I am using a database of students, clustered on schools (there is not the same number of students by school). According to the ICC1 command, the interclass correlation is .44 ICC1(anova1) [1] 0.4414491 However, I cannot get the same ICC from the lmer output: anova2 - lmer(math ~ 1 + (1|schoolid), data=nels88) summary(anova2 - lmer(math ~ 1 + (1|schoolid), data=nels88)) Linear mixed model fit by REML Formula: math ~ 1 + (1 | schoolid) Data: nels88 AIC BIC logLik deviance REMLdev 1878 1888 -935.8 18751872 Random effects: Groups NameVariance Std.Dev. schoolid (Intercept) 34.011 5.8319 Residual 72.256 8.5003 Number of obs: 260, groups: schoolid, 10 Fixed effects: Estimate Std. Error t value (Intercept) 48.861 1.927 25.36 The intercept random effect is 34.011. If I compute the ICC manually I get: 34.011/(34.011+72.256) [1] 0.3200523 According to my Anova analysis, the between-component variance is 59.004. Does anyone know what the problem is? How can I get the 59.004 figure using R? -- Sebastián Daza sebastian.d...@gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] beanplot
On 02/06/2011 09:12 AM, Paul Ossenbruggen wrote: Hello, sample1, sample2 and sample3 are vectors of numbers. The boxplot works properly but beanplot gives the following error for following run: boxplot(sample1, sample2, sample3, ylim = ylim, main = boxplot, names = 1:3) beanplot(sample1, sample2, sample3, ylim = ylim, main = beanplot, col = c(#CAB2D6, #33A02C, #B2DF8A), border = #CAB2D6, names=c(A,B,C)) log=y selected Warning message: In plot.window(xlim = c(0.5, 3.5), ylim = c(0, 100), log = y) : nonfinite axis limits [GScale(-inf,2,2, .); log=1] I cannot duplicate the error and warning with the following set of commands: par(mfrow = c(1, 2), mai = c(0.5, 0.5, 0.5, 0.1)) mu- 2 si- 0.6 c- 500 mu1 = 60 bimodal1- c(rnorm(100, 60, 10), rnorm(900, 30, 10)) bimodal2-c(rnorm(999, 60, 10), rnorm(1, 30, 10)) bimodal3-c(rnorm(500, 60, 10), rnorm(500, 30, 10)) ylim- c(0, 100) boxplot(bimodal1, bimodal2, bimodal3, ylim = ylim, main = boxplot, names = 1:3) beanplot(bimodal1, bimodal2, bimodal3, ylim = ylim, main = beanplot, col = c(#CAB2D6, #33A02C, #B2DF8A), border = #CAB2D6, names=c(A,B,C)) Here, bimodal1, bimodal2 and bimodal3 are vectors with numbers with the same range of values as above for sample1,sample2 and sample3. Please let me know how to fix the error. Hi Paul, Looks like you have a zero value in one of your sample vectors and not in your bimodal vectors. Jim __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] error with source(): invalid 'times' value
This is somewhat fixed now in R-patched and R-devel, as of revision 54235. It won't die with an error, but it still might not be perfect. The problem is that the line #line 516 VolStocksDec2010.Rnw is taken as a statement by you that the next few lines are copied from line 516 and following in the VolStocksDec2010.Rnw file. If that's not true (e.g. that file doesn't exist any more, or has changed) you might not get what you want in the echoed code. I may still make some more changes: either stop Stangle from including those lines, or add an option to source() to get it to ignore them. The trouble is that those lines are often useful: they're how errors are reported relative to the original Rnw file, rather than the intermediate tangled file. I've added a note to ?source to point out that there might be a problem; I may just stop with that. Duncan Murdoch On 24/01/2011 7:04 PM, Duncan Murdoch wrote: On 11-01-24 5:09 PM, mat wrote: Le 24. 01. 11 20:43, Duncan Murdoch a écrit : On 11-01-24 12:07 PM, Matthieu Stigler wrote: hi I am seeing a strange behavior I can't understand... doing: source(/tmp/RFile.r,echo=TRUE) Error in rep.int(c(prompt.echo, continue.echo), c(leading, length(dep) - : invalid 'times' value traceback() 3: rep.int(c(prompt.echo, continue.echo), c(leading, length(dep) - leading)) 2: paste(rep.int(c(prompt.echo, continue.echo), c(leading, length(dep) - leading)), dep, sep = , collapse = \n) 1: source(/tmp/RFile.r, echo = TRUE) But the file I am trying to source is very simple... see: $ more /tmp/RFile.r ### ### chunk number 1: ### #line 516 VolStocksDec2010.Rnw path-~/Dropbox/FAO/Papers/Volatility only pathMarkov-~/Dropbox/FAO/Markov Model/ library(zoo) Any idea where it can come from? It works fine when echo=FALSE... I am using R 2.12, on Ubuntu Linux 10.4 (R from CRAN), full session info below. Should I rather send this to r-devel? There is no such version, but this looks like a bug that was fixed in 2.12.1. Are you using 2.12.0? (I might be wrong about the timing of the fix; if you're using 2.12.1, try 2.12.1-patched.) Indeed 2.12.1, sorry for imprecision! I will give a try to 2.12.1-patched, although I am not so sure how I can install it (should I compile) on linux... Bill Dunlap has already confirmed that this is not what was fixed (or what was fixed never made it into the sources). I'll get to it, but not for a couple of weeks. Duncan Murdoch thanks!! Duncan Murdoch Thanks a lot Matthieu sessionInfo() R version 2.12.1 (2010-12-16) Platform: i486-pc-linux-gnu (32-bit) locale: [1] LC_CTYPE=fr_CH.utf8 LC_NUMERIC=C [3] LC_TIME=fr_CH.utf8LC_COLLATE=fr_CH.utf8 [5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=fr_CH.utf8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=fr_CH.utf8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices datasets utils methods base loaded via a namespace (and not attached): [1] grid_2.12.1 lattice_0.19-17 Matrix_0.999375-45 [4] nnet_7.3-1 tsDyn_0.7-40tseries_0.10-23 [7] tseriesChaos_0.1-11 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 list within a list?
Hello, I am planning of building a list of lists specifically, my first list is some what of the sort: lidta - list(m, p, r, s, q, A, B) where A and B are matrices that may be of different number of rows . The number of rows in matrix A and matrix B depends on the the values of m. The question is I don;t know how to put all the 1000 or so of these lists into a 'mega' list. Can you help me? -- Thanks, Jim. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Lightweight store plots as images.
I would use a pdf file, see ?pdf for examples. This plots a bunch of plots to
a single file. You can also use tools like png (see ?png for examples) which
will create 1 file per plot. You can either specify a name each time, or set a
name pattern and have it fill in automatically. Either way the plots go
straight to files without plotting on the screen.
--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
project.org] On Behalf Of Alaios
Sent: Saturday, February 05, 2011 2:22 PM
To: R-help@r-project.org
Subject: [R] Lightweight store plots as images.
Dear all I would like to store a few hundrends images to my hard disk
so I need your comment that is the most efficient way in R to do
for i in c(1:100){
A. create plot with title # I do not want to give output to screen\
B. store it as an image using plot's title as filename
C. Destroy that 'object'
}
I do not have any experience so can I have some advice from your side?
Best Regards
Alex.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 list within a list?
How are your lists being created? You can add each list to a mega list in a for loop, or use lapply to run a function multiple times which outputs a list each time and these will automatically be put together into a mega list. If these don't work for you then tell us more about how you are creating the inner lists (sample code helps). -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare greg.s...@imail.org 801.408.8111 -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- project.org] On Behalf Of Jim Silverton Sent: Saturday, February 05, 2011 7:39 PM To: r-help@r-project.org Subject: Re: [R] A list within a list? Hello, I am planning of building a list of lists specifically, my first list is some what of the sort: lidta - list(m, p, r, s, q, A, B) where A and B are matrices that may be of different number of rows . The number of rows in matrix A and matrix B depends on the the values of m. The question is I don;t know how to put all the 1000 or so of these lists into a 'mega' list. Can you help me? -- Thanks, Jim. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 list within a list?
On 6/02/2011 3:38 p.m., Jim Silverton wrote:
Hello,
I am planning of building a list of lists specifically, my first list is
some what of the sort:
lidta- list(m, p, r, s, q, A, B)
where A and B are matrices that may be of different number of rows . The
number of rows in matrix A and matrix B depends on the the values of m.
The question is I don;t know how to put all the 1000 or so of these lists
into a 'mega' list.
Can you help me?
I use the following for this sort of thing.
megaList - vector(list, length = 1000)
testList - list(x=1:3, y=c(a,b))
for (i in 1:1000){
megaList[[i]] - testList
}
head(megaList)
David Scott
_
David Scott Department of Statistics
The University of Auckland, PB 92019
Auckland 1142,NEW ZEALAND
Phone: +64 9 923 5055, or +64 9 373 7599 ext 85055
Email: d.sc...@auckland.ac.nz, Fax: +64 9 373 7018
Director of Consulting, Department of Statistics
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] exact logistic regression
Dear Brad Dear Jinko Dear David Sorry for the noise out of nothing. It was because of my ignorance and misunderstanding of algorithms of elrm. elrm is great for the assessment of multivariable model where Stata simply runs out of memory fails, while Stata can make exact calculations on one variable. What scared on the first place was variation in p-values when changing 'iter' option. Now I work with small dataset of 73 trials with 59 successes and trying to handle with categorical independent variables. Could you please suggest which number of iterations and burnIn is better to choose? Could you please comment how to interpret p-value for joint effect? For instance, I couldn't get the example in Zamar D, McNeney B and Graham J. elrm: Software Implementing Exact-like Inference for Logistic Regression Models. Journal of Statistical Software 2007, 21(3). on page 9, where p value of 0.76555 is consistent with the model used. And finally, could you please tell if I have to cite numbers number of iterations and burnIn in publication? With best regards Denis P.S. Sorry for naive questions У Пят, 04/02/2011 у 14:52 -0800, Brad McNeney піша: You might get a more useful response if you CC the package maintainer, David Zamar, who I don't think is a member of this list. I wonder if you could elaborate on the reason (and provide a test data set) for your distrust of the output. Brad - Original Message - From: Den d.kazakiew...@gmail.com To: R-help r-help@r-project.org Sent: Friday, 4 February, 2011 8:44:04 AM Subject: [R] exact logistic regression Hate to say that, but it looks like Stata is way above R, considering exact logistic regression. To use elrm() I have to aggregate my data,which is really time consuming when I look for the way out through many variables. But the worst thing is that I am not not sure if I can trust to p-values in output. I would be happy,however, if someone could contradict me on this issue. With best regards Denis __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Gompert-Makeham and Siler models
Hello R users, I am absolutely new with R. I would like to ask could you help me to build (basic ideas) Gompertz-Makeham and Siler's mortality models in R? I now that there are some information about that in M. J. Crawley R book conerning survival analysis. Could you suggest more literature how to build above mentioned mortality models in R? Best regards, Sarunas -- View this message in context: http://r.789695.n4.nabble.com/Gompert-Makeham-and-Siler-models-tp3262317p3262317.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Creating covariance matrices for simple and complex factor structure
Hello all, I'm trying to create covariance matrices that, when processed via CFA methods (in the sem package) will produce exact fit with simple structure down to poor fit with cross loadings. What is the best way to do this? I don't really need to have the exact loop code, but maybe an explanation of how to make a few of the matrices or an explanation of the rationale behind performing this task would be of benefit. Thanks in advance. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] exact logistic regression
Dear Łukas Thank you very much ...FOR EACH ROW... That's cool data2elrm-cbind(mydata,n=1) With best regards Denis У Пят, 04/02/2011 у 22:16 +0100, Łukasz Ręcławowicz піша: 2011/2/4 Den d.kazakiew...@gmail.com To use elrm() I have to aggregate my data,which is really time consuming when I look for the way out through many variables. You don't have to do that. One exception is that the binomial response should be specified as success/trials, where success gives the number of successes and trials gives the number of binomial trials for each row of dataset. So... when one row = one trial: data2elrm-cbind(mydata,n=rep(1,dim(mydata)[1])) But the worst thing is that I am not not sure if I can trust to p-values in output. It depends, but it's always a guess. -- Miłego dnia __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Multiplying elements of vectors
Sorry I just sent a message, but I forgot to say thanks for any suggestion and help!! Mariana __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Multiplying elements of vectors
Hi guys: Sorry if this question is very basic. I’m learning basic matrix and vectors multiplication to develop a population matrix model for plants. I’m trying to multiply the elements of two vectors (each of the “x” values by each of the “y” values) to obtain a square matrix of xy values. f.e. x-seq(5,205) y-seq(5,20,5) stages-c(“Sdl”, “Juv”, “Ad1”, “Ad2”) If I just multiply xy as a matrix xy-matrix(x,y,nrow=4,ncol=4,dimnames=list(stages,stages)) I obtain this xy Sdl Juv A1 A2 Sdl 5 10 15 20 Juv 5 10 15 20 A15 10 15 20 A25 10 15 20 but what I want to obtain is this matrix SdlJuv A1A2 Sdl 2550 75 100 Juv 50100150 200 A17550 225 300 A2100 200300 400 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Passing Bablok
Hi Benjamin,
I was frustrated by the same thing. I pulled the papers last week and wrote
this little function. Produces the same results as MethVal, a commercial
clinical chemistry package. I have not tested it thoroughly but I hope that
helps you.
Dan Holmes, MD, Vancouver
#
# R code for Passing Bablok Regression with
#
# confidence intervals
#
# Ref: J Clin Chem Biochem Vol 26, 1988 pp 783-790 #
# Written by Daniel Holmes,MD
#
# Department of Pathology and Laboratory Medicine #
# University of British Columbia
#
# St. Paul's Hospital
#
# 1081 Burrard St.
#
# Vancouver, BC
#
#
isodd - function(N){
if (N%%2==0) {
ans-FALSE
} else if (N%%2==1) {
ans-TRUE
} else {
ans-Neither!
}
return(ans)
}
PB.reg-function(data){
x-data[,1]
y-data[,2]
#then we proceed as usual
lx-length(x)
l-choose(lx,2)
S-matrix(1:l,nrow=1,ncol=l)
for (i in 1:(lx-1)) {
for (j in (i+1):lx) {
# This expression generates the natural numbers from
# i and j avoiding the use of another counter
index-(j-i)+(i-1)*(lx-i/2)
S[index]-(y[i]-y[j])/(x[i]-x[j])
}
}
S.sort-sort(S)
S.sort-subset(S.sort,S.sort!=is.na(S.sort))
N-length(S.sort)
neg-length(subset(S.sort,S.sort0))
K-floor(neg/2)
if (isodd(N)) {
b-S.sort[(N+1)/2+neg/2]
} else {
b-sqrt(S.sort[N/2+K]*S.sort[N/2+K+1])
}
a-median(y-b*x)
#CI of b
C.gamma-qnorm(0.975)*sqrt(lx*(lx-1)*(2*lx+5)/18)
M1-round((N-C.gamma)/2)
M2-N-M1+1
b.lower-S.sort[M1+K]
b.upper-S.sort[M2+K]
#CI of a
a.lower-median(y-b.upper*x)
a.upper-median(y-b.lower*x)
result-list(intercept=a,intercept.CI=c(a.lower,a.upper),slope=b,slope.CI=c(b.lower,b.upper))
return(result)
}
#Example application to mock data
x-seq(0:50)
y-5*x+rnorm(50,0,10)
data-as.data.frame(cbind(x,y))
reg-PB.reg(data)
reg
#Plot the regression
plot(x,y,pch=21,bg=gray,main=Passing Bablok Regression)
abline(reg$intercept,reg$slope,col=blue)
abline(reg$intercept.CI[1],reg$slope.CI[1],col=red,lty=2)
abline(reg$intercept.CI[2],reg$slope.CI[2],col=red,lty=2)
correl-paste(R = ,round(cor.test(x,y)$estimate,3),sep=)
slope-paste(Slope = ,round(reg$slope,2),
[,round(reg$slope.CI[1],2),,,round(reg$slope.CI[2],2),],sep=)
intercept-paste(Intercept = ,round(reg$intercept,2),
[,round(reg$intercept.CI[1],2),,,round(reg$intercept.CI[2],2),],sep=)
text-c(correl,slope,intercept)
legend(topleft,text,title=Regression Data,inset = .05)
--
View this message in context:
http://r.789695.n4.nabble.com/Passing-Bablok-tp886121p3262090.html
Sent from the R help mailing list archive at Nabble.com.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Multiplying elements of vectors
If the seq(5,205) was a typo, and should have been seq(5,20,5), then what you're looking for is the outer product of x and y: x = seq(5,20,5) y = seq(5,20,5) x %o% y [,1] [,2] [,3] [,4] [1,] 25 50 75 100 [2,] 50 100 150 200 [3,] 75 150 225 300 [4,] 100 200 300 400 outer(x,y) [,1] [,2] [,3] [,4] [1,] 25 50 75 100 [2,] 50 100 150 200 [3,] 75 150 225 300 [4,] 100 200 300 400 The outer() function will accepts a FUN= argument which defaults to '*'. - Phil Spector Statistical Computing Facility Department of Statistics UC Berkeley spec...@stat.berkeley.edu On Sat, 5 Feb 2011, Mariana Martinez-Morales wrote: Hi guys: Sorry if this question is very basic. I’m learning basic matrix and vectors multiplication to develop a population matrix model for plants. I’m trying to multiply the elements of two vectors (each of the “x” values by each of the “y” values) to obtain a square matrix of xy values. f.e. x-seq(5,205) y-seq(5,20,5) stages-c(“Sdl”, “Juv”, “Ad1”, “Ad2”) If I just multiply xy as a matrix xy-matrix(x,y,nrow=4,ncol=4,dimnames=list(stages,stages)) I obtain this xy Sdl Juv A1 A2 Sdl 5 10 15 20 Juv 5 10 15 20 A15 10 15 20 A25 10 15 20 but what I want to obtain is this matrix SdlJuv A1A2 Sdl 2550 75 100 Juv 50100150 200 A17550 225 300 A2100 200300 400 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] very basic HLM question
2011/2/5 Sebastián Daza sebastian.d...@gmail.com:
Hi everyone,
I need to get a between-component variance (e.g. random effects Anova), but
using lmer I don't get the same results (variance component) than using
random effects Anova. I am using a database of students, clustered on
schools (there is not the same number of students by school).
According to the ICC1 command, the interclass correlation is .44
ICC1(anova1)
[1] 0.4414491
If you don't tell us exactly what model you are calculating in
anova1, how would we guess if there is something wrong?
Similarly, I get this
ICC1
Error: object 'ICC1' not found
so it must mean you've loaded a package or written a function, which
you've not shown us.
I googled my way to a package called multilevel that has ICC1, and
its code for ICC1 shows a formula that does not match the one you used
to calculate ICC from lmer.
function (object)
{
MOD - summary(object)
MSB - MOD[[1]][1, 3]
MSW - MOD[[1]][2, 3]
GSIZE - (MOD[[1]][2, 1] + (MOD[[1]][1, 1] + 1))/(MOD[[1]][1,
1] + 1)
OUT - (MSB - MSW)/(MSB + ((GSIZE - 1) * MSW))
return(OUT)
}
I'm not saying that's right or wrong, just not obviously identical to
the formula you proposed.
However, I cannot get the same ICC from the lmer output:
anova2 - lmer(math ~ 1 + (1|schoolid), data=nels88)
summary(anova2 - lmer(math ~ 1 + (1|schoolid), data=nels88))
Instead, do this (same thing, fits model only once):
anova2 - lmer(math ~ 1 + (1|schoolid), data=nels88)
summary(anova2)
Note that lmer is going to estimate a normally distributed random
effect for each school, as well as an individual observation random
effect (usual error term) that is assumed independent of the
school-level effect. What is anova1 estimating?
--
Paul E. Johnson
Professor, Political Science
1541 Lilac Lane, Room 504
University of Kansas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Creating covariance matrices for simple and complex factor structure
At 6:00 PM + 2/5/11, LOUDERMILK, BRANDON wrote: Hello all, I'm trying to create covariance matrices that, when processed via CFA methods (in the sem package) will produce exact fit with simple structure down to poor fit with cross loadings. What is the best way to do this? I don't really need to have the exact loop code, but maybe an explanation of how to make a few of the matrices or an explanation of the rationale behind performing this task would be of benefit. Thanks in advance. [[alternative HTML version deleted]] Brandon, See some of the simulation functions in the psych package. e.g., sim.congeneric, sim.structure, sim. circumplex, sim.simplex, etc. Bill __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] different results in MASS's mca and SAS's corresp
On Sat, Feb 5, 2011 at 9:19 AM, David Winsemius dwinsem...@comcast.net wrote: On Feb 4, 2011, at 7:06 PM, Gong-Yi Liao wrote: Dear list: I have tried MASS's mca function and SAS's PROC corresp on the farms data (included in MASS, also used as mca's example), the results are different: R: mca(farms)$rs: 1 2 1 0.059296637 0.0455871427 2 0.043077902 -0.0354728795 3 0.059834286 0.0730485572 4 0.059834286 0.0730485572 [snip] And in SAS's corresp output: Row Coordinates Dim1 Dim2 1 1.0607 -0.8155 2 0.7706 0.6346 3 1.0703 -1.3067 4 1.0703 -1.3067 5 0.2308 0.9000 [snip] Is MASS's mca developed with different definition to SAS's corresp ? No, it's just that the values can only be defined up to a scaling factor (the same situation as with eigenvector decompostion). Take a look at the two dimensions, when each is put on the same scale: cbind(scale(rmca$D1),scale(smca$Dim1) ) [,1] [,2] [1,] 1.2824421 1.28242560 [2,] 0.9316703 0.93168561 [3,] 1.2940701 1.29403231 [4,] 1.2940701 1.29403231 [5,] 0.2789996 0.27905048 cbind(scale(rmca$D2),scale(smca$Dim2) ) [,1] [,2] [1,] 1.06673426 -1.06677626 [2,] -0.83006158 0.83012474 [3,] 1.70932841 -1.70932351 [4,] 1.70932841 -1.70932351 [5,] -1.17729983 1.17729909 David Winsemius, MD West Hartford, CT When I came to David's comment, I understood the theory, but not the numbers in his answer. I wanted to see the MASS mca answers match up with SAS, and the example did not (yet). Below see that if you scale the mca output, and then multiply column 1 of the scaled results by 0.827094, then you DO reproduce the SAS column 1 results exactly. Just rescale item 1 in mca's first column to match the SAS output. Repeat same with column 2, multiply -0.7644828, and you reproduce column 2. rmca - mca(farms) scalermca - scale(rmca$rs) scalermca[1,] 12 1.282442 1.066734 1.0607/1.282442 [1] 0.827094 -0.8155/1.06673426 [1] -0.7644828 cbind(scalermca[,1] * 0.827094, scalermca[,2] * -0.7644828) [,1][,2] 1 1.06070017 -0.8154 2 0.77057891 0.63456780 3 1.07031764 -1.30675217 4 1.07031764 -1.30675217 5 0.23075886 0.90002547 6 0.6943 0.60993995 7 0.10530240 0.78445402 8 -0.27026650 0.44225049 9 0.13426089 1.15670532 10 0.11861965 0.64778456 11 0.23807570 1.21775202 12 1.01156703 -0.01927226 13 0.28051938 -0.59805897 14 -1.17343686 -0.27122981 15 -0.83838041 -0.64003061 16 -0.05453708 -0.22925816 17 -0.91732401 -0.49899374 18 -0.92694148 -0.00774156 19 -1.30251038 -0.34994509 20 -1.30251038 -0.34994509 So, that does reproduce SAS exactly. And I'm a little frustrated I can't remember the matrix command to get that multiplication done without cbinding the 2 columns together that way. Question: I don't use mca, but to people who do, how are results supposed to be scaled? Is there a community accepted method or is every user on his/her own to fiddle up the numbers however? -- Paul E. Johnson Professor, Political Science 1541 Lilac Lane, Room 504 University of Kansas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] boostrap an nls regression
Hello! On Thu, Feb 3, 2011 at 6:27 AM, m234 mhairialexan...@gmail.com wrote: II functional response for 2 data sets: nls(eaten~(a*suppl)/(1+a*h*suppl) where eaten is the number of prey eaten by a predator and suppl is the number of prey initially supplied to the same predator. I have parameter estimates of 'a' and 'h' for the two populations studied and would like to know if there is a significant different in the estimates between the two i.e. a1 vs a2, h1 vs h2. I would like to bootstap the data to get multiple (~1000) estimates and compare then via a ttest or equivalent. Is it possible to do this and obtain multiple estimations which Hi, Please read the posting guide--a complete example of your code a link to a data frame will get you more helpful answers. I can tell you how to do bootstrap re-sampling on a model that is estimated from one data set. I wrote up some notes on that last summer (go to middle of this http://pj.freefaculty.org/R/SummerCamp2010/PJ-Lectures/functions1.pdf). I have other intro R material floating about under the R part of that link. But I'm a bit stumped by your question for statistical reasons-- nothing to do with R. From a statistical point of view, how do you test the difference between coefficients from 2 different fitted nls models, which are based on separate data samples? I don't think that's a trivial stat question, even if you have large samples and you only need to calculate that estimated difference one time. Then I wonder, how would a person do re-sampling when there are 2 data sets involved. If you know of a cite on that, I'd like to see it. One avenue I'd consider is to stack the 2 data sets together and rewrite my nls to estimate one equation with indicator variables (dummy variables) to separate the estimates for the 2 separate equations. But there would be some pooling tests you'd have to run first, to justify the idea that the 2 sets of data belong in the same analysis in the first place. Know what I mean? Suppose eaten is one long column, for both sets combined, but create suppl1 and suppl2 that are 0 for the other data set's cases, but suppl for the right one. Fit this: combmod - nls(eaten~(a1*suppl1)/(1+a1*h1*suppl1) + (a2*suppl2)/(1+a2*h2*suppl2)) This would conceivably allow comparison of a1 and a2. I think. I'm trying to remember sampling theory on nls. Well, in summary, I think you've got a harder stat problem than you have R problem. If you write out the code you use for a whole exercise to do this 1 time, we might see what to do. But remember to post the full working example--as much as you have, anyway. -- Paul E. Johnson Professor, Political Science 1541 Lilac Lane, Room 504 University of Kansas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Legend outside the plot? xpd?
Hi All, BG: Will try be brief. I'd like 3 graphs on a page (below each other mfrow=c(3,1)), saved to pdf. The three plot data on the same subject so I'm having one legend, to the right of the center graph. I'm using mar=c(5,15,4,15) to bring the sides in so that the graphs are square and not stretched wide. To have the graph to the side I'm thinking xpd=T. Each graph has a number of points and lines overlaid, so plot ();lines(),lines(),lines() etc. The data is somewhat of a subset of a larger set, so I'm limiting what is being displayed. The lines plot trends of the wider data. P: With xpd=T, the lines are going right out of the graph boxes to the outer limit of the plot boundaries, as would be the intended behaviour of xpd. Goes without saying that this is undesirable. Q: Is there a better way to achieve what I want? At this stage I have xpd=T in the par(), and xpd=F in the plot() commands and all the line commands, just so I can get the legend to actually print... Aside: Given the way I've done this, the lines() are all clipped RIGHT at the limit of the plot box. So given the plot is called first these lines leave little coloured dashes on the black of the plot box. To sort this I've called box() at the end. This all seems very redundant? Thanks Matt -- mattcst...@gmail.com [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
