[R] plotting wind rose data
Hi List, I am trying to create a spatial representation of some wind data. I have the season, frequency, strength and direction of the wind from 10 different locations, the coverage of the area that I am interested in is not 100% there are small gaps in my coverage due to the location of the weather stations. I am trying to create a series of wind maps e.g. the Prevailing Winds, the maximum seasonal wind, etc. Could any body recommend any R-packages that would cover this type problem/issue? Any views expressed in this message are those of the individual sender, except where the sender specifically states them to be the views of the Pinan Software __ 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] Unix batch to different nodes
Hello,
I am struggling with computing nodes in Unix.
I have the use of a Unix server that has 30 nodes and I would like to
batch scripts.
Here is an R example that results in 72 repeated tasks based on the 2
loops. If I wanted to send these out to the different nodes, each node
has 1 task and the remaining 42 tasks are put in a queue, what would I
do?
#Example:
proj.dir-C:/test/
setwd(proj.dir)
group.vars -
c(biol1,biol4,biol5,biol6,biol12,biol15,biol16,biol17)
species - c(aa, bb, cc, dd, ee, ff, gg, hh, ii)
for (gv in 1:length(group.vars)){
for (sp in 1:length(species)){
zzz - file(paste(group.vars[gv],species[sp],.TXT,sep=),w)
cat(paste(group.vars[gv]),file=zzz)
cat(paste(species[sp]),file=zzz)
close(zzz)
}
}
I know it is a silly script test but the division is similar to my
real task which is repeated 459 times and each run takes 18 hours
At this point I am confused what should be written in Unix and then at
what point should I call R. I have read an abundance of things but I
feel like I am missing something essential. Or perhaps I have read all
the wrong things.
Thanks,
Daisy
--
Daisy Englert Duursma
Room E8C156
Dept. Biological Sciences
Macquarie University NSW 2109
Australia
__
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] interactive session
Hi guys,
My concern is to create an automated process from the beginning to the end. I
want to copy all my code together in one piece and paste it to R console and
sit
back and relax :) except one moment in which the program should ask me to enter
a number.. and only then (only after getting a valid number from me) it should
continue to read and process the rest of the code. I tried lots of things
(readline, readLines, scan, interactive, ask, switch,...) and read manuals and
searched help forums.. I found several similar questions but never a satisfying
answer.. so now it became a challenge.. any idea how to overcome this problem
(R
2.11.1 for Windows)? (an example code is below)
Best,
Fatih
library(gtools)
library(YaleToolkit)
library(xts)
### start of my wrong function!
f-function(w){
w-readline(which data? )
w-as.numeric(w)
ifelse(is.numeric(w)==TRUE, w, f())
}
f()
# end of my wrong function
v- ## and output of my function should be a v for example which I can use it
in the next line (v-w or something like that??)
##the rest works fine
p-paste(t, v, .txt, sep = )
t-read.table(p, header=FALSE, sep=\t, dec=,,
blank.lines.skip=FALSE)
rownames(t)-as.Date(t[,1],%d.%m.%Y)
colnames(t)-c(date,start,high,low,end,w.average,lot,
volume)
x-as.xts(t)
whatis(x)
.
.
[[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] R-squared in Robust Regression
May I know how to find the R-squared for robust regression model? Thank you. Hock Ann [[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] Several Lattice plots in one Plot
Hi all, I've been trying for hours, but I do not find a Solution. I want to plot 12 variables over time in separate diagrams in one plot/window using lattice. Two columns, six rows. I used print with the split command, but the graphics are getting really small. Can someone please help me. Following data example: dta = data.frame( day=c(1,2,3,4,5,6,7), var11=c(1,2,2,4,5,3,2), var12=c(1,2,2,4,5,3,2), var13=c(1,2,2,4,5,3,2), var14=c(1,2,2,4,5,3,2), var15=c(1,2,2,4,5,3,2), var16=c(1,2,2,4,5,3,2), var17=c(1,2,2,4,5,3,2), var18=c(1,2,2,4,5,3,2), var19=c(1,2,2,4,5,3,2), var10=c(1,2,2,4,5,3,2), var11=c(1,2,2,4,5,3,2), var12=c(1,2,2,4,5,3,2)) Any idea how I can plot varXX over day like xyplot(var11 ~ day, data=dta, type='b', scales=list(cex=0.5), xlab=NULL, ylab=list(cex=0.5)) and get readable graphs? Thanks in advance. Marcus __ 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] Installing sp and rgdal libraries
Hi, I am still a bit of a beginner with R but I had it all working spgrass6, RMySQL, sp, etc... Then I wanted to install spatstat, this required updating R, all was going fine until I tried to connect to grass gis using spgrass6 library, it wouldn't let be install the package due to a dependency with rgdal. rgdal will not install as the sp package is 1:0.9-56-2 and sp requires a higher version. I can not find a higher version. Here is the error upon install.packages(rgdal) ** preparing package for lazy loading Error : package 'sp' 0.9-56 was found, but = 0.9.60 is required by 'rgdal' ERROR: lazy loading failed for package ârgdalâ * removing â/home/gary/R/x86_64-pc-linux-gnu-library/2.10/rgdalâ The downloaded packages are in â/tmp/RtmpCtm6wH/downloaded_packagesâ Warning message: In install.packages(rgdal) : installation of package 'rgdal' had non-zero exit status I would really like to learn and use R as it has so much time saving potential as well as some great innovative functions which I can't find anywhere else. Especially as it can link to grassgis. Thank you for any help you can offer Gary -- View this message in context: http://r.789695.n4.nabble.com/Installing-sp-and-rgdal-libraries-tp2720313p2720313.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] Hashing a set
Lorenzo Isella lorenzo.isella at gmail.com writes: Dear All, I am given a time series such at, at every time t_i, I am given a set of data (every element of the set is just an integer number). What I need is an injective function able to map every set into a number (possibly an integer number, but that is not engraved in the stone). Does anybody know how to achieve that? In set theory you learn about the function that assigns to a set of integers (i1, ..., in) the integer p1^i1 * ... * pn^in, where p1, ... is the sequence of primes. In practice, unfortunately, this will lead to too large numbers. I would recommend the 'digest' package that provides hashing functions such as 'md5', etc., thanks to Dirk Eddelbuettel. It is injective enough and returns character strings of fixed length. Hans Werner Cheers Lorenzo __ 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] Several Lattice plots in one Plot
Marcus Drescher drescher at tum.de writes: I've been trying for hours, but I do not find a Solution. I want to plot 12 variables over time in separate diagrams in one plot/window using lattice. Two columns, six rows. I used print with the split command, but the graphics are getting really small. Can someone please help me. Following data example: dta = data.frame( day=c(1,2,3,4,5,6,7), var11=c(1,2,2,4,5,3,2), var12=c(1,2,2,4,5,3,2), var13=c(1,2,2,4,5,3,2), var14=c(1,2,2,4,5,3,2), var15=c(1,2,2,4,5,3,2), var16=c(1,2,2,4,5,3,2), var17=c(1,2,2,4,5,3,2), var18=c(1,2,2,4,5,3,2), var19=c(1,2,2,4,5,3,2), var10=c(1,2,2,4,5,3,2), var11=c(1,2,2,4,5,3,2), var12=c(1,2,2,4,5,3,2)) Any idea how I can plot varXX over day like xyplot(var11 ~ day, data=dta, type='b', scales=list(cex=0.5), xlab=NULL, ylab=list(cex=0.5)) and get readable graphs? Thanks in advance. Marcus There are surely several ways to do this but how about DTA - cbind(day = dta$day, stack(dta[, -1])) xyplot(values ~ day | ind, DTA, type = b, layout = c(2, 6)) for which you can add additional annotations as desired. By the way, do you realize that you have repeated column names in your data frame? HTH, Ken -- Ken Knoblauch Inserm U846 Stem-cell and Brain Research Institute Department of Integrative Neurosciences 18 avenue du Doyen Lépine 69500 Bron France tel: +33 (0)4 72 91 34 77 fax: +33 (0)4 72 91 34 61 portable: +33 (0)6 84 10 64 10 http://www.sbri.fr/members/kenneth-knoblauch.html __ 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] matrix plot
On 09/29/2010 08:55 PM, hairryharry wrote: Hi, Fairly new to R - have done basic plots but now faced with plotting a matrix/table of results -I know what I want but cannot find out how to do it. Basically have individual questions ( x) to which an organization can rate themselves 1-10 (y) what I want to show is a matrix/density type plot (like the matrix corrolation plots I have seen on R graph site) showing for eac question (x) a circle or shape of varying size and/or colour to represent the number of organisations rating themselves for each value of y. Hope this makes sense, any advice would be gratefully received. Hi Harry, The balloonplot function will probably give you the circles you want. If you want a matrix of colors, color2D.matplot might help, and if you just want a display of counts, try count.overplot. 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] plotting wind rose data
Hi:
Try this:
library(sos) # install first if you don't have it
findFn('wind rose')
You'll find there are several packages with functions to produce wind rose
plots.
HTH,
Dennis
On Wed, Sep 29, 2010 at 11:12 PM, David Potts dave.po...@pinan.co.ukwrote:
Hi List,
I am trying to create a spatial representation of some wind data.
I have the season, frequency, strength and direction of the wind from 10
different locations, the coverage of the area that I am interested in is
not 100% there are small gaps in my coverage due to the location of the
weather stations.
I am trying to create a series of wind maps e.g. the Prevailing Winds, the
maximum seasonal wind, etc.
Could any body recommend any R-packages that would cover this type
problem/issue?
Any views expressed in this message are those of the individual sender,
except where the sender specifically states them to be the views of the
Pinan Software
__
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.
[[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] plotting wind rose data
Thanks Denis, looks useful.
I not trying to generate a wind rose per data point,
I am trying to generate a raster map/shape file showing for example, the
Prevailing Winds at a point midway between to wind roses.
Hi:
Try this:
library(sos) # install first if you don't have it
findFn('wind rose')
You'll find there are several packages with functions to produce wind rose
plots.
HTH,
Dennis
On Wed, Sep 29, 2010 at 11:12 PM, David Potts
dave.po...@pinan.co.ukwrote:
Hi List,
I am trying to create a spatial representation of some wind data.
I have the season, frequency, strength and direction of the wind from 10
different locations, the coverage of the area that I am interested in is
not 100% there are small gaps in my coverage due to the location of the
weather stations.
I am trying to create a series of wind maps e.g. the Prevailing Winds,
the
maximum seasonal wind, etc.
Could any body recommend any R-packages that would cover this type
problem/issue?
Any views expressed in this message are those of the individual sender,
except where the sender specifically states them to be the views of the
Pinan Software
__
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.
Any views expressed in this message are those of the individual sender,
except where the sender specifically states them to be the views of the
Pinan Software
__
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) Change of email address
On 09/30/2010 03:27 AM, Kevin E. Thorpe wrote: Ted is well aware of how to change his list email. He was advising people on the list who who have his old email address in their address books to remove it. How can we be sure, Kevin? This may be a desperate cry for help from Ted, cast adrift from the security of familiar names and overlearned routines. Should we abandon our dear e-correspondent in his hour of need? No! We must rally round, in a sort of symbolic way, and send instructions. We're still here, Ted, and you're not alone. 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.
[R] igraph / eigenvector centrality score
Hi to all, I have two graphs with the same number of nodes but with different connectivities and also with a different number of clusters. The two graphs represent the same system under different conditions and then there is a one-to-one correspondence between a given node in the two graphs. It is correct to use the eigenvector centrality score as a measure of the relevance assumed by a given node in the two graphs ? If not, which is the better way to perform this comparison ? Thank you in advance -- View this message in context: http://r.789695.n4.nabble.com/igraph-eigenvector-centrality-score-tp2720401p2720401.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] calculating mean and s.d. from a two-column table
On 2010-09-27 15:20, Joshua Wiley wrote: Hi, Peter's suggestion is more general, but for just the weighted mean, there is a built in function you can use (I do not know of any basic weighted standard deviation or variance functions). dat- data.frame(age = 1:5, no = c(21, 31, 9, 12, 6)) weighted.mean(x = dat$age, w = dat$no) Best regards, Josh Hmisc has wtd.mean() and wtd.var() as well as a few other weighted stats. -Peter Ehlers On Mon, Sep 27, 2010 at 9:34 AM, Jonas Josefsson jo...@runtimerecords.net wrote: I have a two-column table as follows where age is in the 1st column and the number of individuals is in the 2nd. age;no 1;21 2;31 3;9 4;12 5;6 Can I use mean() and sd() to calculate the mean and standard deviation from this or do I have to manually multiplicate 21*1+31*2 etc. / N? __ 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] (OT) Change of email address
On 30-Sep-10 09:02:18, Jim Lemon wrote: On 09/30/2010 03:27 AM, Kevin E. Thorpe wrote: Ted is well aware of how to change his list email. He was advising people on the list who who have his old email address in their address books to remove it. How can we be sure, Kevin? This may be a desperate cry for help from Ted, cast adrift from the security of familiar names and overlearned routines. Should we abandon our dear e-correspondent in his hour of need? No! We must rally round, in a sort of symbolic way, and send instructions. We're still here, Ted, and you're not alone. Jim I am overwhelmed!!! Just for the record: Kevin got it exactly right, and any instructions for changing my subscribed addresses in r-help, etc., had already been executed by me. I was indeed broadcasting a notification, to any to whom it might apply (but whom I might not be explicitly aware of myself), that they should update their records. So, Jim, is this enough to make it sure? Or are people's depths of concern and worry so intense that there remain doubts as to my well-being which still need to be addressed? With thanks for the solidarity! Ted. E-Mail: (Ted Harding) ted.hard...@wlandres.net Fax-to-email: +44 (0)870 094 0861 Date: 30-Sep-10 Time: 10:36:41 -- XFMail -- __ 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] Installing sp and rgdal libraries
On Thu, 2010-09-30 at 00:40 -0700, Gary Nobles wrote: Hi, I am still a bit of a beginner with R but I had it all working spgrass6, RMySQL, sp, etc... Then I wanted to install spatstat, this required updating R, all was going fine until I tried to connect to grass gis using spgrass6 library, it wouldn't let be install the package due to a dependency with rgdal. rgdal will not install as the sp package is 1:0.9-56-2 and sp requires a higher version. I can not find a higher version. Here is the error upon install.packages(rgdal) ** preparing package for lazy loading Error : package 'sp' 0.9-56 was found, but = 0.9.60 is required by 'rgdal' ERROR: lazy loading failed for package ârgdalâ * removing â/home/gary/R/x86_64-pc-linux-gnu-library/2.10/rgdalâ The downloaded packages are in â/tmp/RtmpCtm6wH/downloaded_packagesâ Warning message: In install.packages(rgdal) : installation of package 'rgdal' had non-zero exit status I would really like to learn and use R as it has so much time saving potential as well as some great innovative functions which I can't find anywhere else. Especially as it can link to grassgis. Thank you for any help you can offer Gary When you updated R (to what version?) did you also update your package library? update.packages(checkBuilt = TRUE) ## add ask = FALSE if you just want it to do it ## update.packages(checkBuilt = TRUE, ask = FALSE) will normally suffice. The error message indicates you already have the sp package installed but the version you have is too old for the version of rgdal you are trying to install. You could also try install.packages(rgdal, depend = TRUE, lib = path/to/library) when doing the installation, but I don't recall now whether that will check the version of installed packages and install [new versions of] them if necessary. G -- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% Dr. Gavin Simpson [t] +44 (0)20 7679 0522 ECRC, UCL Geography, [f] +44 (0)20 7679 0565 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/ UK. WC1E 6BT. [w] http://www.freshwaters.org.uk %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% __ 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] plotting wind rose data
On 09/30/2010 04:12 PM, David Potts wrote: Hi List, I am trying to create a spatial representation of some wind data. I have the season, frequency, strength and direction of the wind from 10 different locations, the coverage of the area that I am interested in is not 100% there are small gaps in my coverage due to the location of the weather stations. I am trying to create a series of wind maps e.g. the Prevailing Winds, the maximum seasonal wind, etc. Could any body recommend any R-packages that would cover this type problem/issue? Hi David, While there are several packages that include plotting routines for wind roses, it looks to me as though you want to define a small number of vectors representing prevailing wind, etc., possible overlaying these on a plot. That might be a job for a circular plotting routine (e.g. polar.plot) rather than a wind rose. 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] Installing sp and rgdal libraries
hi thankyou for the help, I upgraded Ubuntu and it disabled some repositories, all fixed now -- View this message in context: http://r.789695.n4.nabble.com/Installing-sp-and-rgdal-libraries-tp2720313p2720493.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] Controling R from MS Access
Dear Felipe, In the past I have done something similar. I had Access to run an R script in batch mode. The script reads data from Access, processes the data and puts data back in Access tables. The user hit a button in Access and see a blank console window popping up. When the console window disappear, the script has completed. Best regards, Thierry ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek team Biometrie Kwaliteitszorg Gaverstraat 4 9500 Geraardsbergen Belgium Research Institute for Nature and Forest team Biometrics Quality Assurance Gaverstraat 4 9500 Geraardsbergen Belgium tel. + 32 54/436 185 thierry.onkel...@inbo.be www.inbo.be To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher The plural of anecdote is not data. ~ Roger Brinner The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey -Oorspronkelijk bericht- Van: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] Namens Felipe Carrillo Verzonden: woensdag 29 september 2010 22:44 Aan: rco...@mailman.csd.univie.ac.at; r-h...@stat.math.ethz.ch Onderwerp: [R] Controling R from MS Access HI: I've seen a few threads about this topic but still can't find a straightforward way on this. Is there a package that can control R within an access form. For example, I want to send a query to R, perform some statistics in R and send the output or summary back to Access and display it on a form. I can do this using Excel as the 'middleman' but was wondering if it can be done straight from access. I can open R and send data to Access but I prefer to do it the other way around so that I can automate the process with VBA. It is real easy to control R from Excel with RExcel but can't find any examples where Access is used. Thanks for any advice. Felipe D. Carrillo Supervisory Fishery Biologist Department of the Interior US Fish Wildlife Service California, USA __ 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. Druk dit bericht a.u.b. niet onnodig af. Please do not print this message unnecessarily. Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document. The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document. __ 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] Opening a .R file with R (Windows)
Hi One workaround in Windows You shall have SendTo folder located somewhere in system (depending on Windows version). Put a shortcut to R programme in that folder. Than you can click with right mouse button on *.r file, select send to and choose appropriate shortcut (in that case R) but you can use it for other programs too. Regards Petr r-help-boun...@r-project.org napsal dne 28.09.2010 19:22:22: Hi Kye, I have never gotten .R files to work quite like other types (e.g., double-clicking a .PDF) in Windows. AFAIK there is no simple way to do it, because you do not edit scripts directly in R (I am happy to be corrected if someone knows better). For general use, I would just open R first and then open the file, or if you just want to run the file, you can use R's batch mode from the Windows command prompt. Best regards, Josh On Tue, Sep 28, 2010 at 10:11 AM, Kye Gilder kye.gil...@gmail.com wrote: I am new to using R. I installed R on my computer (Windows) and everything things appears to be just fine. However, I have a simple script RTest.R that does a few simple calculations. When I double-click the RTest.R icon, I get an Information dialong box which says, ARGUMENT 'C:\Documents and Settings\kgilder\Desktop\RTest.R' __ignored__ . When I choose OK, R then opens, but it does not open or display the script. Any help or suggestions? When I open R and use File Open New Script and path to the file, it opens just fine. Regards, Kye [[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. -- Joshua Wiley Ph.D. Student, Health Psychology University of California, Los Angeles http://www.joshuawiley.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-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] ARIMA models - estimation prediction$se
Dear All Can anyone help me? fit - arima(USAccDeaths, order = c(0,1,1),seasonal = list(order=c(0,1,1 fit$sigma2 [1] 99346.89 So, the standard error for my first step prediction is sqrt(fit$sigma2)=315.1934 like predict(fit, n.ahead = 6)$se[1] predict(fit, n.ahead = 6) $pred Jan Feb Mar Apr May Jun 1979 8336.061 7531.829 8314.644 8616.869 9488.912 9859.757 $se Jan Feb Mar Apr May Jun 1979 315.4481 363.0056 405.0168 443.0622 478.0896 510.7203 And now, How can I calculate the standard errors for the 5 steps ahead? Thanks in advance. Ana __ 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] Fitting a half-ellipse curve
Hello Niklaus,
I'm not sure if the following is the sort of thing you are looking for (?)
You can fit an ellipse to your data using a deterministic least
squares method. The following is a function that I use to do this...
fit.ellipse - function (x, y = NULL)
{
# Least squares fitting of an ellipse to point data
# using the algorithm described in:
# Radim Halir Jan Flusser. 1998.
# Numerically stable direct least squares fitting of ellipses.
# Proceedings of the 6th International Conference in Central Europe
# on Computer Graphics and Visualization. WSCG '98, p. 125-132
#
# Adapted from the original Matlab code by Michael Bedward
# michael.bedw...@gmail.com
#
# Arguments:
# x, y - the x and y coordinates of the data points
#
# Returns a list with the following elements:
#
# coef - coefficients of the ellipse as described by the general
#quadratic: ax^2 + bxy + cy^2 + dx + ey + f = 0
#
# center - center x and y
#
# major - major semi-axis length
#
# minor - minor semi-axis length
#
EPS - 1.0e-8
dat - xy.coords(x, y)
D1 - cbind(dat$x * dat$x, dat$x * dat$y, dat$y * dat$y)
D2 - cbind(dat$x, dat$y, 1)
S1 - t(D1) %*% D1
S2 - t(D1) %*% D2
S3 - t(D2) %*% D2
T - -solve(S3) %*% t(S2)
M - S1 + S2 %*% T
M - rbind(M[3,] / 2, -M[2,], M[1,] / 2)
evec - eigen(M)$vec
cond - 4 * evec[1,] * evec[3,] - evec[2,]^2
a1 - evec[, which(cond 0)]
f - c(a1, T %*% a1)
names(f) - letters[1:6]
# calculate the center and lengths of the semi-axes
b2 - f[2]^2 / 4
center - c((f[3] * f[4] / 2 - b2 * f[5]), (f[1] * f[5] / 2 - f[2] *
f[4] / 4)) / (b2 - f[1] * f[3])
names(center) - c(x, y)
num - 2 * (f[1] * f[5]^2 / 4 + f[3] * f[4]^2 / 4 + f[6] * b2 -
f[2]*f[4]*f[5]/4 - f[1]*f[3]*f[6])
den1 - (b2 - f[1]*f[3])
den2 - sqrt((f[1] - f[3])^2 + 4*b2)
den3 - f[1] + f[3]
semi.axes - sqrt(c( num / (den1 * (den2 - den3)), num / (den1 *
(-den2 - den3)) ))
# calculate the angle of rotation
term - (f[1] - f[3]) / f[2]
angle - atan(1 / term) / 2
list(coef=f, center = center, major = max(semi.axes), minor =
min(semi.axes), angle = unname(angle))
}
There are quite a few functions / packages to draw ellipses in R, but
the following is will work directly with the output of the above
function...
get.ellipse - function ( fit, n=360 )
{
# Calculate points on an ellipse described by
# the fit argument as returned by fit.ellipse
#
# n is the number of points to render
tt - seq(0, 2*pi, length=n)
sa - sin(fit$angle)
ca - cos(fit$angle)
ct - cos(tt)
st - sin(tt)
x - fit$center[1] + fit$maj * ct * ca - fit$min * st * sa
y - fit$center[2] + fit$maj * ct * sa + fit$min * st * ca
cbind(x=x, y=y)
}
So if your data were in a matrix or data.frame X...
efit - fit.ellipse( X )
e - get.ellipse( efit )
plot(X)
lines(e, col=red)
Hope this helps,
Michael
On 29 September 2010 23:45, Niklaus Hurlimann niklaus.hurlim...@unil.ch wrote:
Dear mailing list,
I have following array:
X2 Y2
[1,] 422.7900 6.0
[2,] 469.8007 10.5
[3,] 483.9428 11.0
[4,] 532.4917 25.5
[5,] 596.1942 33.5
[6,] 630.8496 40.5
[7,] 733.2996 45.0
[8,] 946.4779 32.0
[9,] 996.8068 35.5
[10,] 1074.3310 23.0
I do afterwards the following:
plot.new()
plot.window(xlim=c(min(X1)-50,max(X1)+50),
ylim=c(min(Y1)-25,max(Y1)+25))
axis(2, cex.axis=1.2)
axis(1, cex.axis=1.2)
points(X2, Y2, type=p, pch=0, cex=1.2, col=black)
title(main=Dyke Geometry Along Strike, cex.main=1.2, font.main=4)
title(xlab=distance [m], cex.lab=1.2)
title(ylab=half-thickness [m], cex.lab=1.2)
box()
I would like curve fitting where I proceeded with a non
linear-regression with the according function below:
nls(formula = Y2 ~ abs(b*abs(1-X2^2/a^2)^(1/2)), start = list( a=282,
b=40))
The formula should give the y-positive part of an ellipse in my opinion
or I might be completely wrong.
In the end I would like to fit a curve of half an ellipse through the
data. I might could do this as well with a 2nd order polynomial fit. I
am grateful for any suggestions and comments to my problem.
Cheers
--
Niklaus Hürlimann
Doctorant-PhD
Université de Lausanne
Institut de Minéralogie et Géochimie
L'Anthropole
CH-1015 Lausanne
Suisse
E-mail: niklaus.hurlim...@unil.ch
Tel:+41(0)21 692 4452
[[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] Unix batch to different nodes
This, I think, all depends on how the server is set up.
On Thu, Sep 30, 2010 at 1:21 AM, Daisy Englert Duursma
daisy.duur...@gmail.com wrote:
Hello,
I am struggling with computing nodes in Unix.
I have the use of a Unix server that has 30 nodes and I would like to
batch scripts.
Here is an R example that results in 72 repeated tasks based on the 2
loops. If I wanted to send these out to the different nodes, each node
has 1 task and the remaining 42 tasks are put in a queue, what would I
do?
#Example:
proj.dir-C:/test/
setwd(proj.dir)
group.vars -
c(biol1,biol4,biol5,biol6,biol12,biol15,biol16,biol17)
species - c(aa, bb, cc, dd, ee, ff, gg, hh, ii)
for (gv in 1:length(group.vars)){
for (sp in 1:length(species)){
zzz - file(paste(group.vars[gv],species[sp],.TXT,sep=),w)
cat(paste(group.vars[gv]),file=zzz)
cat(paste(species[sp]),file=zzz)
close(zzz)
}
}
I know it is a silly script test but the division is similar to my
real task which is repeated 459 times and each run takes 18 hours
At this point I am confused what should be written in Unix and then at
what point should I call R. I have read an abundance of things but I
feel like I am missing something essential. Or perhaps I have read all
the wrong things.
Thanks,
Daisy
--
Daisy Englert Duursma
Room E8C156
Dept. Biological Sciences
Macquarie University NSW 2109
Australia
__
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.
--
Stephen Sefick
| Auburn University |
| Department of Biological Sciences |
| 331 Funchess Hall |
| Auburn, Alabama |
| 36849 |
|___|
| sas0...@auburn.edu |
| http://www.auburn.edu/~sas0025 |
|___|
Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods. We are mammals, and have not exhausted the
annoying little problems of being mammals.
-K. Mullis
A big computer, a complex algorithm and a long time does not equal science.
-Robert Gentleman
__
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] Regular expressions: offsets of groups
Ok, we decided to have a shot at modifying gregexpr. Let's see how it works out. If anybody is interested in discussing this please contact me. R-help doesn't seem like the right place for further discussion. Is there a default place for discussing things like that? Thanks everybody for your responses! Titus __ 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] conditional assignment of colors in xyplot()
On Mon, Sep 27, 2010 at 7:11 AM, nathan pellegrin
nathan.pelleg...@gmail.com wrote:
# Dear R Community,
# I have this data frame:
df1 - data.frame(
F1 = factor( c( rep(D1,12),rep(D2,12),rep(D3,12) ) ),
F2 = factor( rep( rep( paste(O,1:6,sep=), rep(2,6) ), 3) ),
F3 = factor( rep( c(V1,V2), 18 ) ),
S1 =
c(8.840955e-02,2.546822e-01,7.562658e-01,8.801181e-01,6.041766e-02,2.172731e-01,6.542187e-98,
7.067840e-04,1.430933e-01,9.764401e-01,1.380848e-05,1.192620e-01,9.107259e-01,1.235232e-01,
3.021847e-01,1.331351e-01,5.272103e-01,3.663577e-01,1.539690e-38,2.224451e-01,2.873052e-01,
5.110313e-01,7.840802e-01,8.210762e-10,1.553356e-01,4.173335e-01,5.964021e-01,4.955694e-01,
8.849240e-02,5.739598e-01,1.879075e-17,1.071003e-03,7.298928e-01,6.347287e-01,8.884034e-01,
4.460295e-11),
S2 =
c(1.32249139,1.02831831,-0.09650252,-0.05454486,2.62105492,2.00310250,8.07269064,
-1.55397883,1.77390551,0.04161954,7.14188540,-2.98033713,-0.49904251,-0.74309058,
-0.49904251,-0.74309058,1.22492623,-1.79003492,7.60003121,-0.74549596,2.53418936,
-1.60112296,0.67131380,-15.31744351,-0.18380339,0.28628435,-0.18380339,0.28628435,
2.96108998,1.18267783,5.78419118,2.70861763,0.66287857,1.10397741,0.27160971,
-15.37506924) )
# Two of the factors are used to define an array
# of panels with the third showing up as bars along the x axis.
# How do I color the bars according to their sign (red 0, blue 0)
# conditional on the value of S2 .025 - in which case color them gray?
# Initially, I tried to pass a character vector specifying colors,
# which does not achieve the desired result:
library(lattice)
library(latticeExtra)
df1$barcols - ifelse(df1$S1 .025, gray, ifelse( df1$S2 0,
blue,red))
ctp - xyplot(S2 ~ F2 | F1 + F3,
data=df1, as.table=TRUE, ylim=c(-10,10),
panel = function(x,y){
panel.barchart(x,y,horizontal=FALSE, origin=0,
col=df1$barcols ) } )
useOuterStrips(ctp)
# Any help you can provide would be greatly appreciated.
# Thank you!
Itis better to keep the idea of grouping and the actual colors distinct.
df1 -
transform(df1,
colgrp = factor(ifelse(S1 .025, A,
ifelse( df1$S2 0, B, C)),
levels = c(A, B, C)))
barchart(S2 ~ F2 | F1 + F3, stack = TRUE,
data=df1, as.table=TRUE, ylim=c(-10,10),
groups = colgrp,
par.settings = simpleTheme(col = c(grey, blue, red)))
-Deepayan
__
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] how to make list() return a list of *named* elements
If I combine elements into a list b - c(22.4, 12.2, 10.9, 8.5, 9.2) my.c - sample.int(round(2*mean(b)), 5) my.list - list(b, my.c) the names of the elements seems to get lost in the process: str(my.list) List of 2 $ : num [1:5] 22.4 12.2 10.9 8.5 9.2 $ : int [1:5] 11 8 6 9 20 If I explicitly name the elements at list-creation, I get what I want: my.list - list(b=b, my.c=my.c) str(my.list) List of 2 $ b : num [1:5] 22.4 12.2 10.9 8.5 9.2 $ my.c: int [1:5] 11 8 6 9 20 Now, is there a way to get list() (or some other function) to automatically name the elements? I often use list() in return(), and I am getting tired of having to repeat myself. -- Hans Ekbrand (http://sociologi.cjb.net) h...@sociologi.cjb.net A. Because it breaks the logical sequence of discussion Q. Why is top posting bad? signature.asc Description: Digital signature __ 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] cor() alternative for huge data set
Peter, Many thank for suggesting me this package. I very much believe that this will help me. But I was trying to correlate all probes(correlation between entities not variables) to calculate differentially coexpressed gene sets using package coXpress in R. I could not reduce the number on the basis of intensity, since most of the genes are down regulated and upregulated in treated conditions, so they are of my interest and cannot be removed from control samples(since I have to compare both). can you further suggest me an alternative for differentially coexpression analysis, since this is what I need to know the most-- the sets which are behaving differently across conditions. Has any one ever used this package--coXpress?? Regards .. Jyotasana - Original Message - From: Peter Langfelder peter.langfel...@gmail.com To: Jyotasana Gulati jgul...@ice.mpg.de Cc: r-help@r-project.org Sent: Thursday, September 30, 2010 4:05:44 AM Subject: Re: [R] cor() alternative for huge data set On Wed, Sep 29, 2010 at 1:27 PM, Jyotasana Gulati jgul...@ice.mpg.de wrote: Hi, I am have a data set of around 43000 probes(rows), and have to calculate correlation matrix. When I run cor function in R, its throwing an error message of RAM shortage which was obvious for such huge number of rows. I am not getting a logical way to cut off this huge number of entities, is there an alternative to pearson correlation or with other dist() methods calculation(euclidean) that can be run on such a huge data set?? Every help will be appreciated. Hmm... Are you calculating a correlation of 43000 probes, or of some number of samples across 43000 probes? If the former, read below. If the latter, I'm surprised you are running out of memory. Issuing garbage collection (gc()) before the calculation, closing all other programs, removing all other large objects from the R workspace etc. may help. If you really need the 43k times 43k correlation matrix of your 43k probes, read on. [Disclosure: this is a shameless plug for the package WGCNA (Weighted Gene Co-expression Network Analysis, also known as Weighted Correlation Network Analysis), from the package author, namely me.] First, since the distance matrix will be huge, you will not gain using other distance methods either. Second, depending on what you want to do with the 43k probes, the package WGCNA may help you. It has methods for creating correlation networks among a large number of probes. The idea is to pre-cluster the probes using what I call projective K-means, function projectiveKMeans. The pre-clustering will return what we call blocks of probes (or genes). We assume (this is a big assumption) that correlations among probes belonging to different blocks can be neglected. Then we treat each block separately for network construction (or, in your case, possibly simple calculation of correlation). Although this isn't strictly an R topic but rather microarray analysis issue, in my experience it is often useful to filter out probes before actually calculating and interpreting large correlation matrices. In conjunction with filtering, it can be advantageous to only keep one probe per gene (presumably there is more than one probe per gene in you data set). The filtering criterion varies from analysis to analysis, but if your data represent intensities, it is often a good idea to throw away probes whose intensity is always low, because such signals are mostly noise. If you decide to check out WGCNA, look at http://www.genetics.ucla.edu/labs/horvath/CoexpressionNetwork/Rpackages/WGCNA/. Peter __ 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] how to make list() return a list of *named* elements
You should try:
eapply(.GlobalEnv, I)[c('b', 'my.c')]
On Thu, Sep 30, 2010 at 8:49 AM, Hans Ekbrand
hans.ekbr...@sociology.gu.sewrote:
If I combine elements into a list
b - c(22.4, 12.2, 10.9, 8.5, 9.2)
my.c - sample.int(round(2*mean(b)), 5)
my.list - list(b, my.c)
the names of the elements seems to get lost in the process:
str(my.list)
List of 2
$ : num [1:5] 22.4 12.2 10.9 8.5 9.2
$ : int [1:5] 11 8 6 9 20
If I explicitly name the elements at list-creation, I get what I want:
my.list - list(b=b, my.c=my.c)
str(my.list)
List of 2
$ b : num [1:5] 22.4 12.2 10.9 8.5 9.2
$ my.c: int [1:5] 11 8 6 9 20
Now, is there a way to get list() (or some other function) to
automatically name the elements?
I often use list() in return(), and I am getting tired of having to
repeat myself.
--
Hans Ekbrand (http://sociologi.cjb.net) h...@sociologi.cjb.net
A. Because it breaks the logical sequence of discussion
Q. Why is top posting bad?
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.9 (GNU/Linux)
iEYEARECAAYFAkykeUwACgkQfCyHKnBQYU4J+ACgrdjMSoyr/Uzt9fpTsietde3n
d8UAnRskbOM7mDhDexiS7T9LkhRs287P
=MsDG
-END PGP SIGNATURE-
__
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.
--
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40 S 49° 16' 22 O
[[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] interactive session
Hi Fatih
I believe that readLines(n=1) will do the job. It works
fine from the Windows RGui, but I noticed that it hangs my Aquamacs/ESS
when R runs from there, and a C-g was needed (which may be completely
irrelevant to you).
Best, Niels
On 30/09/10 08.55, Pam wrote:
Hi guys,
My concern is to create an automated process from the beginning to the end. I
want to copy all my code together in one piece and paste it to R console and sit
back and relax :) except one moment in which the program should ask me to enter
a number.. and only then (only after getting a valid number from me) it should
continue to read and process the rest of the code. I�tried�lots of things
(readline, readLines, scan, interactive, ask, switch,...) and read manuals and
searched help forums.. I found several similar questions but never a satisfying
answer.. so now it became a challenge.. any idea how�to�overcome this problem (R
2.11.1 for Windows)? (an example�code is below)�
Best,
Fatih
library(gtools)
library(YaleToolkit)
library(xts)
�
### start of my�wrong�function!
f-function(w){
�w-readline(which data? )
�w-as.numeric(w)
�ifelse(is.numeric(w)==TRUE, w, f())
�}
f()
# end of my�wrong function
v- ## and output of my function should be a v�for example which I can use it
in the next line�(v-w� or something like that??)
�
##the rest works fine
p-paste(t, v, .txt, sep = )
t-read.table(p, header=FALSE, sep=\t, dec=,,
blank.lines.skip=FALSE)
rownames(t)-as.Date(t[,1],%d.%m.%Y)
colnames(t)-c(date,start,high,low,end,w.average,lot,
volume)
x-as.xts(t)
whatis(x)���
.
.
[[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.
--
Niels Richard Hansen Web: www.math.ku.dk/~richard
Associate Professor Email: niels.r.han...@math.ku.dk
Department of Mathematical Sciences nielsrichardhan...@gmail.com
University of Copenhagen Skype: nielsrichardhansen.dk
Universitetsparken 5 Phone: +45 353 20783 (office)
2100 Copenhagen Ø +45 2859 0765 (mobile)
Denmark
__
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] Opening a .R file with R (Windows)
Here's a thread that has some good discussion on the topic
http://r.789695.n4.nabble.com/Running-script-with-double-click-td1579014.html
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Bert Gunter
Sent: Tuesday, September 28, 2010 1:48 PM
To: Joshua Wiley
Cc: r-help@r-project.org
Subject: Re: [R] Opening a .R file with R (Windows)
Well, try following the correct conventions ...
1. Double click on an .Rdata file, which is produced by saving from R, and it
will open.
2. Drag and drop a .R or any text file into an open R window and it will source
the contents.
This is probably documented somewhere .. maybe in the RW FAQ.
-- Bert
On Tue, Sep 28, 2010 at 10:22 AM, Joshua Wiley jwiley.ps...@gmail.com wrote:
Hi Kye,
I have never gotten .R files to work quite like other types (e.g.,
double-clicking a .PDF) in Windows. AFAIK there is no simple way to
do it, because you do not edit scripts directly in R (I am happy to be
corrected if someone knows better). For general use, I would just
open R first and then open the file, or if you just want to run the
file, you can use R's batch mode from the Windows command prompt.
Best regards,
Josh
On Tue, Sep 28, 2010 at 10:11 AM, Kye Gilder kye.gil...@gmail.com wrote:
I am new to using R. I installed R on my computer (Windows) and
everything things appears to be just fine. However, I have a simple
script RTest.R that does a few simple calculations. When I
double-click the RTest.R icon, I get an Information dialong box which
says, ARGUMENT 'C:\Documents and Settings\kgilder\Desktop\RTest.R'
__ignored__ . When I choose OK, R then opens, but it does not open or
display the script. Any help or suggestions?
When I open R and use File Open New Script and path to the file,
it opens just fine.
Regards,
Kye
[[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.
--
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.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.
--
Bert Gunter
Genentech Nonclinical Biostatistics
__
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.
===
P Please consider the environment before printing this e-mail
Cleveland Clinic is ranked one of the top hospitals
in America by U.S.News World Report (2009).
Visit us online at http://www.clevelandclinic.org for
a complete listing of our services, staff and
locations.
Confidentiality Note: This message is intended for use\...{{dropped:13}}
__
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] How to get a proportion of a Vector Member
Hi or prop.table(table(foo)) Regards Petr r-help-boun...@r-project.org napsal dne 30.09.2010 03:45:30: On Wed, Sep 29, 2010 at 6:43 PM, Gundala Viswanath gunda...@gmail.com wrote: I have a vector that looks like this: foo [1] o o o x o o o o o x x o x How can we find the percentage of o and x in that vector in R? table(foo)/length(foo) Peter __ 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] how to make list() return a list of *named* elements
On Thu, Sep 30, 2010 at 7:49 AM, Hans Ekbrand hans.ekbr...@sociology.gu.se wrote: If I combine elements into a list b - c(22.4, 12.2, 10.9, 8.5, 9.2) my.c - sample.int(round(2*mean(b)), 5) my.list - list(b, my.c) the names of the elements seems to get lost in the process: str(my.list) List of 2 $ : num [1:5] 22.4 12.2 10.9 8.5 9.2 $ : int [1:5] 11 8 6 9 20 If I explicitly name the elements at list-creation, I get what I want: my.list - list(b=b, my.c=my.c) str(my.list) List of 2 $ b : num [1:5] 22.4 12.2 10.9 8.5 9.2 $ my.c: int [1:5] 11 8 6 9 20 Now, is there a way to get list() (or some other function) to automatically name the elements? I often use list() in return(), and I am getting tired of having to repeat myself. A data frame is a list in which every component (i.e. every column) must have the same length (i.e. the same number of rows). data.frame() does preserve names: data.frame(b, my.c) b my.c 1 22.48 2 12.29 3 10.9 15 4 8.51 5 9.2 14 -- Statistics Software Consulting GKX Group, GKX Associates Inc. tel: 1-877-GKX-GROUP email: ggrothendieck at 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] What's the meaning of Species ~ . in IRIS data
See ?lm and, more usefully, ?formula
Gundala Viswanath gunda...@gmail.com 29/09/2010 15:51
I am refering to a function call like this:
data(iris)
x - svmlight(Species ~ ., data = iris)
I tried to see the content of it by typing:
Species ~ .
but it gives nothing. How can I see it's content ?
- P.Dubois
__
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.
***
This email and any attachments are confidential. Any use...{{dropped:8}}
__
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] String split and concatenation
Steven,
I noticed that I hit the send button too early. I answered the question you
asked, but probably not the one that you should be asking. Why do you want to
do this? It looks suspiciously like you are trying to create code to then
evaluate. That is usually like printing a document, scanning it back in, an
importing it as a graphic in another document so they can change the size (yes
this happens), instead of just changing the font size and margins.
If you tell us more of your ultimate goal, we can probably be more helpful.
--
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-boun...@r-
project.org] On Behalf Of Greg Snow
Sent: Wednesday, September 29, 2010 12:32 PM
To: Steven Kang; bill.venab...@csiro.au
Cc: r-help@r-project.org
Subject: Re: [R] String split and concatenation
paste( '(', paste( ', rep(letters[1:3],2), ', sep=,
collapse=','), ')', sep= )
[1] ('a','b','c','a','b','c')
If you need the space after the comma then just change ',' to ', '.
The outer paste can be replaced with sprintf (and that may be more
readable).
--
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-boun...@r-
project.org] On Behalf Of Steven Kang
Sent: Wednesday, September 29, 2010 2:16 AM
To: bill.venab...@csiro.au
Cc: r-help@r-project.org
Subject: Re: [R] String split and concatenation
x - rep(letters[1:3], 2)
Are there any ways to transform assign the above as the one shown
below
to an object? (in exact format; i.e length of 1 class of
character),
i.e
x
('a', 'b', 'c', 'a', 'b', 'c')
Highly appreciate for any advice.
On Wed, Sep 29, 2010 at 3:33 PM, bill.venab...@csiro.au wrote:
dump(x, file = x.R)
file.show(x.R)
will get you most of the way.
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
project.org]
On Behalf Of Steven Kang
Sent: Wednesday, 29 September 2010 3:11 PM
To: r-help@r-project.org
Subject: [R] String split and concatenation
Hi R users,
I desire to transform the following vector consisting of repeated
characters
x - rep(letters, 3)
into this exact format (i.e a single string containing each
characters in
quotation mark separated by comma between each; al ).
(a, b, c, d, a, b, c, d, ..., a,
b,
c,
d, .z)
Any advice would be much appreciated.
--
Steven
[[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.htmlhttp://www.r-
project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Steven
[[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-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] Unix batch to different nodes
stephen sefick ssefick at gmail.com writes: This, I think, all depends on how the server is set up. On Thu, Sep 30, 2010 at 1:21 AM, Daisy Englert Duursma daisy.duursma at gmail.com wrote: Hello, I am struggling with computing nodes in Unix. I have the use of a Unix server that has 30 nodes and I would like to batch scripts. [snip] Do you mean a cluster with 30 nodes? Or a single workstation with 30 CPUs/cores? I don't know the exact answer (I'm trying to figure it out for myself on my local cluster), but see the High Performance Computing task view: depending on what software you have or can get installed on the cluster, some combination of the foreach, snow, Rmpi ... packages will be helpful). __ 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] Opening a .R file with R (Windows)
On Tue, Sep 28, 2010 at 1:11 PM, Kye Gilder kye.gil...@gmail.com wrote: I am new to using R. I installed R on my computer (Windows) and everything things appears to be just fine. However, I have a simple script RTest.R that does a few simple calculations. When I double-click the RTest.R icon, I get an Information dialong box which says, ARGUMENT 'C:\Documents and Settings\kgilder\Desktop\RTest.R' __ignored__ . When I choose OK, R then opens, but it does not open or display the script. Any help or suggestions? When I open R and use File Open New Script and path to the file, it opens just fine. There is a batch file called #Rscript.bat at http://batchfiles.googlecode.com . If you (1) place #Rscript.bat anywhere on your path and then (2) place this line at the top of your R file: #Rscript %0 %* and (3) rename your R file to end in .bat then it can be used as both a bat file and as an R file. If you double click it then it will run the script since its a bat file yet you can also source() it into R like any R file.It does produce one garbage line of output at the top so it will only be usable in situations where that is ok. Note that #Rscript.bat finds R from the registry so if you upgrade R you will not need to make any changes to your scripts that use it. -- Statistics Software Consulting GKX Group, GKX Associates Inc. tel: 1-877-GKX-GROUP email: ggrothendieck at 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] R crashes when loading rgl package before minqa package
Gaspard Lequeux Gaspard.Lequeux at biomath.ugent.be writes:
Calling newuoa (from the minqa package) makes R crash when
the package rgl
is loaded first. This however only on certain selected data.
The data used for testing (saved to 'bugs.R'):
xvals =
c(1,2,4,5,7,8,9,10,11,12,14,15,16,18,19,21,22,23,24,25,26,
27,28,29,30,31,32,33,34,35,36)
yvals =
c(857.7597,975.8624,978.2655,979.3034,965.5919,983.8946,992.2512,
992.1178,979.5379,974.4269,968.4113,991.5210,977.3361,985.7800,
975.5220,974.6880,973.8102,980.7295,982.0034,984.7993,978.4948,
970.4351,969.0718,983.7892,976.3637,980.7833,987.1665,976.6000,
975.1332,971.0757,989.4693)
initpar = c(-5.1471384, -3861.8905839, 979.2616002, 0.2572355, 27.5705764)
optimft - function(x) {
yft = x[2] + (x[3] - x[2])/((1 + exp(x[1] * (log(xvals) -
log(x[4]^x[5])
return(sum((yvals - yft)^2))
}
Sequence of commands needed to make the bug appear:
Start R
source('bugs.R')
library(minqa)
library(rgl)
newuoa(initpar, optimft)
= OK
Start R
source('bugs.R')
library(rgl)
library(minqa)
newuoa(initpar, optimft)
= Crash: segfault: address 0x18, cause 'memory not mapped'
I found the bug using the package qpcR, where rgl is loaded when loading
qpcR while minqa is only loaded later, when needed.
Running on Debian squeeze 64 bit.
R version: R version 2.11.1 (2010-05-31) x86_64-pc-linux-gnu
rgl version: 0.91
minqa version: 1.1.9
Rcpp version: 0.8.6 (loaded by minqa)
Kind regards,
Gaspard Lequeux
Duplicated on Ubuntu 10.04
sessionInfo()
R version 2.11.1 (2010-05-31)
i486-pc-linux-gnu
locale:
[1] LC_CTYPE=en_CA.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_CA.UTF-8LC_COLLATE=en_CA.UTF-8
[5] LC_MONETARY=C LC_MESSAGES=en_CA.UTF-8
[7] LC_PAPER=en_CA.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] minqa_1.1.9 Rcpp_0.8.6 rgl_0.91.787
Running valgrind:
bol...@ubuntu-10:~/R/misc$ R -d valgrind --vanilla
==26985== Memcheck, a memory error detector
==26985== Copyright (C) 2002-2009, and GNU GPL'd, by Julian Seward et al.
==26985== Using Valgrind-3.6.0.SVN-Debian and LibVEX; rerun with -h for
copyright info
==26985== Command: /usr/lib/R/bin/exec/R --vanilla
==26985==
R version 2.11.1 (2010-05-31)
[snip startup info]
xvals =
c(1,2,4,5,7,8,9,10,11,12,14,15,16,18,19,21,22,23,
24,25,26,27,28,29,30,31,32,33,34,35,36)
yvals =
c(857.7597,975.8624,978.2655,979.3034,965.5919,
983.8946,992.2512,992.1178,979.5379,974.4269,968.4113,
991.5210,977.3361,985.7800,975.5220,974.6880,973.8102,
980.7295,982.0034,984.7993,978.4948,970.4351,969.0718,
983.7892,976.3637,980.7833,987.1665,976.6000,975.1332,
971.0757,989.4693)
initpar = c(-5.1471384, -3861.8905839, 979.2616002, 0.2572355, 27.5705764)
optimft - function(x) {
+yft = x[2] + (x[3] - x[2])/((1 + exp(x[1] * (log(xvals) -
log(x[4]^x[5])
+return(sum((yvals - yft)^2))
+ }
library(rgl)
library(minqa)
Loading required package: Rcpp
newuoa(initpar, optimft)
==26985== Invalid read of size 4
==26985==at 0x6763C17: __cxa_allocate_exception (in
/usr/lib/libstdc++.so.6.0.13)
==26985==by 0x6EC026E: calfun_ (minqa.cpp:30)
==26985==by 0x6EC4F56: newuob_ (newuob.f:323)
==26985==by 0x6EC4ADA: newuoa_ (newuoa.f:68)
==26985==by 0x6EC2350: newuoa_cpp__rcpp__wrapper__(Rcpp::Vector14,
Rcpp::Environment, SEXPREC*) (minqa.cpp:120)
==26985==by 0x6EC27D1: newuoa_cpp (minqa.cpp:110)
==26985==by 0x40BB170: ??? (in /usr/lib/R/lib/libR.so)
==26985==by 0x40F20C1: Rf_eval (in /usr/lib/R/lib/libR.so)
==26985==by 0x40F42BF: ??? (in /usr/lib/R/lib/libR.so)
==26985==by 0x40F1E37: Rf_eval (in /usr/lib/R/lib/libR.so)
==26985==by 0x40F5C6F: Rf_applyClosure (in /usr/lib/R/lib/libR.so)
==26985==by 0x40F1CDC: Rf_eval (in /usr/lib/R/lib/libR.so)
==26985== Address 0xc is not stack'd, malloc'd or (recently) free'd
==26985==
*** caught segfault ***
address 0xc, cause 'memory not mapped'
Traceback:
1: .Call(newuoa_cpp, par, ctrl, fn1)
2: newuoa(initpar, optimft)
Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace
Selection:
I strongly suspect that the problem is with minqa, and
that loading rgl is just a way to make the problem surface.
I may take a look at minqa , but you may want to send an e-mail
to the maintainer [ maintainer(minqa) ] as well ...
Ben Bolker
__
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] AIC for tweedie glm
Dear R-users, I'm trying to model some data using a tweedie GLM approach. My response variable is the number of pupae that are the offspring of a subordinate wasp on a wasp's nest. However, they're not count data- for each nest, I only know the mean number of pupae per subordinate, which is continous. The data also contain a high proportion of zeros. I'm not very experienced at statistical modelling, but from reading previous posts, it seems that my data would suit a tweedie approach. I can't use a zero-inflated Poisson model, because my data are not counts. Many of my values are between 0 and 1, so if I rounded to the nearest integer, I'd lose a lot of the variation. Here's my code: out-tweedie.profile(PUPAE_PER_SUB~1,p.vec=seq(1.1,1.9,length=9),method=interpolation,do.ci=TRUE,do.smooth=TRUE,do.plot=TRUE) tweedie1-glm(GSA_TOTAL_DF_PERSUB~GROUP_SIZE+PERIOD+SITE+PERIOD*GROUP_SIZE,family=tweedie(var.power=out$p.max,link.power=0)) This worked fine, and gave results I expected, but I don't know what the best method is to evaluate the fit of the model. I am used to using AIC to compare models. A site search turned up AICtweedie, within the tweedie package, but I get the following message: Error: could not find function AICtweedie when I try to use this command, even though tweedie and statmod are both loaded. I've also read that AIC can be calculated using dtweedie, but I'm a beginner and so, despite lots of searching, I'm not sure how. I'm sorry to ask a basic statistics rather than programming question, but I'm really stuck. Could anyone advise me on the best way to assess goodness-of-fit for this type of model, in order to compare models? Thanks -- View this message in context: http://r.789695.n4.nabble.com/AIC-for-tweedie-glm-tp2720813p2720813.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] conditional assignment of colors in xyplot()
Thank you, Deepayan, for the answer - and for creating amazingly helpful
tools!
Thanks also to Jim Lemon who pointed out a typo in the code: S2
.025 should read S2 .025,
Nathan
On Thu, Sep 30, 2010 at 5:21 AM, Deepayan Sarkar
deepayan.sar...@gmail.comwrote:
On Mon, Sep 27, 2010 at 7:11 AM, nathan pellegrin
nathan.pelleg...@gmail.com wrote:
# Dear R Community,
# I have this data frame:
df1 - data.frame(
F1 = factor( c( rep(D1,12),rep(D2,12),rep(D3,12) ) ),
F2 = factor( rep( rep( paste(O,1:6,sep=), rep(2,6) ), 3) ),
F3 = factor( rep( c(V1,V2), 18 ) ),
S1 =
c(8.840955e-02,2.546822e-01,7.562658e-01,8.801181e-01,6.041766e-02,2.172731e-01,6.542187e-98,
7.067840e-04,1.430933e-01,9.764401e-01,1.380848e-05,1.192620e-01,9.107259e-01,1.235232e-01,
3.021847e-01,1.331351e-01,5.272103e-01,3.663577e-01,1.539690e-38,2.224451e-01,2.873052e-01,
5.110313e-01,7.840802e-01,8.210762e-10,1.553356e-01,4.173335e-01,5.964021e-01,4.955694e-01,
8.849240e-02,5.739598e-01,1.879075e-17,1.071003e-03,7.298928e-01,6.347287e-01,8.884034e-01,
4.460295e-11),
S2 =
c(1.32249139,1.02831831,-0.09650252,-0.05454486,2.62105492,2.00310250,8.07269064,
-1.55397883,1.77390551,0.04161954,7.14188540,-2.98033713,-0.49904251,-0.74309058,
-0.49904251,-0.74309058,1.22492623,-1.79003492,7.60003121,-0.74549596,2.53418936,
-1.60112296,0.67131380,-15.31744351,-0.18380339,0.28628435,-0.18380339,0.28628435,
2.96108998,1.18267783,5.78419118,2.70861763,0.66287857,1.10397741,0.27160971,
-15.37506924) )
# Two of the factors are used to define an array
# of panels with the third showing up as bars along the x axis.
# How do I color the bars according to their sign (red 0, blue 0)
# conditional on the value of S2 .025 - in which case color them gray?
# Initially, I tried to pass a character vector specifying colors,
# which does not achieve the desired result:
library(lattice)
library(latticeExtra)
df1$barcols - ifelse(df1$S1 .025, gray, ifelse( df1$S2 0,
blue,red))
ctp - xyplot(S2 ~ F2 | F1 + F3,
data=df1, as.table=TRUE, ylim=c(-10,10),
panel = function(x,y){
panel.barchart(x,y,horizontal=FALSE, origin=0,
col=df1$barcols ) } )
useOuterStrips(ctp)
# Any help you can provide would be greatly appreciated.
# Thank you!
Itis better to keep the idea of grouping and the actual colors distinct.
df1 -
transform(df1,
colgrp = factor(ifelse(S1 .025, A,
ifelse( df1$S2 0, B, C)),
levels = c(A, B, C)))
barchart(S2 ~ F2 | F1 + F3, stack = TRUE,
data=df1, as.table=TRUE, ylim=c(-10,10),
groups = colgrp,
par.settings = simpleTheme(col = c(grey, blue, red)))
-Deepayan
[[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] R crashes when loading rgl package before minqa package
I get this on Windows (it does not crash):
library(minqa)
library(rgl)
newuoa(initpar, optimft)
Error in newuoa(initpar, optimft) :
non-finite x values not allowed in calfun
In addition: Warning message:
In log(x[4]) : NaNs produced
This tells me that you should be constraining your parameter x[4] (may be even
x[5]) to be non-negative:
Here is what I get with `bobyqa':
bobyqa(initpar, optimft, lower=c(-Inf, -Inf, -Inf, 0, 0))
parameter estimates: -5.311767080681, -3861.89005072333, 979.239647766226,
0.268156271922112, 27.6418856936228
objective: 1457.20987728737
number of function evaluations: 78
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Gaspard Lequeux gaspard.lequ...@biomath.ugent.be
Date: Wednesday, September 29, 2010 11:40 am
Subject: [R] R crashes when loading rgl package before minqa package
To: r-help@r-project.org
Hej,
Calling newuoa (from the minqa package) makes R crash when the
package rgl is loaded first. This however only on certain selected data.
The data used for testing (saved to 'bugs.R'):
xvals =
c(1,2,4,5,7,8,9,10,11,12,14,15,16,18,19,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36)
yvals =
c(857.7597,975.8624,978.2655,979.3034,965.5919,983.8946,992.2512,992.1178,979.5379,974.4269,968.4113,991.5210,977.3361,985.7800,975.5220,974.6880,973.8102,980.7295,982.0034,984.7993,978.4948,970.4351,969.0718,983.7892,976.3637,980.7833,987.1665,976.6000,975.1332,971.0757,989.4693)
initpar = c(-5.1471384, -3861.8905839, 979.2616002, 0.2572355, 27.5705764)
optimft - function(x) {
yft = x[2] + (x[3] - x[2])/((1 + exp(x[1] * (log(xvals) -
log(x[4]^x[5])
return(sum((yvals - yft)^2))
}
Sequence of commands needed to make the bug appear:
Start R
source('bugs.R')
library(minqa)
library(rgl)
newuoa(initpar, optimft)
= OK
Start R
source('bugs.R')
library(rgl)
library(minqa)
newuoa(initpar, optimft)
= Crash: segfault: address 0x18, cause 'memory not mapped'
I found the bug using the package qpcR, where rgl is loaded when
loading qpcR while minqa is only loaded later, when needed.
Running on Debian squeeze 64 bit.
R version: R version 2.11.1 (2010-05-31) x86_64-pc-linux-gnu
rgl version: 0.91
minqa version: 1.1.9
Rcpp version: 0.8.6 (loaded by minqa)
Kind regards,
Gaspard Lequeux
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
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] Ordered logit with polr won't match SPSS output
I think the most common reason to see different parameter estimates with ordinal regression is that programs set up the model differently. For example, check out library(MASS) ?polr We see polr uses: logit P(Y = k | x) = zeta_k - eta and notes that other software packages may use the opposite sign for eta. Also check out that the ordered variables have the same order (reference category). I sometimes have messed that up. On Mon, Sep 27, 2010 at 2:53 PM, Ben Hunter bjameshun...@gmail.com wrote: I am learning R via a textbook that performs analysis with SPSS and SAS. In trying to reproduce the results for an ordinal logit model, I get very similar point estimates for my cut-off points, but the parameters for the covariate q60 do not match. The estimate for q51 also matches. Is this because I need to change a base case for the ordered covariate q60? Can this be done in or is it always the first case? mod-polr(as.ordered(q43j)~as.ordered(q60)+q51, method=logistic) Perhaps a book using R would be a better idea, but it's the content and price (free) of this book that I like. Thanks a lot, Ben [[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] AIC for tweedie glm
eleadbeater e.leadbeater at sussex.ac.uk writes: Dear R-users, I'm trying to model some data using a tweedie GLM approach. My response variable is the number of pupae that are the offspring of a subordinate wasp on a wasp's nest. However, they're not count data- for each nest, I only know the mean number of pupae per subordinate, which is continous. The data also contain a high proportion of zeros. This worked fine, and gave results I expected, but I don't know what the best method is to evaluate the fit of the model. I am used to using AIC to compare models. A site search turned up AICtweedie, within the tweedie package, but I get the following message: Error: could not find function AICtweedie when I try to use this command, even though tweedie and statmod are both loaded. I've also read that AIC can be calculated using dtweedie, but I'm a beginner and so, despite lots of searching, I'm not sure how. I'm sorry to ask a basic statistics rather than programming question, but I'm really stuck. Could anyone advise me on the best way to assess goodness-of-fit for this type of model, in order to compare models? Everything you're saying sounds sensible. The only (!) problem is that your problem isn't reproducible (for me at least). What happens if you run library(tweedie) example(AICtweedie) from a fresh R session (possibly using R --vanilla)? What are the results of sessionInfo() ? Works for me. Ben Bolker __ 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] polr, lrm - ordinal data
Dear List, I have developed a model and am looking to predict a response for 1-6 ( it is ordered i.e the difference between level 1 and 2 is the same as between level 2 and 3 etc. I have used the predict function for a polr model (below) and a lrm model, and both give similar results, however for some reason the outcome are all 1 or 6 i.e no level 2,3,4,5. This is not correct, i am unsure if it is because my model is not good enough to predict to this accuracy or if it is something i am doing wrong? Thanks Chris test - polr(extinction ~ FR*HAB+WO+ALT+BIO+REG,method=logistic) test Call: polr(formula = extinction ~ FR * HAB + WO + ALT + BIO + REG, method = logistic) Coefficients: FRNon_fleshy HABSemi-aquatic HABTerrestrial WOWoody ALTHigh 0.09758543 -0.05988101 -0.29744997 0.32746840 -0.39191606 ALTLow ALTMid BIOBoreal BIOMediterranean BIOSubtropical/Tropical -0.56523156 -0.18979562 -0.22743656 -0.31344233 1.77031824 BIOTemperate REGTwo_plus FRNon_fleshy:HABSemi-aquatic FRNon_fleshy:HABTerrestrial 1.45071627 -0.67654880 -2.06919408 -0.31706541 Intercepts: 1|2 2|3 3|4 4|5 5|6 0.4874828 0.8340901 1.3994091 1.8091463 2.1295630 Residual Deviance: 2471.053 AIC: 2509.053 predict(test) [1] 1 1 1 1 1 1 1 1 1 1 6 1 6 1 1 6 1 1 1 6 1 1 1 6 1 1 1 1 6 1 6 6 1 1 6 1 1 1 1 6 1 1 1 1 6 6 1 1 6 1 1 1 1 1 6 1 1 1 1 6 1 1 6 1 6 1 6 1 1 1 1 1 1 1 [75] 1 1 1 1 1 1 1 6 6 1 1 1 6 6 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 6 6 1 1 1 1 1 1 6 1 1 1 1 1 6 1 1 1 1 1 6 6 1 6 1 1 1 1 1 1 6 1 1 1 1 1 6 1 1 1 1 6 [149] 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 6 1 1 1 1 1 1 6 1 6 1 1 1 1 1 6 1 1 6 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 6 1 1 6 1 1 1 1 1 1 6 1 1 [223] 1 1 1 1 1 1 1 1 1 1 1 6 1 1 6 6 6 1 1 6 6 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 6 1 1 1 1 1 1 6 1 1 1 1 6 1 1 6 6 1 6 1 1 1 1 1 1 1 1 1 1 1 1 6 1 6 6 1 1 [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 6 1 1 6 1 1 6 1 6 1 1 6 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 6 1 1 [371] 1 1 1 1 1 1 1 6 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 6 6 6 1 1 1 1 6 1 1 1 1 1 1 1 1 6 1 1 1 6 1 1 1 1 1 1 1 6 1 1 1 1 6 1 6 1 1 6 1 1 1 1 1 1 1 1 1 1 1 [445] 1 6 1 1 1 1 1 1 6 1 1 1 1 1 6 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 [519] 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 6 1 1 1 1 1 1 1 6 1 1 6 1 6 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 6 1 1 6 1 6 1 1 6 6 1 1 1 1 1 1 1 1 1 1 1 [593] 6 1 1 1 1 1 1 6 1 1 6 1 1 6 6 1 6 1 1 1 6 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 6 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 [667] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 [741] 1 1 1 6 6 1 1 1 6 1 1 1 1 1 6 1 1 1 1 1 1 1 6 1 1 6 1 1 6 1 1 1 6 1 6 1 6 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 6 1 1 1 1 1 6 1 [815] 6 6 1 1 1 1 1 1 1 6 1 6 1 1 1 1 1 1 1 1 Levels: 1 2 3 4 5 6 this is the correct output extinction [1] 3 1 6 4 2 1 6 1 2 6 6 1 5 1 1 1 1 6 2 4 1 4 1 3 6 3 5 1 4 3 3 5 3 3 3 5 1 6 6 3 6 6 1 1 5 1 5 3 3 6 1 1 6 2 6 1 3 1 2 6 1 1 3 4 6 3 3 1 1 3 6 3 1 6 [75] 3 4 1 1 1 3 1 3 6 1 2 1 6 6 3 6 6 5 1 1 4 6 3 1 2 1 1 1 4 6 3 1 2 6 1 3 2 1 1 6 1 6 2 3 6 2 4 6 6 5 2 4 6 6 1 2 6 6 4 1 1 1 5 5 1 1 1 6 6 1 1 1 2 6 [149] 4 5 6 6 6 1 3 1 4 3 1 3 6 1 3 6 1 2 1 3 4 1 3 1 6 3 6 1 6 2 1 3 4 6 6 1 1 1 1 5 1 2 1 6 3 1 6 6 1 1 6 6 1 1 2 2 4 1 4 3 5 4 1 4 3 6 1 2 6 1 1 5 6 1 [223] 1 2 1 6 3 6 1 1 1 6 1 6 3 4 1 5 5 6 3 2 3 6 1 1 1 5 1 6 6 1 4 3 1 4 6 1 3 5 1 4 4 6 2 6 6 1 6 3 2 6 3 1 4 1 5 6 1 6 4 5 4 1 4 1 4 1 1 1 4 1 6 6 6 6 [297] 1 1 1 1 6 1 3 2 1 6 6 6 1 1 3 3 1 4 1 6 2 2 6 3 1 1 6 6 2 6 1 1 6 6 1 2 1 6 4 1 4 1 1 5 1 3 6 6 1 1 3 6 6 6 1 6 1 1 6 3 4 6 1 6 1 3 1 2 1 1 6 6 1 1 [371] 1 6 3 1 1 2 5 1 1 1 1 1 6 1 1 1 1 2 4 6 1 5 1 1 6 6 3 1 6 1 6 6 4 1 1 4 2 6 1 1 3 6 1 2 6 1 2 1 1 4 5 1 6 1 6 1 1 6 1 3 6 1 1 1 3 1 2 6 3 5 2 1 3 1 [445] 3 2 4 1 6 6 1 4 5 2 1 6 1 3 4 3 1 1 1 6 1 3 5 6 1 1 1 6 1 1 2 1 1 1 1 2 4 1 1 1 1 3 1 1 1 1 1 1 6 6 4 1 1 2 1 6 6 6 6 6 6 6 1 3 4 3 1 1 1 1 2 1 5 1 [519] 1 6 1 1 1 6 6 1 3 5 2 1 1 1 1 4 3 6 2 1 3 1 1 1 6 3 3 4 1 5 6 6 1 2 6 3 1 5 6 1 1 6 1 1 1 1 1 3 1 1 1 3 6 1 1 1 6 6 1 3 1 2 6 1 6 2 3 1 6 1 6 6 2 4 [593] 1 1 5 3 1 1 1 6 3 1 6 6 4 4 1 6 6 6 1 1 6 1 6 3 1 1 2 1 6 1 1 1 6 1 3 6 6 1 3 3 2 2 6 1 1 1 1 6 1 3 1 2 1 6 1 4 1 1 1 6 1 4 4 1 6 3 1 6 1 3 1 1 2 6 [667] 5 1 2 1 1
Re: [R] String split and concatenation
On Wed, Sep 29, 2010 at 4:15 AM, Steven Kang stochastick...@gmail.com wrote:
x - rep(letters[1:3], 2)
Are there any ways to transform assign the above as the one shown below
to an object? (in exact format; i.e length of 1 class of character),
i.e
x
('a', 'b', 'c', 'a', 'b', 'c')
Highly appreciate for any advice.
Here are a few variations. They all use paste (or the paste0 wrapper)
and toString. The last one uses sQuote to do the quoting, turning off
fancy quotes so that ordinary single quotes are used.
# 1
paste((, toString(paste(', x, ', sep = )), ), sep = )
#2
library(gsubfn) # paste0
paste0((, toString(paste0(', x, ')), ))
# 3
P - function(x, pre = ', post = pre) paste(pre, x, post, sep = )
P(toString(P(x)), (, ))
# 4
old - options(useFancyQuotes = FALSE)
paste((, toString(sQuote(x)), ), sep = )
options(old)
--
Statistics Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at 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] Fitting a half-ellipse curve
Hello Michael,
Thanks very much for the hint to my problem this helps really a lot.
Cheers
Niklaus
--
Niklaus Hürlimann
Doctorant-PhD
Université de Lausanne
Institut de Minéralogie et Géochimie
L'Anthropole
CH-1015 Lausanne
Suisse
E-mail: niklaus.hurlim...@unil.ch
Tel:+41(0)21 692 4452
---BeginMessage---
Hello Niklaus,
I'm not sure if the following is the sort of thing you are looking for (?)
You can fit an ellipse to your data using a deterministic least
squares method. The following is a function that I use to do this...
fit.ellipse - function (x, y = NULL)
{
# Least squares fitting of an ellipse to point data
# using the algorithm described in:
# Radim Halir Jan Flusser. 1998.
# Numerically stable direct least squares fitting of ellipses.
# Proceedings of the 6th International Conference in Central Europe
# on Computer Graphics and Visualization. WSCG '98, p. 125-132
#
# Adapted from the original Matlab code by Michael Bedward
# michael.bedw...@gmail.com
#
# Arguments:
# x, y - the x and y coordinates of the data points
#
# Returns a list with the following elements:
#
# coef - coefficients of the ellipse as described by the general
#quadratic: ax^2 + bxy + cy^2 + dx + ey + f = 0
#
# center - center x and y
#
# major - major semi-axis length
#
# minor - minor semi-axis length
#
EPS - 1.0e-8
dat - xy.coords(x, y)
D1 - cbind(dat$x * dat$x, dat$x * dat$y, dat$y * dat$y)
D2 - cbind(dat$x, dat$y, 1)
S1 - t(D1) %*% D1
S2 - t(D1) %*% D2
S3 - t(D2) %*% D2
T - -solve(S3) %*% t(S2)
M - S1 + S2 %*% T
M - rbind(M[3,] / 2, -M[2,], M[1,] / 2)
evec - eigen(M)$vec
cond - 4 * evec[1,] * evec[3,] - evec[2,]^2
a1 - evec[, which(cond 0)]
f - c(a1, T %*% a1)
names(f) - letters[1:6]
# calculate the center and lengths of the semi-axes
b2 - f[2]^2 / 4
center - c((f[3] * f[4] / 2 - b2 * f[5]), (f[1] * f[5] / 2 - f[2] *
f[4] / 4)) / (b2 - f[1] * f[3])
names(center) - c(x, y)
num - 2 * (f[1] * f[5]^2 / 4 + f[3] * f[4]^2 / 4 + f[6] * b2 -
f[2]*f[4]*f[5]/4 - f[1]*f[3]*f[6])
den1 - (b2 - f[1]*f[3])
den2 - sqrt((f[1] - f[3])^2 + 4*b2)
den3 - f[1] + f[3]
semi.axes - sqrt(c( num / (den1 * (den2 - den3)), num / (den1 *
(-den2 - den3)) ))
# calculate the angle of rotation
term - (f[1] - f[3]) / f[2]
angle - atan(1 / term) / 2
list(coef=f, center = center, major = max(semi.axes), minor =
min(semi.axes), angle = unname(angle))
}
There are quite a few functions / packages to draw ellipses in R, but
the following is will work directly with the output of the above
function...
get.ellipse - function ( fit, n=360 )
{
# Calculate points on an ellipse described by
# the fit argument as returned by fit.ellipse
#
# n is the number of points to render
tt - seq(0, 2*pi, length=n)
sa - sin(fit$angle)
ca - cos(fit$angle)
ct - cos(tt)
st - sin(tt)
x - fit$center[1] + fit$maj * ct * ca - fit$min * st * sa
y - fit$center[2] + fit$maj * ct * sa + fit$min * st * ca
cbind(x=x, y=y)
}
So if your data were in a matrix or data.frame X...
efit - fit.ellipse( X )
e - get.ellipse( efit )
plot(X)
lines(e, col=red)
Hope this helps,
Michael
On 29 September 2010 23:45, Niklaus Hurlimann niklaus.hurlim...@unil.ch wrote:
Dear mailing list,
I have following array:
X2 Y2
[1,] 422.7900 6.0
[2,] 469.8007 10.5
[3,] 483.9428 11.0
[4,] 532.4917 25.5
[5,] 596.1942 33.5
[6,] 630.8496 40.5
[7,] 733.2996 45.0
[8,] 946.4779 32.0
[9,] 996.8068 35.5
[10,] 1074.3310 23.0
I do afterwards the following:
plot.new()
plot.window(xlim=c(min(X1)-50,max(X1)+50),
ylim=c(min(Y1)-25,max(Y1)+25))
axis(2, cex.axis=1.2)
axis(1, cex.axis=1.2)
points(X2, Y2, type=p, pch=0, cex=1.2, col=black)
title(main=Dyke Geometry Along Strike, cex.main=1.2, font.main=4)
title(xlab=distance [m], cex.lab=1.2)
title(ylab=half-thickness [m], cex.lab=1.2)
box()
I would like curve fitting where I proceeded with a non
linear-regression with the according function below:
nls(formula = Y2 ~ abs(b*abs(1-X2^2/a^2)^(1/2)), start = list( a=282,
b=40))
The formula should give the y-positive part of an ellipse in my opinion
or I might be completely wrong.
In the end I would like to fit a curve of half an ellipse through the
data. I might could do this as well with a 2nd order polynomial fit. I
am grateful for any suggestions and comments to my problem.
Cheers
--
Niklaus Hürlimann
Doctorant-PhD
Université de Lausanne
Institut de Minéralogie et Géochimie
L'Anthropole
CH-1015 Lausanne
Suisse
E-mail: niklaus.hurlim...@unil.ch
Tel:+41(0)21 692 4452
[[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] AIC for tweedie glm
Many thanks! Using a fresh session didn't work, but I upgraded to R.2.11.1 and it runs fine now. bbolker wrote: eleadbeater e.leadbeater at sussex.ac.uk writes: Dear R-users, I'm trying to model some data using a tweedie GLM approach. My response variable is the number of pupae that are the offspring of a subordinate wasp on a wasp's nest. However, they're not count data- for each nest, I only know the mean number of pupae per subordinate, which is continous. The data also contain a high proportion of zeros. This worked fine, and gave results I expected, but I don't know what the best method is to evaluate the fit of the model. I am used to using AIC to compare models. A site search turned up AICtweedie, within the tweedie package, but I get the following message: Error: could not find function AICtweedie when I try to use this command, even though tweedie and statmod are both loaded. I've also read that AIC can be calculated using dtweedie, but I'm a beginner and so, despite lots of searching, I'm not sure how. I'm sorry to ask a basic statistics rather than programming question, but I'm really stuck. Could anyone advise me on the best way to assess goodness-of-fit for this type of model, in order to compare models? Everything you're saying sounds sensible. The only (!) problem is that your problem isn't reproducible (for me at least). What happens if you run library(tweedie) example(AICtweedie) from a fresh R session (possibly using R --vanilla)? What are the results of sessionInfo() ? Works for me. Ben Bolker __ 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. -- View this message in context: http://r.789695.n4.nabble.com/AIC-for-tweedie-glm-tp2720813p2720934.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] R crashes when loading rgl package before minqa package
You data is not good enough for the model that you are trying to fit (or
conversely, the model is not appropriate for the data).
Some of the parameters in your model will not be estimable because there is no
information in the data.
See the following:
xvals =
c(1,2,4,5,7,8,9,10,11,12,14,15,16,18,19,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36)
yvals =
c(857.7597,975.8624,978.2655,979.3034,965.5919,983.8946,992.2512,992.1178,979.5379,974.4269,968.4113,991.5210,977.3361,985.7800,975.5220,974.6880,973.8102,980.7295,982.0034,984.7993,978.4948,970.4351,969.0718,983.7892,976.3637,980.7833,987.1665,976.6000,975.1332,971.0757,989.4693)
initpar = c(-5.1471384, -3861.8905839, 979.2616002, 0.2572355, 27.5705764)
fn - function(x) {
yft = x[2] + (x[3] - x[2])/((1 + exp(x[1] * (log(xvals) - log(x[4]^x[5])
return(yft)
}
plot(xvals, yvals, type=p)
lines(xvals, fn(initpar))
Hope this helps,
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Ravi Varadhan rvarad...@jhmi.edu
Date: Thursday, September 30, 2010 10:15 am
Subject: Re: [R] R crashes when loading rgl package before minqa package
To: Gaspard Lequeux gaspard.lequ...@biomath.ugent.be
Cc: r-help@r-project.org
I get this on Windows (it does not crash):
library(minqa)
library(rgl)
newuoa(initpar, optimft)
Error in newuoa(initpar, optimft) :
non-finite x values not allowed in calfun
In addition: Warning message:
In log(x[4]) : NaNs produced
This tells me that you should be constraining your parameter x[4]
(may be even x[5]) to be non-negative:
Here is what I get with `bobyqa':
bobyqa(initpar, optimft, lower=c(-Inf, -Inf, -Inf, 0, 0))
parameter estimates: -5.311767080681, -3861.89005072333,
979.239647766226, 0.268156271922112, 27.6418856936228
objective: 1457.20987728737
number of function evaluations: 78
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Gaspard Lequeux gaspard.lequ...@biomath.ugent.be
Date: Wednesday, September 29, 2010 11:40 am
Subject: [R] R crashes when loading rgl package before minqa package
To: r-help@r-project.org
Hej,
Calling newuoa (from the minqa package) makes R crash when the
package rgl is loaded first. This however only on certain selected
data.
The data used for testing (saved to 'bugs.R'):
xvals =
c(1,2,4,5,7,8,9,10,11,12,14,15,16,18,19,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36)
yvals =
c(857.7597,975.8624,978.2655,979.3034,965.5919,983.8946,992.2512,992.1178,979.5379,974.4269,968.4113,991.5210,977.3361,985.7800,975.5220,974.6880,973.8102,980.7295,982.0034,984.7993,978.4948,970.4351,969.0718,983.7892,976.3637,980.7833,987.1665,976.6000,975.1332,971.0757,989.4693)
initpar = c(-5.1471384, -3861.8905839, 979.2616002, 0.2572355, 27.5705764)
optimft - function(x) {
yft = x[2] + (x[3] - x[2])/((1 + exp(x[1] * (log(xvals) -
log(x[4]^x[5])
return(sum((yvals - yft)^2))
}
Sequence of commands needed to make the bug appear:
Start R
source('bugs.R')
library(minqa)
library(rgl)
newuoa(initpar, optimft)
= OK
Start R
source('bugs.R')
library(rgl)
library(minqa)
newuoa(initpar, optimft)
= Crash: segfault: address 0x18, cause 'memory not mapped'
I found the bug using the package qpcR, where rgl is loaded when
loading qpcR while minqa is only loaded later, when needed.
Running on Debian squeeze 64 bit.
R version: R version 2.11.1 (2010-05-31) x86_64-pc-linux-gnu
rgl version: 0.91
minqa version: 1.1.9
Rcpp version: 0.8.6 (loaded by minqa)
Kind regards,
Gaspard Lequeux
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
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] panel.pairs in splom
Hello,
I have a customized pairs () fonction as follows that displays correctely my
data.
panel.cor1 - function (x, y, digits=2, prefix=)
{
usr - par(usr); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r - cor(x, y,use=pairwise.complete.obs, method = meth)
if (r0) {alt-greater
} else alt-less
co -cor.test(x,y,method = meth,alternative=alt)$p.value
if (co0.05 ) {colo-red
} else colo-black
txt - format(c(r, 0.123456789), digits=digits)[1]
txt - paste(prefix, txt, sep=)
cex.cor - 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = abs(cex.cor * r),col=colo)
}
panel.line1 - function (x, y, col = par(col), bg = NA, pch =
par(pch), cex = 1)
{
points(x, y, pch = pch, col = colour, bg = bg, cex = cex)
ok - is.finite(x) is.finite(y)
if (any(ok))
abline(lsfit(x,y,intercept = TRUE), col= red )
}
pairs(temp.df,
lower.panel=panel.line1,upper.panel=panel.cor1 )
However I decided to add a green line to my display and I was not able to do
it with the pairs() function, I used insted splom().
But I dont'know how to change panel.pairs inside superpanel in order to use
my panel.line1 and panel.cor1 functions.
splom(~temp.df, aspect=fill, varnames=paste(coord, 1:10, sep=),
xlab=, pscales=0, varname.cex=0.6, cex=0.2,
superpanel = function(...) {
panel.pairs(...)
panel.abline(h =6.5, v = 6.5, col = green, lwd = 4)
})
I very much appreciate if someone could help me.
All the best,
Isabel
--
View this message in context:
http://r.789695.n4.nabble.com/panel-pairs-in-splom-tp2720948p2720948.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] panel.pairs in splom
Sorry, I forget to define temp.df - sapply(1:10, function(i) rnorm(20, 0,1)) Best, Isabel -- View this message in context: http://r.789695.n4.nabble.com/panel-pairs-in-splom-tp2720948p2720968.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] panel.pairs in splom
Indeed, some commands are missing. Sorry.
My function is as follows,
panel.cor1 - function (x, y, digits=2, prefix=)
{
usr - par(usr); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r - cor(x, y,use=pairwise.complete.obs, method = meth)
if (r0) {alt-greater
} else alt-less
co -cor.test(x,y,method = meth,alternative=alt)$p.value
if (co0.05 ) {colo-red
} else colo-black
txt - format(c(r, 0.123456789), digits=digits)[1]
txt - paste(prefix, txt, sep=)
cex.cor - 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = abs(cex.cor * r),col=colo)
}
panel.line1 - function (x, y, col = par(col), bg = NA, pch =
par(pch), cex = 1)
{
points(x, y, pch = pch, col = colour, bg = bg, cex = cex)
ok - is.finite(x) is.finite(y)
if (any(ok))
abline(lsfit(x,y,intercept = TRUE), col= red )
}
temp.df - sapply(1:10, function(i) rnorm(20, 0,1))
meth-kendall
colour-c(rep(blue, 20 ))
pairs(temp.df,
lower.panel=panel.line1,upper.panel=panel.cor1 )
I would like to use panel.line1 and panel.cor1 inside superpanel of splom
function below
require(lattice)
splom(~temp.df, aspect=fill, varnames=paste(coord, 1:10, sep=),
xlab=, pscales=0, varname.cex=0.6, cex=0.2,
superpanel = function(...) {
panel.pairs(...)
panel.abline(h =6.5, v = 6.5, col = green, lwd = 4)
})
Best,
Isabel
--
View this message in context:
http://r.789695.n4.nabble.com/panel-pairs-in-splom-tp2720948p2721009.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] R crashes when loading rgl package before minqa package
Hej,
On Thu, 30 Sep 2010, Ravi Varadhan wrote:
I get this on Windows (it does not crash):
library(minqa)
library(rgl)
newuoa(initpar, optimft)
Error in newuoa(initpar, optimft) :
non-finite x values not allowed in calfun
In addition: Warning message:
In log(x[4]) : NaNs produced
Does it crash when you load first rgl and then only minqa? Like this:
library(rgl)
library(minqa)
newuoa(initpar, optimft)
/Gaspard
This tells me that you should be constraining your parameter x[4] (may be even
x[5]) to be non-negative:
Here is what I get with `bobyqa':
bobyqa(initpar, optimft, lower=c(-Inf, -Inf, -Inf, 0, 0))
parameter estimates: -5.311767080681, -3861.89005072333, 979.239647766226,
0.268156271922112, 27.6418856936228
objective: 1457.20987728737
number of function evaluations: 78
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Gaspard Lequeux gaspard.lequ...@biomath.ugent.be
Date: Wednesday, September 29, 2010 11:40 am
Subject: [R] R crashes when loading rgl package before minqa package
To: r-help@r-project.org
Hej,
Calling newuoa (from the minqa package) makes R crash when the
package rgl is loaded first. This however only on certain selected data.
The data used for testing (saved to 'bugs.R'):
xvals =
c(1,2,4,5,7,8,9,10,11,12,14,15,16,18,19,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36)
yvals =
c(857.7597,975.8624,978.2655,979.3034,965.5919,983.8946,992.2512,992.1178,979.5379,974.4269,968.4113,991.5210,977.3361,985.7800,975.5220,974.6880,973.8102,980.7295,982.0034,984.7993,978.4948,970.4351,969.0718,983.7892,976.3637,980.7833,987.1665,976.6000,975.1332,971.0757,989.4693)
initpar = c(-5.1471384, -3861.8905839, 979.2616002, 0.2572355, 27.5705764)
optimft - function(x) {
yft = x[2] + (x[3] - x[2])/((1 + exp(x[1] * (log(xvals) - log(x[4]^x[5])
return(sum((yvals - yft)^2))
}
Sequence of commands needed to make the bug appear:
Start R
source('bugs.R')
library(minqa)
library(rgl)
newuoa(initpar, optimft)
= OK
Start R
source('bugs.R')
library(rgl)
library(minqa)
newuoa(initpar, optimft)
= Crash: segfault: address 0x18, cause 'memory not mapped'
I found the bug using the package qpcR, where rgl is loaded when
loading qpcR while minqa is only loaded later, when needed.
Running on Debian squeeze 64 bit.
R version: R version 2.11.1 (2010-05-31) x86_64-pc-linux-gnu
rgl version: 0.91
minqa version: 1.1.9
Rcpp version: 0.8.6 (loaded by minqa)
Kind regards,
Gaspard Lequeux
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
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] R crashes when loading rgl package before minqa package
Hej,
On Thu, 30 Sep 2010, Ravi Varadhan wrote:
You data is not good enough for the model that you are trying to fit (or
conversely, the model is not appropriate for the data).
Yes. I agree. However, this should not cause R to crash. (Around 100.000
models are fitted, and it is not possible to look at the data by hand
first. Bogus results are filtered out afterwards.)
/Gaspard
Some of the parameters in your model will not be estimable because there is no
information in the data.
See the following:
xvals =
c(1,2,4,5,7,8,9,10,11,12,14,15,16,18,19,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36)
yvals =
c(857.7597,975.8624,978.2655,979.3034,965.5919,983.8946,992.2512,992.1178,979.5379,974.4269,968.4113,991.5210,977.3361,985.7800,975.5220,974.6880,973.8102,980.7295,982.0034,984.7993,978.4948,970.4351,969.0718,983.7892,976.3637,980.7833,987.1665,976.6000,975.1332,971.0757,989.4693)
initpar = c(-5.1471384, -3861.8905839, 979.2616002, 0.2572355, 27.5705764)
fn - function(x) {
yft = x[2] + (x[3] - x[2])/((1 + exp(x[1] * (log(xvals) - log(x[4]^x[5])
return(yft)
}
plot(xvals, yvals, type=p)
lines(xvals, fn(initpar))
Hope this helps,
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Ravi Varadhan rvarad...@jhmi.edu
Date: Thursday, September 30, 2010 10:15 am
Subject: Re: [R] R crashes when loading rgl package before minqa package
To: Gaspard Lequeux gaspard.lequ...@biomath.ugent.be
Cc: r-help@r-project.org
I get this on Windows (it does not crash):
library(minqa)
library(rgl)
newuoa(initpar, optimft)
Error in newuoa(initpar, optimft) :
non-finite x values not allowed in calfun
In addition: Warning message:
In log(x[4]) : NaNs produced
This tells me that you should be constraining your parameter x[4]
(may be even x[5]) to be non-negative:
Here is what I get with `bobyqa':
bobyqa(initpar, optimft, lower=c(-Inf, -Inf, -Inf, 0, 0))
parameter estimates: -5.311767080681, -3861.89005072333,
979.239647766226, 0.268156271922112, 27.6418856936228
objective: 1457.20987728737
number of function evaluations: 78
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Gaspard Lequeux gaspard.lequ...@biomath.ugent.be
Date: Wednesday, September 29, 2010 11:40 am
Subject: [R] R crashes when loading rgl package before minqa package
To: r-help@r-project.org
Hej,
Calling newuoa (from the minqa package) makes R crash when the
package rgl is loaded first. This however only on certain selected
data.
The data used for testing (saved to 'bugs.R'):
xvals =
c(1,2,4,5,7,8,9,10,11,12,14,15,16,18,19,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36)
yvals =
c(857.7597,975.8624,978.2655,979.3034,965.5919,983.8946,992.2512,992.1178,979.5379,974.4269,968.4113,991.5210,977.3361,985.7800,975.5220,974.6880,973.8102,980.7295,982.0034,984.7993,978.4948,970.4351,969.0718,983.7892,976.3637,980.7833,987.1665,976.6000,975.1332,971.0757,989.4693)
initpar = c(-5.1471384, -3861.8905839, 979.2616002, 0.2572355, 27.5705764)
optimft - function(x) {
yft = x[2] + (x[3] - x[2])/((1 + exp(x[1] * (log(xvals) -
log(x[4]^x[5])
return(sum((yvals - yft)^2))
}
Sequence of commands needed to make the bug appear:
Start R
source('bugs.R')
library(minqa)
library(rgl)
newuoa(initpar, optimft)
= OK
Start R
source('bugs.R')
library(rgl)
library(minqa)
newuoa(initpar, optimft)
= Crash: segfault: address 0x18, cause 'memory not mapped'
I found the bug using the package qpcR, where rgl is loaded when
loading qpcR while minqa is only loaded later, when needed.
Running on Debian squeeze 64 bit.
R version: R version 2.11.1 (2010-05-31) x86_64-pc-linux-gnu
rgl version: 0.91
minqa version: 1.1.9
Rcpp version: 0.8.6 (loaded by minqa)
Kind regards,
Gaspard Lequeux
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
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] getting the output after bootstraping
Thanks to the help of people from this forum I was able to bootstrap my data and then apply a model to it. Thanks for all your help. Everything worked out well, but I am having a difficult time getting the new parameter values. I bootstrapped the data 300 times and I want to get the 300 sets of parameter estimates and plot them in Excel. Here is my code: par-list(Linf=700, K=0.3, to=-0.1)#starting values for parameter estimates vb-nls(length~Linf*(1-exp(-K*(age-to))), start=par, data=small) #von Bertalanffy growth mode boo- nlsBoot(vb, niter=300) #This bootstraps my data 300 times and applies the growth model I can get R to plot the 300 parameter estimates if I do the following function: plot(boo) However, I have tried different print function options but have not had any luck getting the 300 parameter estimates. Any advice would be greatly appreciated. Thanks. Mike [[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] Polar coordinates - contour plots
Like Todor I have been trying to make a contour plot using polar coordinates, unfortunately without success. The problem with converting to a cartesian system and plotting using e.g. filled.contour is that this function requires the same amount and value of x values for each y-value and also the different way around. I have a very simple dataset with 360 values at 10 equidistant radii (say seq(10)) for 36 polar angles (say seq(36)*10/(2*pi)). There must be a simple way to do this in R! I can imagine more people having had the same problem and a solution is at hand.. However I have not found an R-package suitable to do the job. Can anyone help out? Cheers, Koen -- View this message in context: http://r.789695.n4.nabble.com/Polar-coordinates-contour-plots-tp876469p2726219.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] R crashes when loading rgl package before minqa package
No. It still does not crash in Windows.
library(rgl)
library(minqa)
Loading required package: Rcpp
newuoa(initpar, optimft)
Error in newuoa(initpar, optimft) :
non-finite x values not allowed in calfun
In addition: Warning message:
In log(x[4]) : NaNs produced
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Gaspard Lequeux gaspard.lequ...@biomath.ugent.be
Date: Thursday, September 30, 2010 11:43 am
Subject: Re: [R] R crashes when loading rgl package before minqa package
To: r-help@r-project.org
Hej,
On Thu, 30 Sep 2010, Ravi Varadhan wrote:
I get this on Windows (it does not crash):
library(minqa)
library(rgl)
newuoa(initpar, optimft)
Error in newuoa(initpar, optimft) :
non-finite x values not allowed in calfun
In addition: Warning message:
In log(x[4]) : NaNs produced
Does it crash when you load first rgl and then only minqa? Like this:
library(rgl)
library(minqa)
newuoa(initpar, optimft)
/Gaspard
This tells me that you should be constraining your parameter x[4]
(may be even x[5]) to be non-negative:
Here is what I get with `bobyqa':
bobyqa(initpar, optimft, lower=c(-Inf, -Inf, -Inf, 0, 0))
parameter estimates: -5.311767080681, -3861.89005072333,
979.239647766226, 0.268156271922112, 27.6418856936228
objective: 1457.20987728737
number of function evaluations: 78
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Gaspard Lequeux gaspard.lequ...@biomath.ugent.be
Date: Wednesday, September 29, 2010 11:40 am
Subject: [R] R crashes when loading rgl package before minqa package
To: r-help@r-project.org
Hej,
Calling newuoa (from the minqa package) makes R crash when the
package rgl is loaded first. This however only on certain selected
data.
The data used for testing (saved to 'bugs.R'):
xvals =
c(1,2,4,5,7,8,9,10,11,12,14,15,16,18,19,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36)
yvals =
c(857.7597,975.8624,978.2655,979.3034,965.5919,983.8946,992.2512,992.1178,979.5379,974.4269,968.4113,991.5210,977.3361,985.7800,975.5220,974.6880,973.8102,980.7295,982.0034,984.7993,978.4948,970.4351,969.0718,983.7892,976.3637,980.7833,987.1665,976.6000,975.1332,971.0757,989.4693)
initpar = c(-5.1471384, -3861.8905839, 979.2616002, 0.2572355, 27.5705764)
optimft - function(x) {
yft = x[2] + (x[3] - x[2])/((1 + exp(x[1] * (log(xvals) -
log(x[4]^x[5])
return(sum((yvals - yft)^2))
}
Sequence of commands needed to make the bug appear:
Start R
source('bugs.R')
library(minqa)
library(rgl)
newuoa(initpar, optimft)
= OK
Start R
source('bugs.R')
library(rgl)
library(minqa)
newuoa(initpar, optimft)
= Crash: segfault: address 0x18, cause 'memory not mapped'
I found the bug using the package qpcR, where rgl is loaded when
loading qpcR while minqa is only loaded later, when needed.
Running on Debian squeeze 64 bit.
R version: R version 2.11.1 (2010-05-31) x86_64-pc-linux-gnu
rgl version: 0.91
minqa version: 1.1.9
Rcpp version: 0.8.6 (loaded by minqa)
Kind regards,
Gaspard Lequeux
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
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] relevance vector machines for classification
The rvm function from the kernlab package can only be used for regression at the present time. In fact, in the description in the kernlab documentation for the type argument for rvm says, type rvm can only be used for regression at the moment. Are there any R packages that do classification with relevance vector machines? Keith This e-mail and any files transmitted with it are confidential and are intended solely for the use of the individual or entity to whom it is addressed. If you are not the intended recipient or the person responsible for delivering the e-mail to the intended recipient, be advised that you have received this e-mail in error, and that any use, dissemination, forwarding, printing, or copying of this e-mail is strictly prohibited. If you received this e-mail in error, please return the e-mail to the sender at Merrick Bank and delete it from your computer. Although Merrick Bank attempts to sweep e-mail and attachments for viruses, it does not guarantee that either are virus-free and accepts no liability for any damage sustained as a result of viruses. [[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] cor() alternative for huge data set
Hi Jyotasana, if I understand your aim correctly, you want to find correlated sets (clusters) of genes, and then find those clusters that are differentially expressed? You can do that with WGCNA, or you can just use the projectiveKMeans for splitting your probes into blocks and then feed each block into coXpress. The correlations of probes in different blocks will be very small and can be considered zero. Peter On Thu, Sep 30, 2010 at 5:05 AM, Jyotasana Gulati jgul...@ice.mpg.de wrote: Peter, Many thank for suggesting me this package. I very much believe that this will help me. But I was trying to correlate all probes(correlation between entities not variables) to calculate differentially coexpressed gene sets using package coXpress in R. I could not reduce the number on the basis of intensity, since most of the genes are down regulated and upregulated in treated conditions, so they are of my interest and cannot be removed from control samples(since I have to compare both). can you further suggest me an alternative for differentially coexpression analysis, since this is what I need to know the most-- the sets which are behaving differently across conditions. Has any one ever used this package--coXpress?? Regards .. Jyotasana - Original Message - From: Peter Langfelder peter.langfel...@gmail.com To: Jyotasana Gulati jgul...@ice.mpg.de Cc: r-help@r-project.org Sent: Thursday, September 30, 2010 4:05:44 AM Subject: Re: [R] cor() alternative for huge data set On Wed, Sep 29, 2010 at 1:27 PM, Jyotasana Gulati jgul...@ice.mpg.de wrote: Hi, I am have a data set of around 43000 probes(rows), and have to calculate correlation matrix. When I run cor function in R, its throwing an error message of RAM shortage which was obvious for such huge number of rows. I am not getting a logical way to cut off this huge number of entities, is there an alternative to pearson correlation or with other dist() methods calculation(euclidean) that can be run on such a huge data set?? Every help will be appreciated. Hmm... Are you calculating a correlation of 43000 probes, or of some number of samples across 43000 probes? If the former, read below. If the latter, I'm surprised you are running out of memory. Issuing garbage collection (gc()) before the calculation, closing all other programs, removing all other large objects from the R workspace etc. may help. If you really need the 43k times 43k correlation matrix of your 43k probes, read on. [Disclosure: this is a shameless plug for the package WGCNA (Weighted Gene Co-expression Network Analysis, also known as Weighted Correlation Network Analysis), from the package author, namely me.] First, since the distance matrix will be huge, you will not gain using other distance methods either. Second, depending on what you want to do with the 43k probes, the package WGCNA may help you. It has methods for creating correlation networks among a large number of probes. The idea is to pre-cluster the probes using what I call projective K-means, function projectiveKMeans. The pre-clustering will return what we call blocks of probes (or genes). We assume (this is a big assumption) that correlations among probes belonging to different blocks can be neglected. Then we treat each block separately for network construction (or, in your case, possibly simple calculation of correlation). Although this isn't strictly an R topic but rather microarray analysis issue, in my experience it is often useful to filter out probes before actually calculating and interpreting large correlation matrices. In conjunction with filtering, it can be advantageous to only keep one probe per gene (presumably there is more than one probe per gene in you data set). The filtering criterion varies from analysis to analysis, but if your data represent intensities, it is often a good idea to throw away probes whose intensity is always low, because such signals are mostly noise. If you decide to check out WGCNA, look at http://www.genetics.ucla.edu/labs/horvath/CoexpressionNetwork/Rpackages/WGCNA/. Peter __ 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] R crashes when loading rgl package before minqa package
I also cannot reproduce the crash.
sessionInfo()
R version 2.11.1 (2010-05-31)
x86_64-unknown-linux-gnu
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] minqa_1.1.9 Rcpp_0.8.6 rgl_0.91
On Thu, 30 Sep 2010, Ravi Varadhan wrote:
No. It still does not crash in Windows.
library(rgl)
library(minqa)
Loading required package: Rcpp
newuoa(initpar, optimft)
Error in newuoa(initpar, optimft) :
non-finite x values not allowed in calfun
In addition: Warning message:
In log(x[4]) : NaNs produced
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Gaspard Lequeux gaspard.lequ...@biomath.ugent.be
Date: Thursday, September 30, 2010 11:43 am
Subject: Re: [R] R crashes when loading rgl package before minqa package
To: r-help@r-project.org
Hej,
On Thu, 30 Sep 2010, Ravi Varadhan wrote:
I get this on Windows (it does not crash):
library(minqa)
library(rgl)
newuoa(initpar, optimft)
Error in newuoa(initpar, optimft) :
non-finite x values not allowed in calfun
In addition: Warning message:
In log(x[4]) : NaNs produced
Does it crash when you load first rgl and then only minqa? Like this:
library(rgl)
library(minqa)
newuoa(initpar, optimft)
/Gaspard
This tells me that you should be constraining your parameter x[4]
(may be even x[5]) to be non-negative:
Here is what I get with `bobyqa':
bobyqa(initpar, optimft, lower=c(-Inf, -Inf, -Inf, 0, 0))
parameter estimates: -5.311767080681, -3861.89005072333,
979.239647766226, 0.268156271922112, 27.6418856936228
objective: 1457.20987728737
number of function evaluations: 78
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Gaspard Lequeux gaspard.lequ...@biomath.ugent.be
Date: Wednesday, September 29, 2010 11:40 am
Subject: [R] R crashes when loading rgl package before minqa package
To: r-help@r-project.org
Hej,
Calling newuoa (from the minqa package) makes R crash when the
package rgl is loaded first. This however only on certain selected
data.
The data used for testing (saved to 'bugs.R'):
xvals =
c(1,2,4,5,7,8,9,10,11,12,14,15,16,18,19,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36)
yvals =
c(857.7597,975.8624,978.2655,979.3034,965.5919,983.8946,992.2512,992.1178,979.5379,974.4269,968.4113,991.5210,977.3361,985.7800,975.5220,974.6880,973.8102,980.7295,982.0034,984.7993,978.4948,970.4351,969.0718,983.7892,976.3637,980.7833,987.1665,976.6000,975.1332,971.0757,989.4693)
initpar = c(-5.1471384, -3861.8905839, 979.2616002, 0.2572355, 27.5705764)
optimft - function(x) {
yft = x[2] + (x[3] - x[2])/((1 + exp(x[1] * (log(xvals) -
log(x[4]^x[5])
return(sum((yvals - yft)^2))
}
Sequence of commands needed to make the bug appear:
Start R
source('bugs.R')
library(minqa)
library(rgl)
newuoa(initpar, optimft)
= OK
Start R
source('bugs.R')
library(rgl)
library(minqa)
newuoa(initpar, optimft)
= Crash: segfault: address 0x18, cause 'memory not mapped'
I found the bug using the package qpcR, where rgl is loaded when
loading qpcR while minqa is only loaded later, when needed.
Running on Debian squeeze 64 bit.
R version: R version 2.11.1 (2010-05-31) x86_64-pc-linux-gnu
rgl version: 0.91
minqa version: 1.1.9
Rcpp version: 0.8.6 (loaded by minqa)
Kind regards,
Gaspard Lequeux
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
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-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,
Re: [R] Understanding linear contrasts in Anova using R
Hi Professor Howell, I think the issue here is simply in the assumption that the regression coefficients will always be equal to the product of the means and the contrast codes. I tend to think of regression coefficients as the quotient of the covariance of x and y divided by the variance of x, and this definition agrees with the coefficients calculated by lm(). See below for a long-winded example. On Wed, Sep 29, 2010 at 3:42 PM, David Howell david.how...@uvm.edu wrote: #I am trying to understand how R fits models for contrasts in a #simple one-way anova. This is an example, I am not stupid enough to want #to simultaneously apply all of these contrasts to real data. With a few #exceptions, the tests that I would compute by hand (or by other software) #will give the same t or F statistics. It is the contrast estimates that R produces #that I can't seem to understand. # # In searching for answers to this problem, I found a great PowerPoint slide (I think by John Fox). # The slide pointed to the coefficients, said something like these are coeff. that no one could love, and #then suggested looking at the means to understand where they came from. I have stared # and stared at his means and then my means, but can't find a relationship. # The following code and output illustrates the problem. # Various examples of Anova using R dv - c(1.28, 1.35, 3.31, 3.06, 2.59, 3.25, 2.98, 1.53, -2.68, 2.64, 1.26, 1.06, -1.18, 0.15, 1.36, 2.61, 0.66, 1.32, 0.73, -1.06, 0.24, 0.27, 0.72, 2.28, -0.41, -1.25, -1.33, -0.47, -0.60, -1.72, -1.74, -0.77, -0.41, -1.20, -0.31, -0.74, -0.45, 0.54, -0.98, 1.68, 2.25, -0.19, -0.90, 0.78, 0.05, 2.69, 0.15, 0.91, 2.01, 0.40, 2.34, -1.80, 5.00, 2.27, 6.47, 2.94, 0.47, 3.22, 0.01, -0.66) group - factor(rep(1:5, each = 12)) # Use treatment contrasts to compare each group to the first group. options(contrasts = c(contr.treatment,contr.poly)) # The default model2 - lm(dv ~ group) summary(model2) # Summary table is the same--as it should be # Intercept is Group 1 mean and other coeff. are deviations from that. # This is what I would expect. #summary(model1) # Df Sum Sq Mean Sq F value Pr(F) # group 4 62.46 15.6151 6.9005 0.0001415 *** # Residuals 55 124.46 2.2629 #Coefficients: # Estimate Std. Error t value Pr(|t|) #(Intercept) 1.80250 0.43425 4.151 0.000116 *** #group2 -1.12750 0.61412 -1.836 0.071772 . #group3 -2.71500 0.61412 -4.421 4.67e-05 *** #group4 -1.25833 0.61412 -2.049 0.045245 * #group5 0.08667 0.61412 0.141 0.888288 # Use sum contrasts to compare each group against grand mean. options(contrasts = c(contr.sum,contr.poly)) model3 - lm(dv ~ group) summary(model3) # Again, this is as expected. Intercept is grand mean and others are deviatoions from that. #Coefficients: # Estimate Std. Error t value Pr(|t|) # (Intercept) 0.7997 0.1942 4.118 0.000130 *** # group1 1.0028 0.3884 2.582 0.012519 * # group2 -0.1247 0.3884 -0.321 0.749449 # group3 -1.7122 0.3884 -4.408 4.88e-05 *** # group4 -0.2555 0.3884 -0.658 0.513399 #SO FAR, SO GOOD # IF I wanted polynomial contrasts BY HAND I would use # a(i) = -2 -1 0 1 2 for linear contrast (or some linear function of this ) # Effect = Sum(a(j)M(i)) # where M = mean # Effect(linear) = -2(1.805) -1(0.675) +0(-.912) +1(.544) +2(1.889) = 0.043 # SS(linear) = n*(Effect(linear)^2)/Sum((a(j)^2)) = 12(.043)/10 = .002 # F(linear) = SS(linear)/MS(error) = .002/2.263 = .001 # t(linear) = sqrt(.001) = .031 # To do this in R I would use order.group - ordered(group) model4 - lm(dv~order.group) summary(model4) # This gives: #Coefficients: # Estimate Std. Error t value Pr(|t|) # (Intercept) 0.79967 0.19420 4.118 0.000130 *** # order.group.L 0.01344 0.43425 0.031 0.975422 # order.group.Q 2.13519 0.43425 4.917 8.32e-06 *** # order.group.C 0.11015 0.43425 0.254 0.800703 # order.group^4 -0.79602 0.43425 -1.833 0.072202 . # The t value for linear is same as I got (as are others) but I don't understand # the estimates. The intercept is the grand mean, but I don't see the relationship # of other estimates to that or to the ones I get by hand. # My estimates are the sum of (coeff times means) i.e. 0 (intercept), .0425, 7.989, .3483, -6.66 # and these are not a linear (or other nice pretty) function of est. from R. # OK, let's break it down Means - tapply(dv, order.group, mean) coefs.by.hand.m4 - c(mean(Means), sum(Means*contrasts(order.group)[,1]), sum(Means*contrasts(order.group)[,2]), sum(Means*contrasts(order.group)[,3]),
[R] Can this code be written more efficiently?
Dear users,
I'm working on binary classification problem using Support Vector
Machines (SVM). My objective is to train a series of SVM models on a
grid of hyperparameters and then select those that maximize the AUC
based on an independent validation sample.
My attempted code is shown below. It runs well on small data sets but
when I use it on a slightly larger sample (e.g., my train data is
composed of about 8,000 observations on each class and 21 inputs), it
takes forever to run (more than 1 day already and still running). I'm
wondering if there's any way I can optimize this code. Thanks in advance
for any help.
I'm using 64-bit R 2.11.1 on Win 7.
Start Code
library(e1071)
library(ROCR)
### Create grid of hyperparameters
Gseq - seq(-15,3,2); G - rep(2, length(Gseq)); G - G^Gseq
Cseq - seq(-5,13,2); C - rep(2, length(Cseq)); C - C^Cseq
mygrid - expand.grid(C=C, G=G)
### Train models
svm.models - lapply(1:nrow(mygrid), function(i) {
svm(churn.form, data = mytraindata,
method = C-classification, kernel = radial,
cost = mygrid[i,1], gamma = mygrid[i,2],
probability=TRUE)
})
### Predict on test set
pred.step3 - numeric(length(svm.models))
for (i in 1:length(svm.models)) {
pred.step1 - predict(svm.models[[i]], myvaliddata, decision.values = F,
probability=T)
pred.step2 -
prediction(predictions=attr(pred.step1,probabilities)[,1],
labels=myvaliddata$churn)
pred.step3[i] - performance(pred.step2, auc)@y.values[[1]]
}
pred.step3
End Code
Thanks,
Leo.
___
This e-mail may be privileged and/or confidential, and the sender does not waive
any related rights and obligations. Any distribution, use or copying of this
e-mail or the information
it contains by other than an intended recipient is unauthorized.
If you received this e-mail in error, please advise me (by return e-mail or
otherwise) immediately.
Ce courriel peut contenir des renseignements protégés et confidentiels.
LÂexpéditeur ne renonce pas aux droits et obligations qui sÂy rapportent.
Toute diffusion, utilisation ou copie de ce courriel ou des renseignements
quÂil contient
par une personne autre que le destinataire désigné est interdite.
Si vous recevez ce courriel par erreur, veuillez mÂen aviser immédiatement,
par retour de courriel ou par un autre moyen.
[[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] how to avoid NaN in optim()
hi ,
lik - function(nO, nA, nB, nAB){
loglik - function(par)
{
p=par[1]
q=par[2]
r - 1 - p - q
if (c(p,q,r) rep(0,3) c(p,q,r) rep(1,3) )
{
-(2 * nO * log (r) + nA * log (p^2 + 2 * p * r)
+ nB * log (q^2 + 2 * q * r)
+ nAB * (log(2) +log(p) +log(q)))
}
else
NA
}
loglik
}
i want to maximize this likelihood function over the range (0,0,0) to
(1,1,1). so i tried
optim ( c(0.3,0.3), lik ( 176,182 , 60 ,17) , method = CG)
and obtained the following :
optim(c(0.3,0.3), fn, method=CG)
$par
[1] 0.26444187 0.09316946
$value
[1] 492.5353
$counts
function gradient
130 19
$convergence
[1] 0
$message
NULL
Warning messages:
1: In log(q^2 + 2 * q * r) : NaNs produced
2: In log(q) : NaNs produced
3: In log(q^2 + 2 * q * r) : NaNs produced
4: In log(q) : NaNs produced
please help ...
--
Arindam Fadikar
M.Stat
Indian Statistical Institute.
New Delhi, India
[[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] Understanding linear contrasts in Anova using R
These two resources might also help: http://cran.r-project.org/doc/contrib/Faraway-PRA.pdf http://cran.r-project.org/web/packages/contrast/vignettes/contrast.pdf Max On Thu, Sep 30, 2010 at 1:33 PM, Ista Zahn iz...@psych.rochester.edu wrote: Hi Professor Howell, I think the issue here is simply in the assumption that the regression coefficients will always be equal to the product of the means and the contrast codes. I tend to think of regression coefficients as the quotient of the covariance of x and y divided by the variance of x, and this definition agrees with the coefficients calculated by lm(). See below for a long-winded example. On Wed, Sep 29, 2010 at 3:42 PM, David Howell david.how...@uvm.edu wrote: #I am trying to understand how R fits models for contrasts in a #simple one-way anova. This is an example, I am not stupid enough to want #to simultaneously apply all of these contrasts to real data. With a few #exceptions, the tests that I would compute by hand (or by other software) #will give the same t or F statistics. It is the contrast estimates that R produces #that I can't seem to understand. # # In searching for answers to this problem, I found a great PowerPoint slide (I think by John Fox). # The slide pointed to the coefficients, said something like these are coeff. that no one could love, and #then suggested looking at the means to understand where they came from. I have stared # and stared at his means and then my means, but can't find a relationship. # The following code and output illustrates the problem. # Various examples of Anova using R dv - c(1.28, 1.35, 3.31, 3.06, 2.59, 3.25, 2.98, 1.53, -2.68, 2.64, 1.26, 1.06, -1.18, 0.15, 1.36, 2.61, 0.66, 1.32, 0.73, -1.06, 0.24, 0.27, 0.72, 2.28, -0.41, -1.25, -1.33, -0.47, -0.60, -1.72, -1.74, -0.77, -0.41, -1.20, -0.31, -0.74, -0.45, 0.54, -0.98, 1.68, 2.25, -0.19, -0.90, 0.78, 0.05, 2.69, 0.15, 0.91, 2.01, 0.40, 2.34, -1.80, 5.00, 2.27, 6.47, 2.94, 0.47, 3.22, 0.01, -0.66) group - factor(rep(1:5, each = 12)) # Use treatment contrasts to compare each group to the first group. options(contrasts = c(contr.treatment,contr.poly)) # The default model2 - lm(dv ~ group) summary(model2) # Summary table is the same--as it should be # Intercept is Group 1 mean and other coeff. are deviations from that. # This is what I would expect. #summary(model1) # Df Sum Sq Mean Sq F value Pr(F) # group 4 62.46 15.6151 6.9005 0.0001415 *** # Residuals 55 124.46 2.2629 #Coefficients: # Estimate Std. Error t value Pr(|t|) #(Intercept) 1.80250 0.43425 4.151 0.000116 *** #group2 -1.12750 0.61412 -1.836 0.071772 . #group3 -2.71500 0.61412 -4.421 4.67e-05 *** #group4 -1.25833 0.61412 -2.049 0.045245 * #group5 0.08667 0.61412 0.141 0.888288 # Use sum contrasts to compare each group against grand mean. options(contrasts = c(contr.sum,contr.poly)) model3 - lm(dv ~ group) summary(model3) # Again, this is as expected. Intercept is grand mean and others are deviatoions from that. #Coefficients: # Estimate Std. Error t value Pr(|t|) # (Intercept) 0.7997 0.1942 4.118 0.000130 *** # group1 1.0028 0.3884 2.582 0.012519 * # group2 -0.1247 0.3884 -0.321 0.749449 # group3 -1.7122 0.3884 -4.408 4.88e-05 *** # group4 -0.2555 0.3884 -0.658 0.513399 #SO FAR, SO GOOD # IF I wanted polynomial contrasts BY HAND I would use # a(i) = -2 -1 0 1 2 for linear contrast (or some linear function of this ) # Effect = Sum(a(j)M(i)) # where M = mean # Effect(linear) = -2(1.805) -1(0.675) +0(-.912) +1(.544) +2(1.889) = 0.043 # SS(linear) = n*(Effect(linear)^2)/Sum((a(j)^2)) = 12(.043)/10 = .002 # F(linear) = SS(linear)/MS(error) = .002/2.263 = .001 # t(linear) = sqrt(.001) = .031 # To do this in R I would use order.group - ordered(group) model4 - lm(dv~order.group) summary(model4) # This gives: #Coefficients: # Estimate Std. Error t value Pr(|t|) # (Intercept) 0.79967 0.19420 4.118 0.000130 *** # order.group.L 0.01344 0.43425 0.031 0.975422 # order.group.Q 2.13519 0.43425 4.917 8.32e-06 *** # order.group.C 0.11015 0.43425 0.254 0.800703 # order.group^4 -0.79602 0.43425 -1.833 0.072202 . # The t value for linear is same as I got (as are others) but I don't understand # the estimates. The intercept is the grand mean, but I don't see the relationship # of other estimates to that or to the ones I get by hand. # My estimates are the sum of (coeff times means) i.e. 0 (intercept), .0425, 7.989, .3483, -6.66 # and these are not a linear (or other nice pretty) function of est. from R. # OK, let's break it down Means - tapply(dv, order.group,
[R] (no subject)
Hello, I am having a problem figuring out how to model a continuous outcome (y) given a continuous predictor (x1) and two levels of nested categorical predictors (x3 nested in x2). The data are observational, not from a designed experiment. There are about 15 levels of x2 and between 3 and 14 levels of x3 nested within each level of x2. There are between 6 and 50 x1, y observations for each unique x3 (i.e. the data are unbalanced). I would like to get a prediction equation for y based on x1 and the levels of x2 and x3, and be able to test for significant effects of the levels of x2 and x3. The variables x2 and x3 are drawn from a fixed population. [[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] nested unbalanced regression analysis
Hello, I am having a problem figuring out how to model a continuous outcome (y) given a continuous predictor (x1) and two levels of nested categorical predictors (x3 nested in x2). The data are observational, not from a designed experiment. There are about 15 levels of x2 and between 3 and 14 levels of x3 nested within each level of x2. There are between 6 and 50 x1,y observations for each unique x3 (i.e. the data are unbalanced) . I would like to get a prediction equation for y based on x1 and the levels of x2 and x3, and be able to test for significant effects of the levels of x2 and x3. The variables x2 and x3 are drawn from a fixed population. [[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] how to avoid NaN in optim()
Change the `else NA' to `else Inf' in your loglik function and then run the
following:
library(BB)
p0 - runif(2)
spg(p0, lik (176,182 , 60 ,17) , lower=0, upper=1)
This will give you correct results.
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: arindam fadikar arindam.fadi...@gmail.com
Date: Thursday, September 30, 2010 2:17 pm
Subject: [R] how to avoid NaN in optim()
To: r-help@r-project.org
hi ,
lik - function(nO, nA, nB, nAB){
loglik - function(par)
{
p=par[1]
q=par[2]
r - 1 - p - q
if (c(p,q,r) rep(0,3) c(p,q,r) rep(1,3) )
{
-(2 * nO * log (r) + nA * log (p^2 + 2 * p * r)
+ nB * log (q^2 + 2 * q * r)
+ nAB * (log(2) +log(p) +log(q)))
}
else
NA
}
loglik
}
i want to maximize this likelihood function over the range (0,0,0) to
(1,1,1). so i tried
optim ( c(0.3,0.3), lik ( 176,182 , 60 ,17) , method = CG)
and obtained the following :
optim(c(0.3,0.3), fn, method=CG)
$par
[1] 0.26444187 0.09316946
$value
[1] 492.5353
$counts
function gradient
130 19
$convergence
[1] 0
$message
NULL
Warning messages:
1: In log(q^2 + 2 * q * r) : NaNs produced
2: In log(q) : NaNs produced
3: In log(q^2 + 2 * q * r) : NaNs produced
4: In log(q) : NaNs produced
please help ...
--
Arindam Fadikar
M.Stat
Indian Statistical Institute.
New Delhi, India
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
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] how to avoid NaN in optim()
I forgot to mention:
You actually got correct results with using optim and `CG'. The warning
messages were just telling you that there were some log(negative number)
operations during the iterative process.
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: arindam fadikar arindam.fadi...@gmail.com
Date: Thursday, September 30, 2010 2:17 pm
Subject: [R] how to avoid NaN in optim()
To: r-help@r-project.org
hi ,
lik - function(nO, nA, nB, nAB){
loglik - function(par)
{
p=par[1]
q=par[2]
r - 1 - p - q
if (c(p,q,r) rep(0,3) c(p,q,r) rep(1,3) )
{
-(2 * nO * log (r) + nA * log (p^2 + 2 * p * r)
+ nB * log (q^2 + 2 * q * r)
+ nAB * (log(2) +log(p) +log(q)))
}
else
NA
}
loglik
}
i want to maximize this likelihood function over the range (0,0,0) to
(1,1,1). so i tried
optim ( c(0.3,0.3), lik ( 176,182 , 60 ,17) , method = CG)
and obtained the following :
optim(c(0.3,0.3), fn, method=CG)
$par
[1] 0.26444187 0.09316946
$value
[1] 492.5353
$counts
function gradient
130 19
$convergence
[1] 0
$message
NULL
Warning messages:
1: In log(q^2 + 2 * q * r) : NaNs produced
2: In log(q) : NaNs produced
3: In log(q^2 + 2 * q * r) : NaNs produced
4: In log(q) : NaNs produced
please help ...
--
Arindam Fadikar
M.Stat
Indian Statistical Institute.
New Delhi, India
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
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] how to avoid NaN in optim()
You also need the constrain that par[1] + par[2] 1 in order to avoid NaNs.
You can do this using the `projectLinear' argument in `spg'.
library(BB)
?spg
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Ravi Varadhan rvarad...@jhmi.edu
Date: Thursday, September 30, 2010 2:54 pm
Subject: Re: [R] how to avoid NaN in optim()
To: arindam fadikar arindam.fadi...@gmail.com
Cc: r-help@r-project.org
Change the `else NA' to `else Inf' in your loglik function and then
run the following:
library(BB)
p0 - runif(2)
spg(p0, lik (176,182 , 60 ,17) , lower=0, upper=1)
This will give you correct results.
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: arindam fadikar arindam.fadi...@gmail.com
Date: Thursday, September 30, 2010 2:17 pm
Subject: [R] how to avoid NaN in optim()
To: r-help@r-project.org
hi ,
lik - function(nO, nA, nB, nAB){
loglik - function(par)
{
p=par[1]
q=par[2]
r - 1 - p - q
if (c(p,q,r) rep(0,3) c(p,q,r) rep(1,3) )
{
-(2 * nO * log (r) + nA * log (p^2 + 2 * p * r)
+ nB * log (q^2 + 2 * q * r)
+ nAB * (log(2) +log(p) +log(q)))
}
else
NA
}
loglik
}
i want to maximize this likelihood function over the range (0,0,0)
to
(1,1,1). so i tried
optim ( c(0.3,0.3), lik ( 176,182 , 60 ,17) , method = CG)
and obtained the following :
optim(c(0.3,0.3), fn, method=CG)
$par
[1] 0.26444187 0.09316946
$value
[1] 492.5353
$counts
function gradient
130 19
$convergence
[1] 0
$message
NULL
Warning messages:
1: In log(q^2 + 2 * q * r) : NaNs produced
2: In log(q) : NaNs produced
3: In log(q^2 + 2 * q * r) : NaNs produced
4: In log(q) : NaNs produced
please help ...
--
Arindam Fadikar
M.Stat
Indian Statistical Institute.
New Delhi, India
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
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] how to avoid NaN in optim()
thanks Ravi yes , we were getiing the correct results but we thought
there should be a way to avoid those NaN values ... and now we are done
the if condition was written wrongly there ... instead of that if we do
p 0 q 0 r 0 p 1 q 1 r 1
its perfectly fine .. but is there a smart way to write this condition ?
On Fri, Oct 1, 2010 at 12:30 AM, Ravi Varadhan rvarad...@jhmi.edu wrote:
I forgot to mention:
You actually got correct results with using optim and `CG'. The warning
messages were just telling you that there were some log(negative number)
operations during the iterative process.
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: arindam fadikar arindam.fadi...@gmail.com
Date: Thursday, September 30, 2010 2:17 pm
Subject: [R] how to avoid NaN in optim()
To: r-help@r-project.org
hi ,
lik - function(nO, nA, nB, nAB){
loglik - function(par)
{
p=par[1]
q=par[2]
r - 1 - p - q
if (c(p,q,r) rep(0,3) c(p,q,r) rep(1,3) )
{
-(2 * nO * log (r) + nA * log (p^2 + 2 * p * r)
+ nB * log (q^2 + 2 * q * r)
+ nAB * (log(2) +log(p) +log(q)))
}
else
NA
}
loglik
}
i want to maximize this likelihood function over the range (0,0,0) to
(1,1,1). so i tried
optim ( c(0.3,0.3), lik ( 176,182 , 60 ,17) , method = CG)
and obtained the following :
optim(c(0.3,0.3), fn, method=CG)
$par
[1] 0.26444187 0.09316946
$value
[1] 492.5353
$counts
function gradient
130 19
$convergence
[1] 0
$message
NULL
Warning messages:
1: In log(q^2 + 2 * q * r) : NaNs produced
2: In log(q) : NaNs produced
3: In log(q^2 + 2 * q * r) : NaNs produced
4: In log(q) : NaNs produced
please help ...
--
Arindam Fadikar
M.Stat
Indian Statistical Institute.
New Delhi, India
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
and provide commented, minimal, self-contained, reproducible code.
--
Arindam Fadikar
M.Stat
Indian Statistical Institute.
New Delhi, India
[[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] how to avoid NaN in optim()
Here is how you do it:
library(BB)
Amat - matrix(c(1,0,0,1,-1,-1), 3, 2, byrow=TRUE)
b - c(0, 0, -1)
p0 - c(0.5, 0.4)
spg(p0, lik ( 176,182 , 60 ,17) , project=projectLinear,
projectArgs=list(A=Amat, b=b, meq=0))
Hope this helps,
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Ravi Varadhan rvarad...@jhmi.edu
Date: Thursday, September 30, 2010 3:04 pm
Subject: Re: [R] how to avoid NaN in optim()
To: Ravi Varadhan rvarad...@jhmi.edu
Cc: arindam fadikar arindam.fadi...@gmail.com, r-help@r-project.org
You also need the constrain that par[1] + par[2] 1 in order to avoid
NaNs.
You can do this using the `projectLinear' argument in `spg'.
library(BB)
?spg
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: Ravi Varadhan rvarad...@jhmi.edu
Date: Thursday, September 30, 2010 2:54 pm
Subject: Re: [R] how to avoid NaN in optim()
To: arindam fadikar arindam.fadi...@gmail.com
Cc: r-help@r-project.org
Change the `else NA' to `else Inf' in your loglik function and
then
run the following:
library(BB)
p0 - runif(2)
spg(p0, lik (176,182 , 60 ,17) , lower=0, upper=1)
This will give you correct results.
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: arindam fadikar arindam.fadi...@gmail.com
Date: Thursday, September 30, 2010 2:17 pm
Subject: [R] how to avoid NaN in optim()
To: r-help@r-project.org
hi ,
lik - function(nO, nA, nB, nAB){
loglik - function(par)
{
p=par[1]
q=par[2]
r - 1 - p - q
if (c(p,q,r) rep(0,3) c(p,q,r) rep(1,3) )
{
-(2 * nO * log (r) + nA * log (p^2 + 2 * p * r)
+ nB * log (q^2 + 2 * q * r)
+ nAB * (log(2) +log(p) +log(q)))
}
else
NA
}
loglik
}
i want to maximize this likelihood function over the range
(0,0,0)
to
(1,1,1). so i tried
optim ( c(0.3,0.3), lik ( 176,182 , 60 ,17) , method = CG)
and obtained the following :
optim(c(0.3,0.3), fn, method=CG)
$par
[1] 0.26444187 0.09316946
$value
[1] 492.5353
$counts
function gradient
130 19
$convergence
[1] 0
$message
NULL
Warning messages:
1: In log(q^2 + 2 * q * r) : NaNs produced
2: In log(q) : NaNs produced
3: In log(q^2 + 2 * q * r) : NaNs produced
4: In log(q) : NaNs produced
please help ...
--
Arindam Fadikar
M.Stat
Indian Statistical Institute.
New Delhi, India
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
PLEASE do read the posting guide
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] how to avoid NaN in optim()
On Thu, Sep 30, 2010 at 12:04 PM, arindam fadikar
arindam.fadi...@gmail.com wrote:
thanks Ravi yes , we were getiing the correct results but we thought
there should be a way to avoid those NaN values ... and now we are done
the if condition was written wrongly there ... instead of that if we do
p 0 q 0 r 0 p 1 q 1 r 1
its perfectly fine .. but is there a smart way to write this condition ?
combine p, q, and r. Here imagine x = c(p, q, r)
x - c(.5, .5, .5)
all(x 0 x 1)
x2 - c(.5, .5, 1.01)
all(x2 0 x2 1)
In this case you will not want to use because x has more than one
element. all() basically just checks that everything is TRUE.
Josh
On Fri, Oct 1, 2010 at 12:30 AM, Ravi Varadhan rvarad...@jhmi.edu wrote:
I forgot to mention:
You actually got correct results with using optim and `CG'. The warning
messages were just telling you that there were some log(negative number)
operations during the iterative process.
Ravi.
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
- Original Message -
From: arindam fadikar arindam.fadi...@gmail.com
Date: Thursday, September 30, 2010 2:17 pm
Subject: [R] how to avoid NaN in optim()
To: r-help@r-project.org
hi ,
lik - function(nO, nA, nB, nAB){
loglik - function(par)
{
p=par[1]
q=par[2]
r - 1 - p - q
if (c(p,q,r) rep(0,3) c(p,q,r) rep(1,3) )
{
-(2 * nO * log (r) + nA * log (p^2 + 2 * p * r)
+ nB * log (q^2 + 2 * q * r)
+ nAB * (log(2) +log(p) +log(q)))
}
else
NA
}
loglik
}
i want to maximize this likelihood function over the range (0,0,0) to
(1,1,1). so i tried
optim ( c(0.3,0.3), lik ( 176,182 , 60 ,17) , method = CG)
and obtained the following :
optim(c(0.3,0.3), fn, method=CG)
$par
[1] 0.26444187 0.09316946
$value
[1] 492.5353
$counts
function gradient
130 19
$convergence
[1] 0
$message
NULL
Warning messages:
1: In log(q^2 + 2 * q * r) : NaNs produced
2: In log(q) : NaNs produced
3: In log(q^2 + 2 * q * r) : NaNs produced
4: In log(q) : NaNs produced
please help ...
--
Arindam Fadikar
M.Stat
Indian Statistical Institute.
New Delhi, India
[[alternative HTML version deleted]]
__
r-h...@r-project.org mailing list
PLEASE do read the posting guide
and provide commented, minimal, self-contained, reproducible code.
--
Arindam Fadikar
M.Stat
Indian Statistical Institute.
New Delhi, India
[[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.
--
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.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] how to avoid NaN in optim()
arindam fadikar wrote:
loglik - function(par)
{
p=par[1]
q=par[2]
r - 1 - p - q
if (c(p,q,r) rep(0,3) c(p,q,r) rep(1,3) )
{
-(2 * nO * log (r) + nA * log (p^2 + 2 * p * r)
+ nB * log (q^2 + 2 * q * r)
+ nAB * (log(2) +log(p) +log(q)))
}
else
NA
}
loglik
}
.
Extending the tests in the if in loglik to
if (c(p,q,r) rep(0,3) c(p,q,r) rep(1,3) (p^2 + 2*p*r)0
(q^2 + 2*q*r)0)
would also help.
/Berend
--
View this message in context:
http://r.789695.n4.nabble.com/how-to-avoid-NaN-in-optim-tp2738093p2746635.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] Nested unbalanced regression analysis
Hello, I am having a problem figuring out how to model a continuous outcome (y) given a continuous predictor (x1) and two levels of nested categorical predictors (x3 nested in x2). The data are observational, not from a designed experiment. There are about 15 levels of x2 and between 3 and 14 levels of x3 nested within each level of x2. There are between 6 and 50 x1 and y observations for each unique x3 (i.e. the data are unbalanced). I would like to get a prediction equation for y based on x1 and the levels of x2 and x3, and be able to test for significant effects of the levels of x2 and x3. The variables x2 and x3 are drawn from a fixed population. I deeply appreciate your assistance if available to assist. R/Bob [[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] Sweave and LaTeX beamer class
I am failing to uncover Sweave chunks step by step using the LaTeX beamer
class.
The following minimal example:
\documentclass{beamer}
\usepackage{Sweave}
\begin{document}
\begin{frame}[fragile]
In the year \uncover2-{25}\uncover3-{\Sexpr{5*5}}
\uncover4-{
echo=TRUE, print=TRUE=
5*5*101
@
}
\end{frame}
\end{document}
leads to an error message when running pdflatex over the *.tex file:
[...]
! FancyVerb Error:
Extraneous input ` 2525 \end {Soutput} \end {Schunk} \bea...@endcovered ' bet
ween \begin{Soutput}[key=value] and line end
.
\...@error ... {FancyVerb Error:
\space \space #1
}
--
Johannes Hüsing There is something fascinating about science.
One gets such wholesale returns of conjecture
mailto:johan...@huesing.name from such a trifling investment of fact.
http://derwisch.wikidot.com (Mark Twain, Life on the Mississippi)
__
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] Problem comparing lme objects with anova using do.call
Hi all, I am running some linear mixed models for longitudinal data using the nlme package. Part of the code I've developed tests for differences among models using the anova function, but I can't simply pass the corresponding lme objects as parameters to the anova function because some of the models do not converge and others are retained for comparison only if certain criteria are met. So, the simple solution seemed to be generating a list of models and using do.call to test models that made the cut, either by converging, or by meeting some of the substantive criteria I set. Unfortunately, the do.call command generates an error. I dug around on the R-help archives and came across some similar problems, but none that helped to resolve this error. The error I receive is: Error in `row.names-.data.frame`(`*tmp*`, value = c(structure(list(modelStruct = structure(list(reStruct = structure(list(, : invalid 'row.names' length The traceback() didn't provide clear feedback about what was wrong, at least to my ignorant eye. Here is a simple reproducible example: library(nlme) dataTest - data.frame(y=c(rnorm(100, 0, 1), rnorm(100,5,1)), time=rep(0:3, each=50), id=rep(1:50, 4)) uncMeansModel - lme(y ~ 1, data=dataTest, random=~1|id, method=ML) uncGrowthModel - lme(y ~ 1 + time, data=dataTest, random=~time|id, method=ML) uncGrowthModel2 - lme(y ~ 1 + time + I(time^2), data=dataTest, random=~time|id, method=ML) models - list(uncMeansModel, uncGrowthModel, uncGrowthModel2) do.call(anova.lme, models) #doesn't work anova.lme(uncMeansModel, uncGrowthModel, uncGrowthModel2) #works And here is my sessionInfo(): sessionInfo() R version 2.11.1 (2010-05-31) x86_64-pc-linux-gnu locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=en_US.UTF-8 [9] LC_ADDRESS=en_US.UTF-8LC_TELEPHONE=en_US.UTF-8 [11] LC_MEASUREMENT=en_US.UTF-8LC_IDENTIFICATION=en_US.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] nlme_3.1-96 rj_0.5.0-5 loaded via a namespace (and not attached): [1] grid_2.11.1 lattice_0.19-11 rJava_0.8-6 tools_2.11.1 Thanks for any help you can provide! Best wishes, Michael Hallquist [[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] interactive session
Thanks Niels but it won't do.. please copy and paste the 2 lines below together to your console in order to see what I mean: cat(?); a-readLines(n=1) b-paste(t,a,sep=) anyone / any idea to overcome this problem? Best, Fatih Niels wrote: Hi Fatih I believe that readLines(n=1) will do the job. It works fine from the Windows RGui, but I noticed that it hangs my Aquamacs/ESS when R runs from there, and a C-g was needed (which may be completely irrelevant to you). Best, Niels On 30/09/10 08.55, Pam wrote: Hi guys, My concern is to create an automated process from the beginning to the end. I want to copy all my code together in one piece and paste it to R console and sit back and relax :) except one moment in which the program should ask me to enter a number.. and only then (only after getting a valid number from me) it should continue to read and process the rest of the code. I�tried�lots of things (readline, readLines, scan, interactive, ask, switch,...) and read manuals and searched help forums.. I found several similar questions but never a satisfying answer.. so now it became a challenge.. any idea how�to�overcome this problem (R 2.11.1 for Windows)? (an example�code is below)� Best, Fatih [[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] Kernel density estimation - problem with matrix created in loop
Hi you all!
I have a very serious problem.
What I am going to do is a kernel density estimation (Nadaraya-Watson) but
without any helping packages. I want to do it on my own and just use R as a
calculator.
The problem is that the matrix which I need does not seem to be calculated
correctly in my loop.
This matrix is created by the multiplication of a vector and its
transponent. But before it becomes multiplicated with a kernel term.
Just let me show this:
A=matrix(0,nrow=22,ncol=22)
temp=matrix(0,nrow=22,ncol=22)
for (i in 1:1473)
{
vectora=c(X_1[i],X_2[i],X_3[i],X_4[i],X_5[i],X_6[i],X_7[i],X_8[i],X_9[i],X_10[i],X_11[i],X_12[i],X_13[i],X_14[i],X_15[i],X_16[i],X_17[i],X_18[i],X_19[i],X_20[i],X_21[i],X_22[i])
vectorb=t(c(X_1[i],X_2[i],X_3[i],X_4[i],X_5[i],X_6[i],X_7[i],X_8[i],X_9[i],X_10[i],X_11[i],X_12[i],X_13[i],X_14[i],X_15[i],X_16[i],X_17[i],X_18[i],X_19[i],X_20[i],X_21[i],X_22[i]))
temp=(3/(4*h))*(1-((i/T - u)/h)^2)*vectora%*%vectorb
Anw=A+temp
}
This is the problem part, the part with the matrix. You see that it is a
22x22 matrix. Each vector has got more than 1000 elements. Thus there are
more than one thousand of those matrices. And temp is the command in which
each of those matrices should be multiplicated with the i-corresponding
calue of (3/(4*h))*(1-((i/T - u)/h)^2) [Epanechnikov kernel].
And with Anw=A+temp finally there should be only a sum. And it is a sum, a
22x22matrix. But I cannot continue with my work because in the formula I
have to use this I have to invert this matrix before. And this does not
work. The determinant is said as being 0.
I am sure that I am wrong with something. Can someone help me??
--
View this message in context:
http://r.789695.n4.nabble.com/Kernel-density-estimation-problem-with-matrix-created-in-loop-tp2720793p2720793.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] R 2.11.1 crashes
Dear R community, I was using R 2.11.1 without internet connection on Windows XP and whenever I type ?mean for example, R would freeze and crash out... Is this something that can be fixed? I would like to use the internel help file if possible... Thanks. ### Assistant Professor Steve Su School of Mathematics and Statistics Faculty of Engineering, Computing and Mathematics M019, 35 Stirling Highway Crawley, 6009, WA, Australia Phone:+6164883369 http://www.uwa.edu.au/people/steve.su CRICOS Provider Code: 00126G [[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] Inserting a plot into another
Hi everyone. I would like to know if it was possible to insert a plot into another one. For example I have a plot and I would like to add a smaller plot in the top right corner. Best regards, Phil -- View this message in context: http://r.789695.n4.nabble.com/Inserting-a-plot-into-another-tp2720936p2720936.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] getting just the x or the y coordinate
I have a data set (from GRASSGIS) I have imported it to R and I can get the coordinates by: xy-coordinates(data) but I would like to have the x coordinates and y coordinates in separate dataframes x and y, does anyone know how to do this? Thankyou Gary --- Gary R. Nobles PhD Candidate Institute of Archaeology Poststraat 6 Groningen University The Netherlands Personal Webpage: Project Website: www.singlegrave.nl g.r.nob...@rug.nl __ 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] more than two NA value names in my data
my data(*.txt) has 1000 observations(numbers with no characters) of 5 variables. quite simple. However, NA values are quite tricky. this observer used more than two names for NA values; . and na and more. 1. If I don't want to manipulate this raw data at all, how can I read this table? (meaning, can I set more than two names for na.strings= in read.table()?) 2. What happens if I don't set any NA value names in read.table()? Does R read .(dot) and na as NA value? 3. If Q1 is not possible, what is the best way to get the values I want? (let's say I want means of all variables. I give conditions on my mean()?) -- View this message in context: http://r.789695.n4.nabble.com/more-than-two-NA-value-names-in-my-data-tp2730161p2730161.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] polr, lrm - ordinal data
Dear List, I have developed a model and am looking to predict a response for 1-6 ( it is ordered i.e the difference between level 1 and 2 is the same as between level 2 and 3 etc. I have used the predict function for a polr model (below) and a lrm model, and both give similar results, however for some reason the outcome are all 1 or 6 i.e no level 2,3,4,5. This is not correct, i am unsure if it is because my model is not good enough to predict to this accuracy or if it is something i am doing wrong? Thanks Chris test - polr(extinction ~ FR*HAB+WO+ALT+BIO+REG,method=logistic) test Call: polr(formula = extinction ~ FR * HAB + WO + ALT + BIO + REG, method = logistic) Coefficients: FRNon_fleshy HABSemi-aquatic HABTerrestrial WOWoody ALTHigh 0.09758543 -0.05988101 -0.29744997 0.32746840 -0.39191606 ALTLow ALTMid BIOBoreal BIOMediterranean BIOSubtropical/Tropical -0.56523156 -0.18979562 -0.22743656 -0.31344233 1.77031824 BIOTemperate REGTwo_plus FRNon_fleshy:HABSemi-aquatic FRNon_fleshy:HABTerrestrial 1.45071627 -0.67654880 -2.06919408 -0.31706541 Intercepts: 1|2 2|3 3|4 4|5 5|6 0.4874828 0.8340901 1.3994091 1.8091463 2.1295630 Residual Deviance: 2471.053 AIC: 2509.053 predict(test) [1] 1 1 1 1 1 1 1 1 1 1 6 1 6 1 1 6 1 1 1 6 1 1 1 6 1 1 1 1 6 1 6 6 1 1 6 1 1 1 1 6 1 1 1 1 6 6 1 1 6 1 1 1 1 1 6 1 1 1 1 6 1 1 6 1 6 1 6 1 1 1 1 1 1 1 [75] 1 1 1 1 1 1 1 6 6 1 1 1 6 6 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 6 6 1 1 1 1 1 1 6 1 1 1 1 1 6 1 1 1 1 1 6 6 1 6 1 1 1 1 1 1 6 1 1 1 1 1 6 1 1 1 1 6 [149] 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 6 1 1 1 1 1 1 6 1 6 1 1 1 1 1 6 1 1 6 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 6 1 1 6 1 1 1 1 1 1 6 1 1 [223] 1 1 1 1 1 1 1 1 1 1 1 6 1 1 6 6 6 1 1 6 6 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 6 1 1 1 1 1 1 6 1 1 1 1 6 1 1 6 6 1 6 1 1 1 1 1 1 1 1 1 1 1 1 6 1 6 6 1 1 [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 6 1 1 6 1 1 6 1 6 1 1 6 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 6 1 1 [371] 1 1 1 1 1 1 1 6 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 6 6 6 1 1 1 1 6 1 1 1 1 1 1 1 1 6 1 1 1 6 1 1 1 1 1 1 1 6 1 1 1 1 6 1 6 1 1 6 1 1 1 1 1 1 1 1 1 1 1 [445] 1 6 1 1 1 1 1 1 6 1 1 1 1 1 6 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 [519] 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 6 1 1 1 1 1 1 1 6 1 1 6 1 6 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 6 1 1 6 1 6 1 1 6 6 1 1 1 1 1 1 1 1 1 1 1 [593] 6 1 1 1 1 1 1 6 1 1 6 1 1 6 6 1 6 1 1 1 6 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 6 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 [667] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 1 1 1 1 [741] 1 1 1 6 6 1 1 1 6 1 1 1 1 1 6 1 1 1 1 1 1 1 6 1 1 6 1 1 6 1 1 1 6 1 6 1 6 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 6 1 1 1 1 1 6 1 [815] 6 6 1 1 1 1 1 1 1 6 1 6 1 1 1 1 1 1 1 1 Levels: 1 2 3 4 5 6 this is the correct output extinction [1] 3 1 6 4 2 1 6 1 2 6 6 1 5 1 1 1 1 6 2 4 1 4 1 3 6 3 5 1 4 3 3 5 3 3 3 5 1 6 6 3 6 6 1 1 5 1 5 3 3 6 1 1 6 2 6 1 3 1 2 6 1 1 3 4 6 3 3 1 1 3 6 3 1 6 [75] 3 4 1 1 1 3 1 3 6 1 2 1 6 6 3 6 6 5 1 1 4 6 3 1 2 1 1 1 4 6 3 1 2 6 1 3 2 1 1 6 1 6 2 3 6 2 4 6 6 5 2 4 6 6 1 2 6 6 4 1 1 1 5 5 1 1 1 6 6 1 1 1 2 6 [149] 4 5 6 6 6 1 3 1 4 3 1 3 6 1 3 6 1 2 1 3 4 1 3 1 6 3 6 1 6 2 1 3 4 6 6 1 1 1 1 5 1 2 1 6 3 1 6 6 1 1 6 6 1 1 2 2 4 1 4 3 5 4 1 4 3 6 1 2 6 1 1 5 6 1 [223] 1 2 1 6 3 6 1 1 1 6 1 6 3 4 1 5 5 6 3 2 3 6 1 1 1 5 1 6 6 1 4 3 1 4 6 1 3 5 1 4 4 6 2 6 6 1 6 3 2 6 3 1 4 1 5 6 1 6 4 5 4 1 4 1 4 1 1 1 4 1 6 6 6 6 [297] 1 1 1 1 6 1 3 2 1 6 6 6 1 1 3 3 1 4 1 6 2 2 6 3 1 1 6 6 2 6 1 1 6 6 1 2 1 6 4 1 4 1 1 5 1 3 6 6 1 1 3 6 6 6 1 6 1 1 6 3 4 6 1 6 1 3 1 2 1 1 6 6 1 1 [371] 1 6 3 1 1 2 5 1 1 1 1 1 6 1 1 1 1 2 4 6 1 5 1 1 6 6 3 1 6 1 6 6 4 1 1 4 2 6 1 1 3 6 1 2 6 1 2 1 1 4 5 1 6 1 6 1 1 6 1 3 6 1 1 1 3 1 2 6 3 5 2 1 3 1 [445] 3 2 4 1 6 6 1 4 5 2 1 6 1 3 4 3 1 1 1 6 1 3 5 6 1 1 1 6 1 1 2 1 1 1 1 2 4 1 1 1 1 3 1 1 1 1 1 1 6 6 4 1 1 2 1 6 6 6 6 6 6 6 1 3 4 3 1 1 1 1 2 1 5 1 [519] 1 6 1 1 1 6 6 1 3 5 2 1 1 1 1 4 3 6 2 1 3 1 1 1 6 3 3 4 1 5 6 6 1 2 6 3 1 5 6 1 1 6 1 1 1 1 1 3 1 1 1 3 6 1 1 1 6 6 1 3 1 2 6 1 6 2 3 1 6 1 6 6 2 4 [593] 1 1 5 3 1 1 1 6 3 1 6 6 4 4 1 6 6 6 1 1 6 1 6 3 1 1 2 1 6 1 1 1 6 1 3 6 6 1 3 3 2 2 6 1 1 1 1 6 1 3 1 2 1 6 1 4 1 1 1 6 1 4 4 1 6 3 1 6 1 3 1 1 2 6 [667] 5 1 2 1 1 4 4 1 1 1 6 4
[R] time in year, month, day, hour ?
Dear R Users, I did not get any reply on my question so I am re-asking. This time I am giving sample data: 160.3162 -13.5993 -0.4353 46.09380.1877 -0.194E-07 260.3713 -13.5992 -0.4423 46.12410.2057 -0.231E-06 360.3430 -13.5981 -1.6163 44.90480.2237 -0.270E-06 460.3227 -13.5970 -2.6258 43.87850.2213 -0.139E-06 560.3352 -13.5961 -2.5874 43.92380.2278 0.141E-06 ... . .. 1918 59.1785 -14.58510.3895 44.9850 -0.0021 0.141E-06 1919 59.1816 -14.56220.3933 44.99720.0155 0.139E-06 1920 59.1637 -14.56660.4420 45.02190.0172 0.138E-06 Column 1, is index of hours on 3 hourly interval and starts from 1 Jan 2009, 00 hrs and run till 240 days of 2009. I want to convert Column 1 of above data in the format, 'year' , 'month', 'day', 'hour' Kindly help. Thanks, Regards, Yogesh My question is below Dear R Users, I have model simulated data for 240 days. The day 1 is Jan 1, 2009, 00:00 hrs and then with 3-hourly interval and so on. The time axis is 1,2,3,4..1920; so the total rows in the data are 1920. How to convert above time axis in year month day hour format Great Thanks, regards, Yogesh -- Yogesh K. Tiwari (Dr.rer.nat), Scientist, Centre for Climate Change Research, Indian Institute of Tropical Meteorology, Homi Bhabha Road, Pashan, Pune-411008 INDIA Phone: 0091-99 2273 9513 (Cell) : 0091-20-25904452 (O) Fax: 0091-20-258 93 825 [[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] xyplot - help please
Dear R Gurus!
I would like to generate a multi-panel plot with grouped and reference lines
in each panel .. Example dataset and code are provide below. Where am I
getting wrong in the code below? Any assistance will be highly appreciated.
Thanks so much!
-Santosh
start dataset ---
X,l,y1,y2,y3,y4,f
-99,15,0.178,0.127,0.152,0.116,a
0,0,0.218,0.12,0.174,0.145,a
6,0.0015561,0.227,0.122,0.175,0.148,a
11,0.018907,0.213,0.12,0.176,0.145,a
17,0.079668,0.204,0.105,0.153,0.16,a
22,0.2368,0.216,0.121,0.175,0.148,a
27,0.58763,0.217,0.122,0.175,0.148,a
32,1.3433,0.223,0.118,0.17,0.145,a
37,2.7883,0.214,0.115,0.172,0.146,a
42,5.6097,0.22,0.121,0.175,0.147,a
47,11.018,0.216,0.12,0.176,0.145,a
-99,15,0.178,0.127,0.152,0.116,b
4,0.0015561,0.159,0.145,0.114,0.543,b
8,0.018907,0.161,0.145,0.116,0.544,b
12,0.079668,0.143,0.139,0.121,0.539,b
17,0.2368,0.161,0.145,0.114,0.549,b
20,0.58763,0.161,0.145,0.116,0.546,b
25,1.3433,0.162,0.144,0.115,0.556,b
29,2.7883,0.162,0.144,0.115,0.549,b
33,5.6097,0.099,0.1,0.302,0.799,b
37,11.018,0.161,0.144,0.114,0.557,b
end dataset ---
library(reshape)
options(scipen=12)
my.settings -
list(plot.symbol = list(col = c(red, blue), pch = 20,alpha=c(0.8,
0.8)),
plot.line = list(col = c(red, blue), lty = c(2,3)))
x1 - read.csv(ex1.csv,as.is=T)
y1 -
melt(x1,measure.vars=names(x1)[3:(length(names(x1))-1)],stringsAsFactors=F)
y1a - y1[order(y1$f,y1$l),] #
ref - y1a[y1a$X0,]
xyplot(value~l|variable,groups=f,main=Compare parameter estimates,
outer=T,allow.mult=T,
data=y1a[y1a$X0,],
strip=strip.custom(strip.names=T,strip.levels=T,var.name=Par),
prepanel=function(x,y,group.number,..., horizontal) {
rngy - range(y,finite=T)
rngx - range(x,finite=T)
list(ylim=rngy,xlim=rngx)
},
panel.superpose,
panel.xyplot,
panel.groups = function(x, y, group.number,subscripts=subscripts,
..., horizontal) {
opar - trellis.par.get()
on.exit(trellis.par.set(opar))
rv -
ref$value[ref$variable%in%sort(unique(y1$variable))[panel.number()]]
trellis.par.set(list(plot.symbol =
Rows(my.settings$plot.symbol,group.number),
plot.line = Rows(my.settings$plot.line,group.number)))
panel.xyplot(x, y, horizontal = horizontal)
panel.abline(h=rv,col='green')
},
xlab=Levels,ylab=Values,
scale=list(x=list(log=T,at=c(0.01,0.1,1,10,15),
lab=c(0.01,0.1,1,10,15),rot=45),rela=free),
auto.key=list(space=bottom,columns=2))
[[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] please help me....
ㅣ recently used laq packages author of this package is jan Ulbricht.. and I must read jan ulbricht's Ph.D thesis variable selection in generalized linear models but I didn't find this thesis on line and there is any method to view this thesis?? If anybody has this thesis file,please send to me... I must this guy's Ph.D. thesis... I searched LMU online library..but this thesis did not exist... please help me...very very very important issue for me.. mu e-mail address:butchma...@gmail.com -- View this message in context: http://r.789695.n4.nabble.com/please-help-me-tp2739350p2739350.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] more than two NA value names in my data
Hi, You were on the right track with na.strings, from ?read.table na.strings: a character vector of strings which are to be interpreted as ‘NA’ values. Blank fields are also considered to be missing values in logical, integer, numeric and complex fields. so, you can just do something like na.strings = c(., na, anotherthing) and so on. If you leave it blank, R will treat NAs as NA, but it is not going to mystically know what values are someone's special term for missing and what are real values. So it would read the data in as character or factor. You could also work to clean the data after you had read it in and then convert everything back to numeric, but it is easier to just specifying what indicates missing values. Hope that helps, Josh On Thu, Sep 30, 2010 at 10:10 AM, JoonGi joo...@hanmail.net wrote: my data(*.txt) has 1000 observations(numbers with no characters) of 5 variables. quite simple. However, NA values are quite tricky. this observer used more than two names for NA values; . and na and more. 1. If I don't want to manipulate this raw data at all, how can I read this table? (meaning, can I set more than two names for na.strings= in read.table()?) 2. What happens if I don't set any NA value names in read.table()? Does R read .(dot) and na as NA value? 3. If Q1 is not possible, what is the best way to get the values I want? (let's say I want means of all variables. I give conditions on my mean()?) -- View this message in context: http://r.789695.n4.nabble.com/more-than-two-NA-value-names-in-my-data-tp2730161p2730161.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. -- Joshua Wiley Ph.D. Student, Health Psychology University of California, Los Angeles http://www.joshuawiley.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] Understanding linear contrasts in Anova using R
On 09/30/2010 08:31 PM, Max Kuhn wrote: These two resources might also help: http://cran.r-project.org/doc/contrib/Faraway-PRA.pdf http://cran.r-project.org/web/packages/contrast/vignettes/contrast.pdf Max They're a tad long though. Let me try and say it shorter: Contrast calculation and contrast parametrization are eachother's (g-)inverses. E.g. m [,1] [,2] [,3] [1,]000 [2,]100 [3,]110 [4,]111 solve(cbind(1,m)) [,1] [,2] [,3] [,4] [1,]1000 [2,] -1100 [3,]0 -110 [4,]00 -11 So you can parametrize with successive differences using m but it is the other matrix that is used to compute successive differences. (There's quite some confusion going on caused by the obsession with ORTHOGONAL contrasts of earlier days. For those kinds of contrasts (only!), the inverse equals the transpose save for scaling factors. In the case of polynomial contrasts, they are even designed to be orthonormal so that the transpose equals the g-inverse) On Thu, Sep 30, 2010 at 1:33 PM, Ista Zahn iz...@psych.rochester.edu wrote: Hi Professor Howell, I think the issue here is simply in the assumption that the regression coefficients will always be equal to the product of the means and the contrast codes. I tend to think of regression coefficients as the quotient of the covariance of x and y divided by the variance of x, and this definition agrees with the coefficients calculated by lm(). See below for a long-winded example. On Wed, Sep 29, 2010 at 3:42 PM, David Howell david.how...@uvm.edu wrote: #I am trying to understand how R fits models for contrasts in a #simple one-way anova. This is an example, I am not stupid enough to want #to simultaneously apply all of these contrasts to real data. With a few #exceptions, the tests that I would compute by hand (or by other software) #will give the same t or F statistics. It is the contrast estimates that R produces #that I can't seem to understand. # # In searching for answers to this problem, I found a great PowerPoint slide (I think by John Fox). # The slide pointed to the coefficients, said something like these are coeff. that no one could love, and #then suggested looking at the means to understand where they came from. I have stared # and stared at his means and then my means, but can't find a relationship. # The following code and output illustrates the problem. # Various examples of Anova using R dv - c(1.28, 1.35, 3.31, 3.06, 2.59, 3.25, 2.98, 1.53, -2.68, 2.64, 1.26, 1.06, -1.18, 0.15, 1.36, 2.61, 0.66, 1.32, 0.73, -1.06, 0.24, 0.27, 0.72, 2.28, -0.41, -1.25, -1.33, -0.47, -0.60, -1.72, -1.74, -0.77, -0.41, -1.20, -0.31, -0.74, -0.45, 0.54, -0.98, 1.68, 2.25, -0.19, -0.90, 0.78, 0.05, 2.69, 0.15, 0.91, 2.01, 0.40, 2.34, -1.80, 5.00, 2.27, 6.47, 2.94, 0.47, 3.22, 0.01, -0.66) group - factor(rep(1:5, each = 12)) # Use treatment contrasts to compare each group to the first group. options(contrasts = c(contr.treatment,contr.poly)) # The default model2 - lm(dv ~ group) summary(model2) # Summary table is the same--as it should be # Intercept is Group 1 mean and other coeff. are deviations from that. # This is what I would expect. #summary(model1) # Df Sum Sq Mean Sq F valuePr(F) # group4 62.46 15.6151 6.9005 0.0001415 *** # Residuals 55 124.46 2.2629 #Coefficients: #Estimate Std. Error t value Pr(|t|) #(Intercept) 1.802500.43425 4.151 0.000116 *** #group2 -1.127500.61412 -1.836 0.071772 . #group3 -2.715000.61412 -4.421 4.67e-05 *** #group4 -1.258330.61412 -2.049 0.045245 * #group5 0.086670.61412 0.141 0.888288 # Use sum contrasts to compare each group against grand mean. options(contrasts = c(contr.sum,contr.poly)) model3 - lm(dv ~ group) summary(model3) # Again, this is as expected. Intercept is grand mean and others are deviatoions from that. #Coefficients: # Estimate Std. Error t value Pr(|t|) # (Intercept) 0.7997 0.1942 4.118 0.000130 *** # group11.0028 0.3884 2.582 0.012519 * # group2 -0.1247 0.3884 -0.321 0.749449 # group3 -1.7122 0.3884 -4.408 4.88e-05 *** # group4 -0.2555 0.3884 -0.658 0.513399 #SO FAR, SO GOOD # IF I wanted polynomial contrasts BY HAND I would use #a(i) = -2 -1 0 1 2 for linear contrast(or some linear function of this ) #Effect = Sum(a(j)M(i))# where M = mean #Effect(linear) = -2(1.805) -1(0.675) +0(-.912) +1(.544) +2(1.889) = 0.043 #SS(linear) = n*(Effect(linear)^2)/Sum((a(j)^2)) = 12(.043)/10 = .002 #F(linear) = SS(linear)/MS(error) = .002/2.263 = .001 #t(linear) = sqrt(.001) = .031 # To do this in R I would use order.group -
Re: [R] time in year, month, day, hour ?
Dear Yogesh, This will create a vector that I believe does what you want. x - as.POSIXct(x = cumsum(c(0, rep(3*60*60, 1919))), origin = 2009-01-01) Let me see if I can explain the logic. In the innermost part, I multiply 3*60*60, or hours*minutes*seconds. You said they were three hour blocks so that is why I used three. This is because POSIXct, which I think is the format you will want to use, measures things in number of seconds since the origin (baseline point). Next, I repeat that number of seconds 1,919 times, and combine it with a 0 at the beginning. This creates a vector of length 1920 with 0, and then number of seconds in 3 hours. Next find the cumulative sum (this gets time moving forwards). Finally specify the origin (traditionally in R 1970-01-01), but we know you want it to be 2009-01-01. The results from this are assigned to x. Hope that helps, Josh On Thu, Sep 30, 2010 at 8:59 AM, Yogesh Tiwari yogesh@googlemail.com wrote: Dear R Users, I did not get any reply on my question so I am re-asking. This time I am giving sample data: 1 60.3162 -13.5993 -0.4353 46.0938 0.1877 -0.194E-07 2 60.3713 -13.5992 -0.4423 46.1241 0.2057 -0.231E-06 3 60.3430 -13.5981 -1.6163 44.9048 0.2237 -0.270E-06 4 60.3227 -13.5970 -2.6258 43.8785 0.2213 -0.139E-06 5 60.3352 -13.5961 -2.5874 43.9238 0.2278 0.141E-06 ... . .. 1918 59.1785 -14.5851 0.3895 44.9850 -0.0021 0.141E-06 1919 59.1816 -14.5622 0.3933 44.9972 0.0155 0.139E-06 1920 59.1637 -14.5666 0.4420 45.0219 0.0172 0.138E-06 Column 1, is index of hours on 3 hourly interval and starts from 1 Jan 2009, 00 hrs and run till 240 days of 2009. I want to convert Column 1 of above data in the format, 'year' , 'month', 'day', 'hour' Kindly help. Thanks, Regards, Yogesh My question is below Dear R Users, I have model simulated data for 240 days. The day 1 is Jan 1, 2009, 00:00 hrs and then with 3-hourly interval and so on. The time axis is 1,2,3,4..1920; so the total rows in the data are 1920. How to convert above time axis in year month day hour format Great Thanks, regards, Yogesh -- Yogesh K. Tiwari (Dr.rer.nat), Scientist, Centre for Climate Change Research, Indian Institute of Tropical Meteorology, Homi Bhabha Road, Pashan, Pune-411008 INDIA Phone: 0091-99 2273 9513 (Cell) : 0091-20-25904452 (O) Fax : 0091-20-258 93 825 [[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. -- Joshua Wiley Ph.D. Student, Health Psychology University of California, Los Angeles http://www.joshuawiley.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] Can this code be written more efficiently?
Have you tried using Rprof to determine where time is being spent in
the current code? Have you looked at how much memory you are using?
Are you paging? Have you run with a size 'x', then '2x' then '4x' to
see what the growth in both CPU time and memory usage is? This is
what I would do if I were trying to debug/optimize one of my scripts.
Before I would run something for a day, I would understand how the
processing time increases with the size of the input file so that I
would have an idea of how long to wait.
On Thu, Sep 30, 2010 at 1:40 PM, Guelman, Leo leo.guel...@rbc.com wrote:
Dear users,
I'm working on binary classification problem using Support Vector
Machines (SVM). My objective is to train a series of SVM models on a
grid of hyperparameters and then select those that maximize the AUC
based on an independent validation sample.
My attempted code is shown below. It runs well on small data sets but
when I use it on a slightly larger sample (e.g., my train data is
composed of about 8,000 observations on each class and 21 inputs), it
takes forever to run (more than 1 day already and still running). I'm
wondering if there's any way I can optimize this code. Thanks in advance
for any help.
I'm using 64-bit R 2.11.1 on Win 7.
Start Code
library(e1071)
library(ROCR)
### Create grid of hyperparameters
Gseq - seq(-15,3,2); G - rep(2, length(Gseq)); G - G^Gseq
Cseq - seq(-5,13,2); C - rep(2, length(Cseq)); C - C^Cseq
mygrid - expand.grid(C=C, G=G)
### Train models
svm.models - lapply(1:nrow(mygrid), function(i) {
svm(churn.form, data = mytraindata,
method = C-classification, kernel = radial,
cost = mygrid[i,1], gamma = mygrid[i,2],
probability=TRUE)
})
### Predict on test set
pred.step3 - numeric(length(svm.models))
for (i in 1:length(svm.models)) {
pred.step1 - predict(svm.models[[i]], myvaliddata, decision.values = F,
probability=T)
pred.step2 -
prediction(predictions=attr(pred.step1,probabilities)[,1],
labels=myvaliddata$churn)
pred.step3[i] - performance(pred.step2, auc)@y.values[[1]]
}
pred.step3
End Code
Thanks,
Leo.
___
This e-mail may be privileged and/or confidential, and the sender does not
waive
any related rights and obligations. Any distribution, use or copying of this
e-mail or the information
it contains by other than an intended recipient is unauthorized.
If you received this e-mail in error, please advise me (by return e-mail or
otherwise) immediately.
Ce courriel peut contenir des renseignements protégés et confidentiels.
L’expéditeur ne renonce pas aux droits et obligations qui s’y rapportent.
Toute diffusion, utilisation ou copie de ce courriel ou des renseignements
qu’il contient
par une personne autre que le destinataire désigné est interdite.
Si vous recevez ce courriel par erreur, veuillez m’en aviser immédiatement,
par retour de courriel ou par un autre moyen.
[[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
Cincinnati, OH
+1 513 646 9390
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.
Re: [R] Nested unbalanced regression analysis
Whoa, easy there: members of the list are volunteers, allow some time before expecting a reply... On Sep 30, 2010, at 3:19 PM, Robert Quinn wrote: Hello, I am having a problem figuring out how to model a continuous outcome (y) given a continuous predictor (x1) and two levels of nested categorical predictors (x3 nested in x2). The data are observational, not from a designed experiment. There are about 15 levels of x2 and between 3 and 14 levels of x3 nested within each level of x2. There are between 6 and 50 x1 and y observations for each unique x3 (i.e. the data are unbalanced). I would like to get a prediction equation for y based on x1 and the levels of x2 and x3, and be able to test for significant effects of the levels of x2 and x3. The variables x2 and x3 are drawn from a fixed population. I deeply appreciate your assistance if available to assist. R/Bob [[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] interactive session
On 09/30/2010 03:33 PM, Pam wrote:
Thanks Niels but it won't do.. please copy and paste the 2 lines below
together
to your console in order to see what I mean:
cat(?); a-readLines(n=1)
b-paste(t,a,sep=)
anyone / any idea to overcome this problem?
Best,
Fatih
You might want to source() a file with those lines rather than pasting
them to the console. There's just no way you can retroactively insert
text between two already submitted lines.
You can do things like this, though
{cat(?); a-readLines(n=1)
hey
b-paste(t,a,sep=)}
However, there's a catch
{cat(?); a-readLines(n=1)
+ hey
+ b-paste(t,a,sep=)}
?ada
b
[1] tada
Notice that the hey doesn't print.
--
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Phone: (+45)38153501
Email: pd@cbs.dk Priv: pda...@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.
[R] can I add line breaks to the paste() function?
Can I add a line break to the paste() function to return the following: 'this is the first line' 'this is the second line' instead of 'this is the first line this is the second line' ? __ 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] interactive session
On 30/09/10 22.23, Peter Dalgaard wrote:
On 09/30/2010 03:33 PM, Pam wrote:
Thanks Niels but it won't do.. please copy and paste the 2 lines below
together
to your console in order to see what I mean:
cat(?); a-readLines(n=1)
b-paste(t,a,sep=)
anyone / any idea to overcome this problem?
Best,
Fatih
You might want to source() a file with those lines rather than pasting
them to the console. There's just no way you can retroactively insert
text between two already submitted lines.
Exactly, it works when you source() the file.
- Niels
You can do things like this, though
{cat(?); a-readLines(n=1)
hey
b-paste(t,a,sep=)}
However, there's a catch
{cat(?); a-readLines(n=1)
+ hey
+ b-paste(t,a,sep=)}
?ada
b
[1] tada
Notice that the hey doesn't print.
--
Niels Richard Hansen Web: www.math.ku.dk/~richard
Associate Professor Email: niels.r.han...@math.ku.dk
Department of Mathematical Sciences nielsrichardhan...@gmail.com
University of Copenhagen Skype: nielsrichardhansen.dk
Universitetsparken 5 Phone: +45 353 20783 (office)
2100 Copenhagen Ø +45 2859 0765 (mobile)
Denmark
__
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] can I add line breaks to the paste() function?
Hello, Although on the surface a simple request, I think you need to be more specific. ?paste returns a character vector. The first question is: do you want a character vector of length 1 or 2? It sounds like you're trying to format text for display on screen or on a graphics device. Perhaps you're hoping to use ?cat instead with the sep argument set to \n ? David LeBauer wrote: Can I add a line break to the paste() function to return the following: 'this is the first line' 'this is the second line' instead of 'this is the first line this is the second line' ? __ 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] interactive session
You might want to use the tcltk package so you can bring up a window
in which to input your data. This is better than trying to read from
the console especially when 'sourcing' in some files or when
cutting/pasting.
On Thu, Sep 30, 2010 at 2:55 AM, Pam fkira...@yahoo.com wrote:
Hi guys,
My concern is to create an automated process from the beginning to the end. I
want to copy all my code together in one piece and paste it to R console and
sit
back and relax :) except one moment in which the program should ask me to
enter
a number.. and only then (only after getting a valid number from me) it should
continue to read and process the rest of the code. I tried lots of things
(readline, readLines, scan, interactive, ask, switch,...) and read manuals and
searched help forums.. I found several similar questions but never a
satisfying
answer.. so now it became a challenge.. any idea how to overcome this problem
(R
2.11.1 for Windows)? (an example code is below)
Best,
Fatih
library(gtools)
library(YaleToolkit)
library(xts)
### start of my wrong function!
f-function(w){
w-readline(which data? )
w-as.numeric(w)
ifelse(is.numeric(w)==TRUE, w, f())
}
f()
# end of my wrong function
v- ## and output of my function should be a v for example which I can use
it
in the next line (v-w or something like that??)
##the rest works fine
p-paste(t, v, .txt, sep = )
t-read.table(p, header=FALSE, sep=\t, dec=,,
blank.lines.skip=FALSE)
rownames(t)-as.Date(t[,1],%d.%m.%Y)
colnames(t)-c(date,start,high,low,end,w.average,lot,
volume)
x-as.xts(t)
whatis(x)
.
.
[[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
Cincinnati, OH
+1 513 646 9390
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.
Re: [R] can I add line breaks to the paste() function?
Try using cat instead. Then \n is the new line character. E.g. cat(1st line\n2nd line\n) Jeremy On 30 September 2010 13:30, David LeBauer dleba...@gmail.com wrote: Can I add a line break to the paste() function to return the following: 'this is the first line' 'this is the second line' instead of 'this is the first line this is the second line' ? __ 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. -- Jeremy Miles Psychology Research Methods Wiki: www.researchmethodsinpsychology.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] Inserting a plot into another
Look at the subplot function in the TeachingDemos package. -- 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-boun...@r- project.org] On Behalf Of Filoche Sent: Thursday, September 30, 2010 9:01 AM To: r-help@r-project.org Subject: [R] Inserting a plot into another Hi everyone. I would like to know if it was possible to insert a plot into another one. For example I have a plot and I would like to add a smaller plot in the top right corner. Best regards, Phil -- View this message in context: http://r.789695.n4.nabble.com/Inserting- a-plot-into-another-tp2720936p2720936.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-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 the output after bootstraping
Michael Larkin wrote: Thanks to the help of people from this forum I was able to bootstrap my data and then apply a model to it. Thanks for all your help. Everything worked out well, but I am having a difficult time getting the new parameter values. I bootstrapped the data 300 times and I want to get the 300 sets of parameter estimates and plot them in Excel. Here is my code: par-list(Linf=700, K=0.3, to=-0.1)#starting values for parameter estimates vb-nls(length~Linf*(1-exp(-K*(age-to))), start=par, data=small) #von Bertalanffy growth mode boo- nlsBoot(vb, niter=300) #This bootstraps my data 300 times and applies the growth model I can get R to plot the 300 parameter estimates if I do the following function: plot(boo) However, I have tried different print function options but have not had any luck getting the 300 parameter estimates. Any advice would be greatly appreciated. It's probably as simple as looking at: vb$t # but this is untested in the absence of a reproducible example. -- David. -- View this message in context: http://r.789695.n4.nabble.com/getting-the-output-after-bootstraping-tp2721024p2758846.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] please help me....
Hi, I am not entirely certain how you think we can help. Many theses/dissertations are not widely published, but there is often a copy stored at the local institution's library. I would also assume that Jan Ulbricht has a copy, which if it is not available online and it is impossible for you to visit the library, is probably your best bet. Here is his website: http://www.statistik.lmu.de/~ulbricht/index.html and I think this is the address of the statistics department at the university: Institut für Statistik Ludwig-Maximilians-Universität München Ludwigstr. 33 80539 München Best of luck to you, Josh On Thu, Sep 30, 2010 at 11:28 AM, butchman has...@hanmail.net wrote: ㅣ recently used laq packages author of this package is jan Ulbricht.. and I must read jan ulbricht's Ph.D thesis variable selection in generalized linear models but I didn't find this thesis on line and there is any method to view this thesis?? If anybody has this thesis file,please send to me... I must this guy's Ph.D. thesis... I searched LMU online library..but this thesis did not exist... please help me...very very very important issue for me.. mu e-mail address:butchma...@gmail.com -- View this message in context: http://r.789695.n4.nabble.com/please-help-me-tp2739350p2739350.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. -- Joshua Wiley Ph.D. Student, Health Psychology University of California, Los Angeles http://www.joshuawiley.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] interactive session
If you don't mind the prompt of 1:, I think
scan will do what you want:
a = scan(n=1,what='',quiet=TRUE);b = paste(t,a,sep='');
1: ada
b
[1] tada
- Phil Spector
Statistical Computing Facility
Department of Statistics
UC Berkeley
spec...@stat.berkeley.edu
On Thu, 30 Sep 2010, Peter Dalgaard wrote:
On 09/30/2010 03:33 PM, Pam wrote:
Thanks Niels but it won't do.. please copy and paste the 2 lines below
together
to your console in order to see what I mean:
cat(?); a-readLines(n=1)
b-paste(t,a,sep=)
anyone / any idea to overcome this problem?
Best,
Fatih
You might want to source() a file with those lines rather than pasting
them to the console. There's just no way you can retroactively insert
text between two already submitted lines.
You can do things like this, though
{cat(?); a-readLines(n=1)
hey
b-paste(t,a,sep=)}
However, there's a catch
{cat(?); a-readLines(n=1)
+ hey
+ b-paste(t,a,sep=)}
?ada
b
[1] tada
Notice that the hey doesn't print.
--
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Phone: (+45)38153501
Email: pd@cbs.dk Priv: pda...@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.
__
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 the output after bootstraping
Hi: To amplify on David's point, the help page of nlsBoot (package nlstools) states that the function returns the following: nlsBoot returns a list of three objects: coefboot contains the bootstrap parameter estimates bootCI contains the bootstrap medians and the bootstrap 95% confidence intervals rse is the vector of bootstrap residual errors On Thu, Sep 30, 2010 at 2:00 PM, David Winsemius dwinsem...@comcast.netwrote: Michael Larkin wrote: Thanks to the help of people from this forum I was able to bootstrap my data and then apply a model to it. Thanks for all your help. Everything worked out well, but I am having a difficult time getting the new parameter values. I bootstrapped the data 300 times and I want to get the 300 sets of parameter estimates and plot them in Excel. Based the above snippet of nlsBoot's help page, see component coefboot of your object boo below - i.e., boo$coefboot or boo[[1]]. I'll let the comment about plotting these in Excel pass... HTH, Dennis Here is my code: par-list(Linf=700, K=0.3, to=-0.1)#starting values for parameter estimates vb-nls(length~Linf*(1-exp(-K*(age-to))), start=par, data=small) #von Bertalanffy growth mode boo- nlsBoot(vb, niter=300) #This bootstraps my data 300 times and applies the growth model I can get R to plot the 300 parameter estimates if I do the following function: plot(boo) However, I have tried different print function options but have not had any luck getting the 300 parameter estimates. Any advice would be greatly appreciated. It's probably as simple as looking at: vb$t # but this is untested in the absence of a reproducible example. -- David. -- View this message in context: http://r.789695.n4.nabble.com/getting-the-output-after-bootstraping-tp2721024p2758846.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. [[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.
