[R] reliability of R-Forge?
How reliable is R-Forge? http://r-forge.r-project.org/ It is down now (for me). Reporting "R-Forge Could Not Connect to Database: " I have just started to use it for a project. It has been down for several hours (at least) on different occasions over the last couple of days. Is that common? Will it be more stable soon? Apologies if this is not an appropriate question for R-help. Dave Kane __ 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] Specify non-default libiconv, readline, libpng, tiff etc during R compilation
Hi guys, I am trying to compile 64-bit R 2.11.1 on Solaris 10. I mainly follow the guide in here https://www.initworks.com/wiki/pages/viewpage.action?pageId=6521038 The guide suggests that install the customised libiconv, readline under the designated R installation folder and become the private libraries that are exclusive to that R only. This seems to be a great solution, as I need to compile several versions of R in the same environment, libiconv, readline, libpng etc would be different version too. This concept will allow me to do so without contaminating the global default libraries. But I found that, R ./configure does not locate the customised, provate libiconv, readline etc in the designated location. It still looks from the default one like /usr/lib, /usr/local/lib etc (typically, is older version). I am wondering if any R user has done the similar work and can share some experience. Regards, Derrick [[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] Looking for an image (R 64-bit on Linux 64-bit) on Amazon EC2
>> You have a 64 bit Linux? If so... >>Dowload the sources Do you mean download Linux kernel source code and then compile it on Amazon EC2? -- View this message in context: http://r.789695.n4.nabble.com/Looking-for-an-image-R-64-bit-on-Linux-64-bit-on-Amazon-EC2-tp2338938p2339072.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] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
Hi Ravi,
Using days and dividing it by 365 effectively converts the number to years
anyway and allows for the irregular times to be specific to the days.
Also, when I replace dates[1] in your line:
times <- as.numeric(difftime(dates, dates[1], units="days") /
365.24) with "2010-08-24" I think I am getting some irregular results.
Effectively, what I was trying to do was match what Excel produced with its
XIRR function. With the example I gave excel returned an IRR of ~0.37 (or 37%)
I am still in the process of debugging it...
-Original Message-
From: Ravi Varadhan [mailto:rvarad...@jhmi.edu]
Sent: Wednesday, August 25, 2010 7:24 PM
To: Adrian Ng
Cc: r-help@r-project.org
Subject: RE: [R] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
The secant method converges just fine. Your problem might have occurred due
to improper conversion of dates to elapsed time. You want to calculate IRR
using "year" as the time unit, not "days".
Here is the secant function (modified to account for irregular times) and
the results for your example:
ANXIRR <- function (cashFlow, dates, guess, tol=1.e-04){
npv <- function (cashFlow, times, irr) {
n <- length(cashFlow)
sum(cashFlow / (1 + irr)^times)
}
if (guess == 0)stop("Initial guess must be strictly greater than 0")
times <- as.numeric(difftime(dates, dates[1], units="days") /
365.24)
irrprev <- c(0)
irr <- guess
pvPrev <- sum(cashFlow)
pv <- npv(cashFlow, times, irr)
eps <- abs(pv-pvPrev)
while (eps >= tol) {
tmp <- irrprev
irrprev <- irr
irr <- irr - ((irr - tmp) * pv / (pv - pvPrev))
pvPrev <- pv
pv <- npv(cashFlow, times, irr)
eps <- abs(pv - pvPrev)
}
list(irr = irr, npv = pv)
}
CF <- c(-1000,500,500,500,500,500)
dates <-
c("1/1/2001","2/1/2002","3/1/2003","4/1/2004","5/1/2005","6/1/2006")
ANXIRR(CF, dates, guess=0.1)
> ANXIRR(CF, dates, guess=0.1)
$irr
[1] 0.4106115
$npv
[1] 2.984279e-13
Ravi.
-Original Message-
From: Adrian Ng [mailto:a...@hamiltonlane.com]
Sent: Wednesday, August 25, 2010 6:23 PM
To: Ravi Varadhan
Subject: RE: [R] Secant Method Convergence (Method to replicate Excel
XIRR/IRR)
The forum is kind of slow so I'm just re-sending you the message here:
Hi Ravi,
I'm just trying a fairly simple example:
CFs: -1000,500,500,500,500,500
dates<-c("1/1/2001","2/1/2002","3/1/2003","4/1/2004","5/1/2005","6/1/2006")
Thanks a lot for your help.
Adrian
-Original Message-
From: Ravi Varadhan [mailto:rvarad...@jhmi.edu]
Sent: Wednesday, August 25, 2010 5:44 PM
To: Adrian Ng; r-help@r-project.org
Subject: RE: [R] Secant Method Convergence (Method to replicate Excel
XIRR/IRR)
Yes, the secant method (like Newton Raphson) is not guaranteed to converge,
unlike the bisection method, but it has a superlinear convergence (not that
this matters much!). Brent's method, which is used in `uniroot', is a
reliable and fast method, which is why I suggested it in my previous email.
Having said that, I am not sure about the convergence problem that you are
having without seeing the actual example.
Ravi.
-Original Message-
From: Adrian Ng [mailto:a...@hamiltonlane.com]
Sent: Wednesday, August 25, 2010 5:28 PM
To: Ravi Varadhan; r-help@r-project.org
Subject: RE: [R] Secant Method Convergence (Method to replicate Excel
XIRR/IRR)
Hi Ravi,
Thanks for the responses. I was actually trying to calculate IRR based on
unevenly spaced cash flows, and that's why I decided to use the secant
method. I'm not sure if my answer isn't converging because I have some
careless mistake in the code, or if it's simply because unlike the bisection
method, the secant method doesn't 'sandwich' the desired root.
-Original Message-
From: Ravi Varadhan [mailto:rvarad...@jhmi.edu]
Sent: Wednesday, August 25, 2010 5:24 PM
To: Adrian Ng; r-help@r-project.org
Subject: RE: [R] Secant Method Convergence (Method to replicate Excel
XIRR/IRR)
Another approach is to use `uniroot' to find the zero of the NPV function:
npv <- function (cashFlow, irr) {
n <- length(cashFlow)
sum(cashFlow / (1 + irr)^{0: (n-1)})
}
uniroot(f=npv, interval=c(0,1), cashFlow=cashFlow)
However, there may be situations where there are no real zeros or there are
multiple zeros of the NPV function.
Ravi.
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Adrian Ng
Sent: Wednesday, August 25, 2010 8:39 AM
To: r-help@r-project.org
Subject: [R] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
Hi,
I am new to R, and as a first exercise, I decided to try to implement an
XIRR function using the secant method. I did a quick search and saw another
posting that used the Bisection method but wanted to see if it was possible
using the secant method.
I would input a Cash Flow and Date vector as well as an i
Re: [R] creation package
Dear r-help,
I took your advice into consideration and i tried to solve my errors. But ,
when I use the command R CMD check namepackage, I find an error in this file
C:/Rpack/namepackage.Rcheck/00check.txt :
* using log directory 'C:/Rpack/namepackge.Rcheck'
* using R version 2.10.1 (2009-12-14)
* using session charset: ISO8859-1
* checking for file 'namepackge/DESCRIPTION' ... OK
* checking extension type ... Package
* this is package 'namepackge' version '1.0'
* checking package name space information ... OK
* checking package dependencies ... OK
* checking if this is a source package ... OK
* checking for executable files ... OK
* checking whether package 'namepackge' can be installed ... OK
* checking package directory ... OK
* checking for portable file names ... OK
* checking DESCRIPTION meta-information ... OK
* checking top-level files ... OK
* checking index information ... OK
* checking package subdirectories ... OK
* checking R files for non-ASCII characters ... OK
* checking R files for syntax errors ... OK
* checking whether the package can be loaded ... OK
* checking whether the package can be loaded with stated dependencies ... OK
* checking whether the name space can be loaded with stated dependencies ...
OK
* checking for unstated dependencies in R code ... OK
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
* checking R code for possible problems ... OK
* checking Rd files ... OK
* checking Rd metadata ... OK
* checking Rd cross-references ... OK
* checking for missing documentation entries ... OK
* checking for code/documentation mismatches ... OK
* checking Rd \usage sections ... OK
* checking data for non-ASCII characters ... OK
* checking examples ... ERROR
Running examples in 'namepackge-Ex.R' failed.
The error most likely occurred in:
> ### * RT
>
> flush(stderr()); flush(stdout())
>
> ### Name: RT
> ### Title: modeling
> ### Aliases: RT
> ### Keywords: programming models
>
> ### ** Examples
>
> data(d)
> data <- d
> RT(t(d), c(2,2,2), c(1,2,3),c("V1","V2","V3"))
Error in text.default(x, y, labels = labels, col = label.color, family =
label.family, :
family 'serif' not included in PostScript device
Calls: RT-> plot -> plot.igraph -> text -> text.default
Execution halted
Best Regards,
2010/8/19, Michael Dewey :
>
> At 13:57 19/08/2010, anderson nuel wrote:
>
>> Dear r-help,
>>
>> I don't use namespace.
>>
>
> Well, as I said in my original reply, it would be a good idea to do so.
>
>
>> How can I make asia available?
>>
>
> Without knowing where asia is that is quite a tough call. How do you access
> it when you test your code before you try to package it?
>
>
>> I think my problem in creating the package in this: I have a singleglobal
>> function (RT) in my package, but inside RT I need to call several other
>> function( comb l,earn_comb, nchoo). When I used package.sekeleton, Iput in
>> lists all the functions(comb l,earn_comb, nchoo,RT),but in the 'man' I
>> left
>> only namepackage-package.Rd and RT.Rd. When I did do this, is it true??
>>
>>
>> Best Regards
>>
>>
>> 2010/8/18, Michael Dewey <
>> i...@aghmed.fsnet.co.uk>:
>>
>> At 10:27 18/08/2010, anderson nuel wrote:
>> Dear r-help,
>>
>> No, I find errors in the file C:/Rp/namepackage.Rcheck/00check.txt :
>>
>> * using log directory 'C:/Rpackage/namepackage.Rcheck'
>> * using R version 2.10.1 (2009-12-14)
>> * using session charset: ISO8859-1
>> * checking for file 'namepackage/DESCRIPTION' ... OK
>> * checking extension type ... Package
>> * this is package 'namepackage' version '1.0'
>> * checking package dependencies ... OK
>> * checking if this is a source package ... OK
>> * checking for executable files ... OK
>> * checking whether package 'namepackage' can be installed ... OK
>> * checking package directory ... OK
>> * checking for portable file names ... OK
>> * checking DESCRIPTION meta-information ... OK
>> * checking top-level files ... OK
>> * checking index information ... OK
>> * checking package subdirectories ... OK
>> * checking R files for non-ASCII characters ... OK
>> * checking R files for syntax errors ... OK
>> * checking whether the package can be loaded ... OK
>> * checking whether the package can be loaded with stated dependencies ...
>> OK
>> * checking for unstated dependencies in R code ... OK
>> * checking S3 generic/method consistency ... OK
>> * checking replacement functions ... OK
>> * checking foreign function calls ... OK
>> * checking R code for possible problems ... NOTE
>> comb: no visible global function definition for 'copy_bloc'
>> comb: no visible global function definition for 'copy_interl'
>> * checking Rd files ... OK
>> * checking Rd metadata ... OK
>> * checking Rd cross-references ... OK
>> * checking for missing documentation entries ... WARNING
>> Undocumented code objects:
>> comb learn_comb nchoo
>> All user-level objects in a package should have d
Re: [R] list of closures
Unless for learning R, you should really consider R.oo or proto packages
that may be more convenient for you (but you don't provide enough
context to tell).
Best,
Philippe
On 26/08/10 06:28, Gabor Grothendieck wrote:
On Thu, Aug 26, 2010 at 12:04 AM, Stephen T. wrote:
Hi, I wanted to create a list of closures. When I use Map(), mapply(),
lapply(), etc., to create this list, it appears that the wrong arguments are
being passed to the main function. For example:
Main function:
adder<- function(x) function(y) x + y
Creating list of closures with Map():
plus<- Map(adder,c(one=1,two=2))> plus$one(1)[1] 3> plus$two(1)[1] 3
Examining what value was bound to "x":
Map(function(fn) get("x",environment(fn)),plus)$one[1] 2$two[1] 2
This is what I had expected:
plus<- list(one=adder(1),two=adder(2))> plus$one(1)[1] 2> plus$two(1)[1] 3
Map(function(fn) get("x",environment(fn)),plus)$one[1] 1$two[1] 2
Anyone know what's going on? Thanks much!
R uses lazy evaluation of function arguments. Try forcing x:
adder<- function(x) { force(x); function(y) x + y }
plus<- Map(adder,c(one=1,two=2))
plus$one(1)
[1] 2
plus$two(1)
[1] 3
__
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] list of closures
On Thu, Aug 26, 2010 at 12:04 AM, Stephen T. wrote:
>
> Hi, I wanted to create a list of closures. When I use Map(), mapply(),
> lapply(), etc., to create this list, it appears that the wrong arguments are
> being passed to the main function. For example:
> Main function:
>> adder <- function(x) function(y) x + y
> Creating list of closures with Map():
>> plus <- Map(adder,c(one=1,two=2))> plus$one(1)[1] 3> plus$two(1)[1] 3
> Examining what value was bound to "x":
>> Map(function(fn) get("x",environment(fn)),plus)$one[1] 2$two[1] 2
>
> This is what I had expected:
>> plus <- list(one=adder(1),two=adder(2))> plus$one(1)[1] 2> plus$two(1)[1] 3
>> Map(function(fn) get("x",environment(fn)),plus)$one[1] 1$two[1] 2
>
> Anyone know what's going on? Thanks much!
R uses lazy evaluation of function arguments. Try forcing x:
> adder <- function(x) { force(x); function(y) x + y }
> plus <- Map(adder,c(one=1,two=2))
> plus$one(1)
[1] 2
> plus$two(1)
[1] 3
__
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] list of closures
Hi, I wanted to create a list of closures. When I use Map(), mapply(),
lapply(), etc., to create this list, it appears that the wrong arguments are
being passed to the main function. For example:
Main function:
> adder <- function(x) function(y) x + y
Creating list of closures with Map():
> plus <- Map(adder,c(one=1,two=2))> plus$one(1)[1] 3> plus$two(1)[1] 3
Examining what value was bound to "x":
> Map(function(fn) get("x",environment(fn)),plus)$one[1] 2$two[1] 2
This is what I had expected:
> plus <- list(one=adder(1),two=adder(2))> plus$one(1)[1] 2> plus$two(1)[1] 3
> Map(function(fn) get("x",environment(fn)),plus)$one[1] 1$two[1] 2
Anyone know what's going on? Thanks much!
Stephen
[[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] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
I have written a general root-finder using the secant method.
secant <- function(par, fn, tol=1.e-07, itmax = 100, trace=TRUE, ...) {
# par = a starting vector with 2 starting values
# fn = a function whose first argument is the variable of interest
#
if (length(par) != 2) stop("You must specify a starting parameter vector of
length 2")
p.2 <- par[1]
p.1 <- par[2]
f <- rep(NA, length(par))
f[1] <- fn(p.1, ...)
f[2] <- fn(p.2, ...)
iter <- 1
pchg <- abs(p.2 - p.1)
fval <- f[2]
if (trace) cat("par: ", par, "fval: ", f, "\n")
while (pchg >= tol & abs(fval) > tol & iter <= itmax) {
p.new <- p.2 - (p.2 - p.1) * f[2] / (f[2] - f[1])
pchg <- abs(p.new - p.2)
fval <- fn(p.new, ...)
p.1 <- p.2
p.2 <- p.new
f[1] <- f[2]
f[2] <- fval
iter <- iter + 1
if (trace) cat("par: ", p.new, "fval: ", fval, "\n")
}
list(par = p.new, value = fval, iter=iter)
}
Now we can use this function to find the zero of your NPV function as follows:
npv <- function (irr, cashFlow, times) sum(cashFlow / (1 + irr)^times)
CF <- c(-1000,500,500,500,500,500)
dates <- c("1/1/2001","2/1/2002","3/1/2003","4/1/2004","5/1/2005","6/1/2006")
cfDate <- as.Date(cfDate,format="%m/%d/%Y")
times <- as.numeric(difftime(cfDate, cfDate[1], units="days"))/365.24
secant(par=c(0,0.1), fn=npv, cashFlow=CF, times=times)
> secant(par=c(0,0.1), fn=npv, cashFlow=CF, times=times)
par: 0 0.1 fval: 854.2388 1500
par: 0.232284 fval: 334.7318
par: 0.2990093 fval: 156.9595
par: 0.3579227 fval: 30.59229
par: 0.3721850 fval: 3.483669
par: 0.3740179 fval: 0.08815743
par: 0.3740655 fval: 0.0002613245
par: 0.3740656 fval: 1.966778e-08
$par
[1] 0.3740656
$value
[1] 1.966778e-08
$iter
[1] 8
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
Date: Wednesday, August 25, 2010 10:51 pm
Subject: Re: [R] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
To: Adrian Ng
Cc: "r-help@r-project.org"
> You should use cfDate[1] as the time origin. You cannot use
> 08-24-2010 as the time origin, since that will yield negative times.
>
>
> Here is the correct solution.
>
> ANXIRR <- function (cashFlow, dates, guess, tol=1.e-04){
>
> npv <- function (cashFlow, times, irr) {
> n <- length(cashFlow)
> sum(cashFlow / (1 + irr)^times)
> }
>
> if (guess == 0) stop("Initial guess must be strictly greater than 0")
>
>
> cfDate <- as.Date(cfDate,format="%m/%d/%Y")
> times <- as.numeric(difftime(cfDate, cfDate[1], units="days") /
> 365.24)
>
> irrprev <- c(0)
> irr <- guess
> pvPrev <- sum(cashFlow)
> pv <- npv(cashFlow, times, irr)
> eps <- abs(pv-pvPrev)
>
> while (eps >= tol) {
> tmp <- irrprev
> irrprev <- irr
> irr <- irr - ((irr - tmp) * pv / (pv - pvPrev))
> pvPrev <- pv
> pv <- npv(cashFlow, times, irr)
> eps <- abs(pv - pvPrev)
> }
> list(irr = irr, npv = pv)
> }
> CF <- c(-1000,500,500,500,500,500)
>
> dates <-
> c("1/1/2001","2/1/2002","3/1/2003","4/1/2004","5/1/2005","6/1/2006")
>
> ANXIRR(CF, dates, guess=0.1)
>
> > ANXIRR(CF, dates, guess=0.1)
> $irr
> [1] 0.3740656
>
> $npv
> [1] 2.102695e-09
>
>
> 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: Adrian Ng
> Date: Wednesday, August 25, 2010 8:33 pm
> Subject: RE: [R] Secant Method Convergence (Method to replicate Excel
> XIRR/IRR)
> To: Ravi Varadhan
> Cc: "r-help@r-project.org"
>
>
> > Hi Ravi,
> >
> > Using days and dividing it by 365 effectively converts the number
> to
> > years anyway and allows for the irregular times to be specific to
> the
> > days.
> >
> > Also, when I replace dates[1] in your line:
> > times <- as.numeric(difftime(dates, dates[1], units="days") /
> > 365.24) with "2010-08-24" I think I am getting some irregular
> > results.
> >
> > Effectively, what I was trying to do was match what Excel produced
>
> > with its XIRR function. With the example I gave excel returned an
> IRR
> > of ~0.37 (or 37%)
> >
> > I am still in the process of debugging it...
> >
> >
> >
> >
> >
> > -Original Message-
> > From: Ravi Varadhan [
> > Sent: Wednesday, August 25, 2010 7:24 PM
> > To: Adrian Ng
> > Cc: r-help@r-project.org
> > Subject: RE: [R] Secant Method Convergence (Method to replicate
> Excel
> > XIRR/IRR)
> >
> > The secant method converges just fine. Your problem m
Re: [R] How to obtain the graph of fitted values against one variable after estimation?
On Aug 25, 2010, at 10:46 PM, Le Wang wrote: Hi there, I have a question regarding the "plot" command after estimation. Specifically, I estimate a model, say regressing y on x and z. And after estimation, I would like to plot the fitted values against x, but I don't need that for z. The following command always gives two graphs, for both variables x and z. plot.np<-plot(model, plot.errors.method = "asymptotic") My question is, what option should I specify in order to get the graph for x only? Pick a constant value for "z" and vary "x" in a dataframe that you offer to the newdata argument of predict. ?predict Then plot those values versus x. I know this is probably a very simple question, but I searched around for a while without any luck. Thank you for your time. -- David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
You should use cfDate[1] as the time origin. You cannot use 08-24-2010 as the
time origin, since that will yield negative times.
Here is the correct solution.
ANXIRR <- function (cashFlow, dates, guess, tol=1.e-04){
npv <- function (cashFlow, times, irr) {
n <- length(cashFlow)
sum(cashFlow / (1 + irr)^times)
}
if (guess == 0) stop("Initial guess must be strictly greater than 0")
cfDate <- as.Date(cfDate,format="%m/%d/%Y")
times <- as.numeric(difftime(cfDate, cfDate[1], units="days") /
365.24)
irrprev <- c(0)
irr <- guess
pvPrev <- sum(cashFlow)
pv <- npv(cashFlow, times, irr)
eps <- abs(pv-pvPrev)
while (eps >= tol) {
tmp <- irrprev
irrprev <- irr
irr <- irr - ((irr - tmp) * pv / (pv - pvPrev))
pvPrev <- pv
pv <- npv(cashFlow, times, irr)
eps <- abs(pv - pvPrev)
}
list(irr = irr, npv = pv)
}
CF <- c(-1000,500,500,500,500,500)
dates <- c("1/1/2001","2/1/2002","3/1/2003","4/1/2004","5/1/2005","6/1/2006")
ANXIRR(CF, dates, guess=0.1)
> ANXIRR(CF, dates, guess=0.1)
$irr
[1] 0.3740656
$npv
[1] 2.102695e-09
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: Adrian Ng
Date: Wednesday, August 25, 2010 8:33 pm
Subject: RE: [R] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
To: Ravi Varadhan
Cc: "r-help@r-project.org"
> Hi Ravi,
>
> Using days and dividing it by 365 effectively converts the number to
> years anyway and allows for the irregular times to be specific to the
> days.
>
> Also, when I replace dates[1] in your line:
> times <- as.numeric(difftime(dates, dates[1], units="days") /
> 365.24) with "2010-08-24" I think I am getting some irregular
> results.
>
> Effectively, what I was trying to do was match what Excel produced
> with its XIRR function. With the example I gave excel returned an IRR
> of ~0.37 (or 37%)
>
> I am still in the process of debugging it...
>
>
>
>
>
> -Original Message-
> From: Ravi Varadhan [
> Sent: Wednesday, August 25, 2010 7:24 PM
> To: Adrian Ng
> Cc: r-help@r-project.org
> Subject: RE: [R] Secant Method Convergence (Method to replicate Excel
> XIRR/IRR)
>
> The secant method converges just fine. Your problem might have
> occurred due
> to improper conversion of dates to elapsed time. You want to
> calculate IRR
> using "year" as the time unit, not "days".
>
> Here is the secant function (modified to account for irregular times)
> and
> the results for your example:
>
> ANXIRR <- function (cashFlow, dates, guess, tol=1.e-04){
>
> npv <- function (cashFlow, times, irr) {
> n <- length(cashFlow)
> sum(cashFlow / (1 + irr)^times)
> }
>
> if (guess == 0)stop("Initial guess must be strictly greater than 0")
>
>
> times <- as.numeric(difftime(dates, dates[1], units="days") /
> 365.24)
>
> irrprev <- c(0)
>irr <- guess
> pvPrev <- sum(cashFlow)
> pv <- npv(cashFlow, times, irr)
> eps <- abs(pv-pvPrev)
>
> while (eps >= tol) {
> tmp <- irrprev
> irrprev <- irr
> irr <- irr - ((irr - tmp) * pv / (pv - pvPrev))
> pvPrev <- pv
> pv <- npv(cashFlow, times, irr)
> eps <- abs(pv - pvPrev)
> }
> list(irr = irr, npv = pv)
> }
>
> CF <- c(-1000,500,500,500,500,500)
>
> dates <-
> c("1/1/2001","2/1/2002","3/1/2003","4/1/2004","5/1/2005","6/1/2006")
>
>
> ANXIRR(CF, dates, guess=0.1)
>
> > ANXIRR(CF, dates, guess=0.1)
> $irr
> [1] 0.4106115
>
> $npv
> [1] 2.984279e-13
>
>
> Ravi.
>
> -Original Message-
> From: Adrian Ng [
> Sent: Wednesday, August 25, 2010 6:23 PM
> To: Ravi Varadhan
> Subject: RE: [R] Secant Method Convergence (Method to replicate Excel
> XIRR/IRR)
>
> The forum is kind of slow so I'm just re-sending you the message here:
>
> Hi Ravi,
>
> I'm just trying a fairly simple example:
> CFs: -1000,500,500,500,500,500
>
> dates<-c("1/1/2001","2/1/2002","3/1/2003","4/1/2004","5/1/2005","6/1/2006")
>
>
> Thanks a lot for your help.
> Adrian
>
> -Original Message-
> From: Ravi Varadhan [
> Sent: Wednesday, August 25, 2010 5:44 PM
> To: Adrian Ng; r-help@r-project.org
> Subject: RE: [R] Secant Method Convergence (Method to replicate Excel
> XIRR/IRR)
>
> Yes, the secant method (like Newton Raphson) is not guaranteed to converge,
> unlike the bisection method, but it has a superlinear convergence
> (not that
> this matters much!). Brent's method, which is used in `uniroot', is
> a
> reliable and fast method, which is why I suggested it in my previous
> email.
>
> Having said that, I am not sure about the convergence problem that
> you are
> having without seeing the actual exa
[R] How to obtain the graph of fitted values against one variable after estimation?
Hi there, I have a question regarding the "plot" command after estimation. Specifically, I estimate a model, say regressing y on x and z. And after estimation, I would like to plot the fitted values against x, but I don't need that for z. The following command always gives two graphs, for both variables x and z. plot.np<-plot(model, plot.errors.method = "asymptotic") My question is, what option should I specify in order to get the graph for x only? I know this is probably a very simple question, but I searched around for a while without any luck. Thank you for your time. Best, Le __ 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] Rotate x-axis label on log scale
Hi Jim
Thanks so much for this. The updated staxlab function works perfectly!
Tim
-Original Message-
From: Jim Lemon [mailto:j...@bitwrit.com.au]
Sent: Wednesday, August 25, 2010 6:36 PM
To: tesutton
Cc: r-help@r-project.org
Subject: Re: [R] Rotate x-axis label on log scale
On 08/25/2010 01:37 PM, Tim Elwell-Sutton wrote:
> Hi Jim
> Thanks for this. The staxlab function seems very useful. Unfortunately,
the
> rotation option doesn't seem to work for me when the y-axis is on a log
> scale. It will stagger the labels but not rotate them. There's no error
> message. On a linear axis the rotation works nicely. Any ideas?
> The example below works if you omit log='y' or srt=45
>
Hi Tim,
Thanks for letting me know about this problem. I never did get around to
making staxlab work with log axes, but I think this new version:
staxlab<-function(side=1,at,labels,nlines=2,top.line=0.5,
line.spacing=0.8,srt=NA,ticklen=0.03,...) {
if(missing(labels)) labels<-at
nlabels<-length(labels)
if(missing(at)) at<-1:nlabels
if(is.na(srt)) {
linepos<-rep(top.line,nlines)
for(i in 2:nlines) linepos[i]<-linepos[i-1]+line.spacing
linepos<-rep(linepos,ceiling(nlabels/nlines))[1:nlabels]
axis(side=side,at=at,labels=rep("",nlabels))
mtext(text=labels,side=side,line=linepos,at=at,...)
}
else {
xylim<-par("usr")
if(side == 1) {
xpos<-at
if(par("ylog")) ypos<-10^(xylim[3]-ticklen*(xylim[4]-xylim[3]))
else ypos<-xylim[3]-ticklen*(xylim[4]-xylim[3])
}
else {
ypos<-at
if(par("xlog")) xpos<-10^(xylim[1]-ticklen*(xylim[2]-xylim[1]))
else xpos<-xylim[1]-ticklen*(xylim[2]-xylim[1])
}
par(xpd=TRUE)
text(xpos,ypos,labels,srt=srt,adj=1,...)
par(xpd=FALSE)
}
}
will do what you want.
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] Adding column to dataframe
To update on this. I ran the same command on a grid of computers with 32gb ram, and it completed in 15 seconds, compared to the ~20 minutes on my desktop. Simply a ram issue as suspected by a few on the list here. Thanks -- View this message in context: http://r.789695.n4.nabble.com/Adding-column-to-dataframe-tp2330556p2339076.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] lattice help required
Deepayan Sarkar has a function combineLimits() in the development
version of the latticeExtra package (i.e. the version on
r-forge.r-project.org) which will set common scales in each row or
column of your layout. It can also remove the internal axes.
# Felix
On 26 August 2010 04:43, Kay Cichini wrote:
>
> .. thanks again, richard.
> and you swiftly saw the next problem comming up - when using par.settings =
> list(layout.widths = list(axis.panel = c(1, 0))) getting rid of the double
> tick labeling would be natural -
> but i'll leave it at that for today.
>
> many thanks,
> kay
>
>
> Richard M. Heiberger wrote:
>>
>> The multiple y axes are protecting you in this situation.
>>
>>
>> z <- cbind(rnorm(100,c(1,10),1), rnorm(100,c(20,30),1))
>> dotplot(z[,1]+z[,2] ~ facs$Treatment|facs$Sites,
>> outer=TRUE,
>> scales = list(
>> y = list(
>> relation="free")),
>> ylab=c("y1", "y2"),
>> xlab=c("Site 1", "Site 2"),
>> strip=FALSE,
>> main="problem")
>>
>> dotplot(z[,1]+z[,2] ~ facs$Treatment|facs$Sites,
>> outer=TRUE,
>> scales = list(
>> y = list(
>> relation="free",
>> limits=list(c(-5,13),c(-5,13),c(18,32),c(18,32,
>> ylab=c("y1", "y2"),
>> xlab=c("Site 1", "Site 2"),
>> strip=FALSE, main="protected")
>>
>> For more control (such as suppressing the y-tick labels in the right-hand
>> column,
>> I recommend Deepayan Sarkar's book.
>>
>> Rich
>>
>> [[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.
>>
>>
>
>
> -
>
> Kay Cichini
> Postgraduate student
> Institute of Botany
> Univ. of Innsbruck
>
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/lattice-help-required-tp2338382p2338707.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.
>
--
Felix Andrews / 安福立
http://www.neurofractal.org/felix/
__
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] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
The secant method converges just fine. Your problem might have occurred due
to improper conversion of dates to elapsed time. You want to calculate IRR
using "year" as the time unit, not "days".
Here is the secant function (modified to account for irregular times) and
the results for your example:
ANXIRR <- function (cashFlow, dates, guess, tol=1.e-04){
npv <- function (cashFlow, times, irr) {
n <- length(cashFlow)
sum(cashFlow / (1 + irr)^times)
}
if (guess == 0)stop("Initial guess must be strictly greater than 0")
times <- as.numeric(difftime(dates, dates[1], units="days") /
365.24)
irrprev <- c(0)
irr <- guess
pvPrev <- sum(cashFlow)
pv <- npv(cashFlow, times, irr)
eps <- abs(pv-pvPrev)
while (eps >= tol) {
tmp <- irrprev
irrprev <- irr
irr <- irr - ((irr - tmp) * pv / (pv - pvPrev))
pvPrev <- pv
pv <- npv(cashFlow, times, irr)
eps <- abs(pv - pvPrev)
}
list(irr = irr, npv = pv)
}
CF <- c(-1000,500,500,500,500,500)
dates <-
c("1/1/2001","2/1/2002","3/1/2003","4/1/2004","5/1/2005","6/1/2006")
ANXIRR(CF, dates, guess=0.1)
> ANXIRR(CF, dates, guess=0.1)
$irr
[1] 0.4106115
$npv
[1] 2.984279e-13
Ravi.
-Original Message-
From: Adrian Ng [mailto:a...@hamiltonlane.com]
Sent: Wednesday, August 25, 2010 6:23 PM
To: Ravi Varadhan
Subject: RE: [R] Secant Method Convergence (Method to replicate Excel
XIRR/IRR)
The forum is kind of slow so I'm just re-sending you the message here:
Hi Ravi,
I'm just trying a fairly simple example:
CFs: -1000,500,500,500,500,500
dates<-c("1/1/2001","2/1/2002","3/1/2003","4/1/2004","5/1/2005","6/1/2006")
Thanks a lot for your help.
Adrian
-Original Message-
From: Ravi Varadhan [mailto:rvarad...@jhmi.edu]
Sent: Wednesday, August 25, 2010 5:44 PM
To: Adrian Ng; r-help@r-project.org
Subject: RE: [R] Secant Method Convergence (Method to replicate Excel
XIRR/IRR)
Yes, the secant method (like Newton Raphson) is not guaranteed to converge,
unlike the bisection method, but it has a superlinear convergence (not that
this matters much!). Brent's method, which is used in `uniroot', is a
reliable and fast method, which is why I suggested it in my previous email.
Having said that, I am not sure about the convergence problem that you are
having without seeing the actual example.
Ravi.
-Original Message-
From: Adrian Ng [mailto:a...@hamiltonlane.com]
Sent: Wednesday, August 25, 2010 5:28 PM
To: Ravi Varadhan; r-help@r-project.org
Subject: RE: [R] Secant Method Convergence (Method to replicate Excel
XIRR/IRR)
Hi Ravi,
Thanks for the responses. I was actually trying to calculate IRR based on
unevenly spaced cash flows, and that's why I decided to use the secant
method. I'm not sure if my answer isn't converging because I have some
careless mistake in the code, or if it's simply because unlike the bisection
method, the secant method doesn't 'sandwich' the desired root.
-Original Message-
From: Ravi Varadhan [mailto:rvarad...@jhmi.edu]
Sent: Wednesday, August 25, 2010 5:24 PM
To: Adrian Ng; r-help@r-project.org
Subject: RE: [R] Secant Method Convergence (Method to replicate Excel
XIRR/IRR)
Another approach is to use `uniroot' to find the zero of the NPV function:
npv <- function (cashFlow, irr) {
n <- length(cashFlow)
sum(cashFlow / (1 + irr)^{0: (n-1)})
}
uniroot(f=npv, interval=c(0,1), cashFlow=cashFlow)
However, there may be situations where there are no real zeros or there are
multiple zeros of the NPV function.
Ravi.
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Adrian Ng
Sent: Wednesday, August 25, 2010 8:39 AM
To: r-help@r-project.org
Subject: [R] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
Hi,
I am new to R, and as a first exercise, I decided to try to implement an
XIRR function using the secant method. I did a quick search and saw another
posting that used the Bisection method but wanted to see if it was possible
using the secant method.
I would input a Cash Flow and Date vector as well as an initial guess. I
hardcoded today's initial date so I could do checks in Excel. This code
seems to only converge when my initial guess is very close to the correct
IRR.
Maybe I have some basic errors in my coding/logic? Any help would be greatly
appreciated.
The Wikipedia article to secant method and IRR:
http://en.wikipedia.org/wiki/Internal_rate_of_return#Numerical_solution
Thanks!
ANXIRR <- function (cashFlow, cfDate, guess){
cfDate<-as.Date(cfDate,format="%m/%d/%Y")
irrprev <- c(0); irr<- guess
pvPrev<- sum(cashFlow)
pv<-
sum(cashFlow/((1+irr)^(as.numeric(difftime(cfDate,"2010-08-24",units="days")
)/360)))
print(pv)
print("Hi")
while (abs(pv) >= 0.001) {
t<-irrprev; irrprev<- irr;
irr<-irr-((irr-t
Re: [R] how to plot y-axis on the right of x-axis
Yes, I appreciated your answers which well hit my questions. (esp the perfect model parts). About the plot part, one more question. Is it possible to make the two plots (northern and southern richness) sharing the same Y axis (latitude = 0, the equator) ? In other words, southern richness would be on the left side of the Y axis, while northern one on the right side. Elaine Two plot panels on the same device like: layout(1:2) plot(1:100) plot(1:100) layout(1) => Yes that's what I want The second panel for southern species richness, do you mean you want the plot to go like this: plot(1:100, xlim = c(100,1)) => Yes, too. [[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 plot y-axis on the right of x-axis
Thank you.
The answers provides the right direction and a code was written accordingly.
One more request, if the label of axis X wants to be drawn from 5 to 1 (left
to right)
rather than 1 to 5, is it fine to change axis (4, at = NULL) ?
If so, which value should be input ?
Thanks again.
Elaine
code
plot(1:5,axes=FALSE)
axis(1)
axis(4)
box()
On Wed, Aug 25, 2010 at 5:39 PM, Jim Lemon wrote:
> On 08/25/2010 09:12 AM, elaine kuo wrote:
>
>> Dear List,
>>
>> I have a richness data distributing across 20 N to 20 S latitude. (120
>> E-140
>> E longitude).
>>
>>
>> I would like to draw the richness in the north hemisphere and a regression
>> line in the plot
>> (x-axis: latitude, y-axis: richness in the north hemisphere).
>> The above demand is done using plot.
>>
>> Then, south hemisphere richness and regression are required to be
>> generated
>> using
>> the same y-axis above but an x-axis on the left side of the y-axis.
>> (The higher latitude in the south hemisphere, the left it would move)
>>
>> Please kindly share how to design the south plot and regression line for
>> richness.
>> Also, please advise if any more info is in need.
>>
>> Hi Elaine,
> Changing the position of the axes is easily done by changing the "side" and
> "pos" arguments of the "axis" function. If you want to move the y-axis to
> the right (or as you state it, the x axis to the left):
>
> # y axis on the left
> plot(1:5,axes=FALSE)
> axis(1)
> axis(2)
> # add a y axis one third of the way to the right
> xylim<-par("usr")
> axis(2,pos=xylim[1]+diff(xylim[1:2])/3)
> # add another y axis two thirds of the way
> axis(4,pos=xylim[2]-diff(xylim[1:2])/3)
> # add one more on the right
> axis(4)
>
> You can move the x axis up and down using the same tricks.
>
> Jim
>
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] OT: R for iPhone/iPad OS?
Hey, Many thanks to all who responded! I'll pass along the links to my pals. Personally I can't imagine doing R on a virtual keyboard in the first place, but to each his own. Carl __ 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] Looking for an image (R 64-bit on Linux 64-bit) on Amazon EC2
You have a 64 bit Linux? If so... Dowload the sources make sure you have all of the dependencies unpack tarball cd to de-compressed directory issue ./configure make sudo make install or maybe use your distros packages managment tool. Am I missing something? On Wed, Aug 25, 2010 at 4:39 PM, noclue_ wrote: > > I have found an existing image on Amazon EC2 including R. But unfortunately, > it is 32-bit > R on 32-bit Linux. > > Does anybody know if there exists an mage (R 64-bit on Linux 64-bit) on > Amazon EC2? > > Or how can I install 64-bit R on my own Linux instance there? > > Thanks. > -- > View this message in context: > http://r.789695.n4.nabble.com/Looking-for-an-image-R-64-bit-on-Linux-64-bit-on-Amazon-EC2-tp2338938p2338938.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. > -- 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] OT: R for iPhone/iPad OS?
On 8/25/2010 2:23 PM, Carl Witthoft wrote: No, seriously: I've had more than one person at work wonder what math toolset could be loaded onto iOS. So, before Matlab, FreeMat, Mathematia, SciLab, Octave, or numpy (:-) ) produces a version for iPad, any chance someone is working on R for iPad? This has come up before and one conclusion was that there may be issues regarding GPL software and the App Store's restrictions, or at least design issues that would make this somewhere between difficult and impossible. For a much more thorough discussion see the threads: https://stat.ethz.ch/pipermail/r-help/2010-May/240669.html https://stat.ethz.ch/pipermail/r-help/2010-June/240901.html -- Brian Diggs Senior Research Associate, Department of Surgery Oregon Health & Science University __ 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] Looking for an image (R 64-bit on Linux 64-bit) on Amazon EC2
I have found an existing image on Amazon EC2 including R. But unfortunately, it is 32-bit R on 32-bit Linux. Does anybody know if there exists an mage (R 64-bit on Linux 64-bit) on Amazon EC2? Or how can I install 64-bit R on my own Linux instance there? Thanks. -- View this message in context: http://r.789695.n4.nabble.com/Looking-for-an-image-R-64-bit-on-Linux-64-bit-on-Amazon-EC2-tp2338938p2338938.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] OT: R for iPhone/iPad OS?
On Aug 25, 2010, at 4:23 PM, Carl Witthoft wrote: > No, seriously: I've had more than one person at work wonder what math > toolset could be loaded onto iOS. So, before Matlab, FreeMat, Mathematia, > SciLab, Octave, or numpy (:-) ) produces a version for iPad, any chance > someone is working on R for iPad? See this coverage of recent discussions at R-Bloggers: http://www.r-bloggers.com/could-we-run-a-statistical-analysis-on-iphoneipad-using-r/ Unless you are willing to jailbreak the devices, the basic answer is not as a traditional stand alone application. However, using a client/server model based GUI application (think Wolfram|Alpha) or via a web browser connecting to a remote R server, yes you can. HTH, Marc Schwartz __ 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] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
Yes, the secant method (like Newton Raphson) is not guaranteed to converge,
unlike the bisection method, but it has a superlinear convergence (not that
this matters much!). Brent's method, which is used in `uniroot', is a
reliable and fast method, which is why I suggested it in my previous email.
Having said that, I am not sure about the convergence problem that you are
having without seeing the actual example.
Ravi.
-Original Message-
From: Adrian Ng [mailto:a...@hamiltonlane.com]
Sent: Wednesday, August 25, 2010 5:28 PM
To: Ravi Varadhan; r-help@r-project.org
Subject: RE: [R] Secant Method Convergence (Method to replicate Excel
XIRR/IRR)
Hi Ravi,
Thanks for the responses. I was actually trying to calculate IRR based on
unevenly spaced cash flows, and that's why I decided to use the secant
method. I'm not sure if my answer isn't converging because I have some
careless mistake in the code, or if it's simply because unlike the bisection
method, the secant method doesn't 'sandwich' the desired root.
-Original Message-
From: Ravi Varadhan [mailto:rvarad...@jhmi.edu]
Sent: Wednesday, August 25, 2010 5:24 PM
To: Adrian Ng; r-help@r-project.org
Subject: RE: [R] Secant Method Convergence (Method to replicate Excel
XIRR/IRR)
Another approach is to use `uniroot' to find the zero of the NPV function:
npv <- function (cashFlow, irr) {
n <- length(cashFlow)
sum(cashFlow / (1 + irr)^{0: (n-1)})
}
uniroot(f=npv, interval=c(0,1), cashFlow=cashFlow)
However, there may be situations where there are no real zeros or there are
multiple zeros of the NPV function.
Ravi.
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Adrian Ng
Sent: Wednesday, August 25, 2010 8:39 AM
To: r-help@r-project.org
Subject: [R] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
Hi,
I am new to R, and as a first exercise, I decided to try to implement an
XIRR function using the secant method. I did a quick search and saw another
posting that used the Bisection method but wanted to see if it was possible
using the secant method.
I would input a Cash Flow and Date vector as well as an initial guess. I
hardcoded today's initial date so I could do checks in Excel. This code
seems to only converge when my initial guess is very close to the correct
IRR.
Maybe I have some basic errors in my coding/logic? Any help would be greatly
appreciated.
The Wikipedia article to secant method and IRR:
http://en.wikipedia.org/wiki/Internal_rate_of_return#Numerical_solution
Thanks!
ANXIRR <- function (cashFlow, cfDate, guess){
cfDate<-as.Date(cfDate,format="%m/%d/%Y")
irrprev <- c(0); irr<- guess
pvPrev<- sum(cashFlow)
pv<-
sum(cashFlow/((1+irr)^(as.numeric(difftime(cfDate,"2010-08-24",units="days")
)/360)))
print(pv)
print("Hi")
while (abs(pv) >= 0.001) {
t<-irrprev; irrprev<- irr;
irr<-irr-((irr-t)*pv/(pv-pvPrev));
pvPrev<-pv;
pv<-sum(cashFlow/((1+irr)^(as.numeric(difftime(cfDate,"2010-08-24",units="da
ys"))/365)))
print(irr);print(pv)
}
}
Please consider the environment before printing this e-mail.
[[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] OT: R for iPhone/iPad OS?
Hi, On Wed, Aug 25, 2010 at 5:37 PM, Greg Snow wrote: > I believe that the iPad conditions for producing apps are incompatible with > the GPL (gnuGo was removed from the apps store recently for this reason), so > don't hold your breath. > > There may be some possibility for an app on the iPad that would access R on a > server somewhere else, but that is beyond my abilities/interest. This has come up a few times before on R-sig-mac. Here's one: http://thread.gmane.org/gmane.comp.lang.r.mac/4912 You can search that list for others discussions that have popped up to get a feel of what was discussed about this already. If you have more to add to the conversation, I reckon you'd get the most traction if you post directly to that mailing list ... -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ 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: R for iPhone/iPad OS?
Instructions for installing R on jailbroken devices: http://rwiki.sciviews.org/doku.php?id=getting-started:installation:iphone Carl Witthoft wrote: No, seriously: I've had more than one person at work wonder what math toolset could be loaded onto iOS. So, before Matlab, FreeMat, Mathematia, SciLab, Octave, or numpy (:-) ) produces a version for iPad, any chance someone is working on R for iPad? __ 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: R for iPhone/iPad OS?
I believe that the iPad conditions for producing apps are incompatible with the GPL (gnuGo was removed from the apps store recently for this reason), so don't hold your breath. There may be some possibility for an app on the iPad that would access R on a server somewhere else, but that is beyond my abilities/interest. -- 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 Carl Witthoft > Sent: Wednesday, August 25, 2010 3:24 PM > To: r-help@r-project.org > Subject: [R] OT: R for iPhone/iPad OS? > > No, seriously: I've had more than one person at work wonder what math > toolset could be loaded onto iOS. So, before Matlab, FreeMat, > Mathematia, SciLab, Octave, or numpy (:-) ) produces a version for > iPad, > any chance someone is working on R for iPad? > > __ > 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] OT: R for iPhone/iPad OS?
No, seriously: I've had more than one person at work wonder what math toolset could be loaded onto iOS. So, before Matlab, FreeMat, Mathematia, SciLab, Octave, or numpy (:-) ) produces a version for iPad, any chance someone is working on R for iPad? __ 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] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
Another approach is to use `uniroot' to find the zero of the NPV function:
npv <- function (cashFlow, irr) {
n <- length(cashFlow)
sum(cashFlow / (1 + irr)^{0: (n-1)})
}
uniroot(f=npv, interval=c(0,1), cashFlow=cashFlow)
However, there may be situations where there are no real zeros or there are
multiple zeros of the NPV function.
Ravi.
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Adrian Ng
Sent: Wednesday, August 25, 2010 8:39 AM
To: r-help@r-project.org
Subject: [R] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
Hi,
I am new to R, and as a first exercise, I decided to try to implement an
XIRR function using the secant method. I did a quick search and saw another
posting that used the Bisection method but wanted to see if it was possible
using the secant method.
I would input a Cash Flow and Date vector as well as an initial guess. I
hardcoded today's initial date so I could do checks in Excel. This code
seems to only converge when my initial guess is very close to the correct
IRR.
Maybe I have some basic errors in my coding/logic? Any help would be greatly
appreciated.
The Wikipedia article to secant method and IRR:
http://en.wikipedia.org/wiki/Internal_rate_of_return#Numerical_solution
Thanks!
ANXIRR <- function (cashFlow, cfDate, guess){
cfDate<-as.Date(cfDate,format="%m/%d/%Y")
irrprev <- c(0); irr<- guess
pvPrev<- sum(cashFlow)
pv<-
sum(cashFlow/((1+irr)^(as.numeric(difftime(cfDate,"2010-08-24",units="days")
)/360)))
print(pv)
print("Hi")
while (abs(pv) >= 0.001) {
t<-irrprev; irrprev<- irr;
irr<-irr-((irr-t)*pv/(pv-pvPrev));
pvPrev<-pv;
pv<-sum(cashFlow/((1+irr)^(as.numeric(difftime(cfDate,"2010-08-24",units="da
ys"))/365)))
print(irr);print(pv)
}
}
Please consider the environment before printing this e-mail.
[[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] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
Hi Adrian,
Your code for the secant method is correct.
I just tweaked it a bit w/o the calendar time feature (assuming that the
cash flows will be available on an annual basis) as follows:
ANXIRR <- function (cashFlow, guess, tol=1.e-04){
npv <- function (cashFlow, irr) {
n <- length(cashFlow)
sum(cashFlow / (1 + irr)^{0: (n-1)})
}
irrprev <- c(0)
irr <- guess
pvPrev <- sum(cashFlow)
pv <- npv(cashFlow, irr)
eps <- abs(pv-pvPrev)
while (eps >= tol) {
tmp <- irrprev
irrprev <- irr
irr <- irr - ((irr - tmp) * pv / (pv - pvPrev))
pvPrev <- pv
pv <- npv(cashFlow, irr)
eps <- abs(pv - pvPrev)
}
list(irr = irr, npv = pv)
}
# example from Wikipedia enntry
cashflow <- c(-4000, 1200, 1410, 1875, 1050)
ANXIRR(cashflow, guess=0.05)
> ANXIRR(cashflow, guess=0.05)
$irr
[1] 0.1429934
$npv
[1] 1.705303e-12
Hope this helps,
Ravi.
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Adrian Ng
Sent: Wednesday, August 25, 2010 8:39 AM
To: r-help@r-project.org
Subject: [R] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
Hi,
I am new to R, and as a first exercise, I decided to try to implement an
XIRR function using the secant method. I did a quick search and saw another
posting that used the Bisection method but wanted to see if it was possible
using the secant method.
I would input a Cash Flow and Date vector as well as an initial guess. I
hardcoded today's initial date so I could do checks in Excel. This code
seems to only converge when my initial guess is very close to the correct
IRR.
Maybe I have some basic errors in my coding/logic? Any help would be greatly
appreciated.
The Wikipedia article to secant method and IRR:
http://en.wikipedia.org/wiki/Internal_rate_of_return#Numerical_solution
Thanks!
ANXIRR <- function (cashFlow, cfDate, guess){
cfDate<-as.Date(cfDate,format="%m/%d/%Y")
irrprev <- c(0); irr<- guess
pvPrev<- sum(cashFlow)
pv<-
sum(cashFlow/((1+irr)^(as.numeric(difftime(cfDate,"2010-08-24",units="days")
)/360)))
print(pv)
print("Hi")
while (abs(pv) >= 0.001) {
t<-irrprev; irrprev<- irr;
irr<-irr-((irr-t)*pv/(pv-pvPrev));
pvPrev<-pv;
pv<-sum(cashFlow/((1+irr)^(as.numeric(difftime(cfDate,"2010-08-24",units="da
ys"))/365)))
print(irr);print(pv)
}
}
Please consider the environment before printing this e-mail.
[[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] package MuMIn
[cc'ing back to r-help: this is good etiquette so that the responses
will be seen by others/
archived for future reference.]
On 10-08-25 04:35 PM, Marino Taussig De Bodonia, Agnese wrote:
> Yes, I meant "MuMIn"
>
> the global formula I introduced was:
>
> rc4.mod<-lm(central$hunting~ central$year + central$gender +
central$hunter + central$k.score + central$seen.in.wild +
central$captivity + central$had.damage +
central$importance.of.being.informed + central$RC1+ central$RC2 +
central$year:central$hunter + central$year:central$had.damage +
central$year:central$seen.in.wild +central$year:central$RC1 +
central$year:central$RC2)
In general it would be much clearer and more generalizable if you
instead used:
rc4.mod<-lm(hunting~ year + gender + hunter + k.score + seen.in.wild +
captivity +
had.damage + importance.of.being.informed + RC1+ RC2 + year:hunter +
year:had.damage +
year:seen.in.wild +year:RC1 + year:RC2, data=central)
or even
rc4.mod <- lm(hunting ~ gender + k.score + captivity +
importance.of.being.informed+
year*(RC1 + RC2 + hunter + had.damage + seen.in.wild), data=central)
I don't know why those spurious interactions are showing up. It
*might* be a bug in MuMIn.
So far, I can't reproduce this behavior with my simple example:
=
library(MuMIn)
dat <- data.frame(y=runif(100),a=factor(sample(1:5,replace=TRUE,size=100)),
b= factor(sample(1:5,replace=TRUE,size=100)),
c=runif(100))
m1 <- lm(y~a+b+c+a:b,data=dat)
gm <- function(x) {
model.avg(get.models(dredge(x),subset=delta
> There are 15 explanatory variables in this model.
>
> my code was:
>
> rc4.mod.Mu <- dredge(rc4.mod, rank = "AICc")
> rc4.mod.Mu
> rc4.mod.Mu.avg<-model.avg(get.models(rc4.mod.Mu, subset = delta < 4))
>
> the output for the command "rc4.mod.Mu.avg$relative.importance" was:
>
> central$hunter1
> central$seen.in.wild1
> central$year0.986960449
> central$had.damage0.89670109
> central$RC20.83866013
> central$k.score0.656654734
> central$had.damage:central$year0.517130185
> central$year:central$hunter0.305988097
> central$hunter:central$year0.212045101
> central$year:central$seen.in.wild0.190520501
> central$seen.in.captivity0.148263242
> central$gender0.119314202
> central$RC10.098234445
> central$importance.of.being.informed0.091842088
> central$year:central$RC20.069501158
> central$seen.in.wild:central$year0.065788243
> central$RC2:central$year0.024221603
>
> There are 17 variables above: central$year:central$RC2 and
central$RC2:central$year are both present, as are
central$hunter:central$year and central$year:central$hunter .
>
> Can you tell me why, if they are the same thing, they are present twice
and have different values?
>
> Thanks a lot,
>
> Agnese
>
> From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] On
Behalf Of Ben Bolker [bbol...@gmail.com]
> Sent: 25 August 2010 20:36
> To: r-h...@stat.math.ethz.ch
> Subject: Re: [R] package MuMI
>
> Marino Taussig De Bodonia, Agnese imperial.ac.uk>
writes:
>
>> I am using the package "MuMI" to run all the possible combinations
>> of variables in my full model, and select
>> my best models. When I enter my variables in the original model I
>> write them like this
>>
>> lm(y~ a +b +c +a:b)
>>
>> However, "MuMI" will also use the variable b:a, which I do not want
>> in my model.
>>
>> How do I stop that from happening?
>>
>
> (1) I think you mean "MuMIn".
> (2) Please send a reproducible example! (Hint: see the posting
> guide that is referred to at the bottom of all R-help messages ...)
> (3) [now proceeding to try to guess what you mean ...]
>
> a:b and b:a have identical meanings in the R formula syntax.
>
> I tried to do something like what I thought you might have meant
> and got what seemed to be reasonable answers.
>
>> library(MuMIn)
>> dat <-
data.frame(y=runif(100),a=factor(sample(1:5,replace=TRUE,size=100)),
> +b= factor(sample(1:5,replace=TRUE,size=100)), c=runif(100))
>> m1 <- lm(y~a+b+c+a:b,data=dat)
>> dredge(m1)
> Global model: lm(formula = y ~ a + b + c + a:b, data = dat)
> ---
> Model selection table
>(Int) a b c a:b k R.sqAdj.R.sq RSS AIC AICc
delta weight
> 4 0.6385 -0.1643 3 0.03019 0.020290 7.386 29.23 29.48
0. 0.420
> 1 0.5530 2 0.0 0.00 7.616 30.29 30.42
0.9392 0.262
> 7 0.6283 + -0.1461 7 0.09397 0.045770 6.900 30.42 31.64
2.1650 0.142
> 3 0.5533 + 6 0.07067 0.031540 7.077 30.96 31.87
2.3890 0.127
> 6 0.6134 + -0.2199 7 0.06700 0.017370 7.105 33.36 34.57
5.0980 0.033
> 2 0.5202 +6 0.01989 -0.021380 7.464 36.28 37.19
7.7100 0.009
> 8 0.6232 + + -0.1941 11 0.12050 0.032540 6.698 35.45 38.45
8.9760 0.005
> 5 0.5352 + + 10 0.08490 0.004455 6.969 37.42 39.89
10.4100 0.002
> 10 0.7433 + + -0.2504 + 27 0.27800 0.034140
Re: [R] Removing inter-bar spaces in barchart
On Aug 24, 2010, at 10:20 PM, Jonathan Greenberg wrote: Rhelpers: I'm trying to make a barchart of a 2-group dataset (barchart(x~y,data=data,groups=z,horizontal=FALSE)). My problem is that I can't, for the life of me, seem to get rid of the inter-bar space -- box.ratio set to 1 doesn't do much. Any ideas? I'd ideally want zero space between the bars. Thanks! You didn't provide any data (nor did you illustrate with one of the available datasets that are used in examples.) Compare these two outputs: barchart(yield ~ year, data = barley, groups = variety, ylab = "Barley Yield (bushels/acre)", ) barchart(yield ~ variety, data = barley, groups = year, ylab = "Barley Yield (bushels/acre)", ) ... and notice that the variables in the "group" have no separation of their bars whereas the rhs variables do. This is the opposite of what I expected, so perhaps you think as I did and reversing the roles of"y" and "z" might help? -- David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] accessing the attr(*,label.table) after importing from spss
Dear all, I just received a file from a colleague in spss. The read.spss could not finish the file due to an error (Unrecognized record type 7, subtype 18 encountered in system file) so instead I converted the file using stat-transfer. Looking at my data I see that most labels are in the attributes and I´d love to access them and assign the pertinent variables to factors without doing the whole factor(levels,labels) thing manually. Is there any remedy to this ?? regards, M > str(dat) 'data.frame': 860 obs. of 19 variables: $ Ag: int 15 15 15 15 15 15 15 15 15 15 ... $ G: int 2 2 2 1 1 1 1 1 1 2 ... $ GCQ: int 15 15 15 15 15 15 15 15 15 15 ... $ Amn: int 2 2 2 2 2 2 1 1 1 1 ... $ HI : int 1 1 1 1 1 1 2 2 2 2 ... $ Hos: int 2 2 2 2 2 2 2 2 2 2 ... $ Risk : int 2 2 2 2 2 2 2 2 2 2 ... $ CTO : int 2 2 2 2 2 1 1 1 1 1 ... $ pat : int NA NA NA NA NA 2 2 2 2 2 ... $ Day: int 7 7 7 5 4 6 5 7 7 5 ... $ Ho : int NA NA NA NA NA NA NA NA NA NA ... $ coh : int 1 1 1 1 1 1 1 1 1 1 ... $ Comp : int 1 1 1 1 1 1 1 1 1 1 ... $ Ethan: int 2 2 2 2 2 2 2 2 2 2 ... $ Pro : num NA NA NA NA NA NA NA NA NA NA ... $ Ye : int 1 1 1 3 3 3 1 1 1 3 ... $ Ageg: int 1 1 1 1 1 1 1 1 1 1 ... $ BAC: int 0 0 0 0 0 0 0 0 0 0 ... - attr(*, "val.labels")= chr "" "VL_Gender" "" "VL_Amnesia" ... - attr(*, "var.labels")= chr "Age (years)" "Gender" "GCQSSA" "Amnesty" ... - attr(*, "label.table")=List of 19 ..$ : NULL ..$ : Named num 1 2 .. ..- attr(*, "names")= chr "Male" "Female" ..$ : NULL ..$ : Named num 1 2 .. ..- attr(*, "names")= chr "Yes" "No" ..$ : NULL ..$ : Named num 1 2 .. ..- attr(*, "names")= chr "Yes" "No" ..$ : Named num 1 2 .. ..- attr(*, "names")= chr "Yes" "No" ..$ : Named num 1 2 .. ..- attr(*, "names")= chr "Yes" "No" ..$ : Named num 1 2 .. ..- attr(*, "names")= chr "Yes" "No" ..$ : Named num 1 2 3 4 5 6 7 .. ..- attr(*, "names")= chr "Monday" "Tuesday" "Wednesday" "Thursday" ... ..$ : Named num 1 2 .. ..- attr(*, "names")= chr "Yes" "No" ..$ : Named num 1 2 .. ..- attr(*, "names")= chr "Yes" "No" ..$ : Named num 1 2 3 4 5 6 7 .. ..- attr(*, "names")= chr "Yes" "Overtriage " "Undertriage with admission" "Overtriage with pos" ... ..$ : Named num 1 2 .. ..- attr(*, "names")= chr "Yes" "No" ..$ : NULL ..$ : NULL ..$ : Named num 1 2 3 4 .. ..- attr(*, "names")= chr "15 - 24" "25-39" "40-59" "60-" ..$ : Named num 1 2 3 .. ..- attr(*, "names")= chr "0.10-0.99" "1.00-1.99" "2.00-" ..$ : Named num 0 1 2 3 4 .. ..- attr(*, "names")= chr "sorry" "undeweight" "0.10-0.99" "1.00-1.99" ... __ 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] Draw a perpendicular line?
Thank you all for your suggestions and especially the codes. Now I am able to solve it by finding the intercepts and slopes for both lines first, and then I can find the coordinates of x and y at the intersection. At first I forgot how to use C to get the intercept of line CD... Thanks again. -- View this message in context: http://r.789695.n4.nabble.com/Draw-a-perpendicular-line-tp2335882p2338864.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] Removing inter-bar spaces in barchart
Hi: On Wed, Aug 25, 2010 at 11:28 AM, Greg Snow wrote: > Using the barplot function in base graphics you just set space=0, but that > function does not have a box.ratio argument which would imply that you are > using something else. If you let us know which function (and which package > it is in) then it is easier (possible) for us to help you, even better is to > give the information asked for at the bottom of every post and the posting > guide. > > It looks like he's using barchart() in lattice, which does have a box.ratio argument. > -- > 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 Jonathan Greenberg > > Sent: Tuesday, August 24, 2010 8:21 PM > > To: r-help > > Subject: [R] Removing inter-bar spaces in barchart > > > > Rhelpers: > > > > I'm trying to make a barchart of a 2-group dataset > > (barchart(x~y,data=data,groups=z,horizontal=FALSE)). My problem is > > that I can't, for the life of me, seem to get rid of the inter-bar > > space -- box.ratio set to 1 doesn't do much. Any ideas? I'd > > ideally want zero space between the bars. Thanks! > I'm more interested in why you'd 'ideally want zero space between the bars'. By default, bar chart functions in R have spaces between bars because the bars represent different categories of a discrete/factor variable and the space is a visual device meant to connote to the viewer that the categories are not to be interpreted as an underlying continuum. In contrast, there are no spaces in between bars of a histogram precisely because the x-scale is meant to be continuous. Since you're using barchart(), why do you want to remove the space between bars? Cheers, Dennis > > > --j > > > > -- > > Jonathan A. Greenberg, PhD > > Assistant Project Scientist > > Center for Spatial Technologies and Remote Sensing (CSTARS) > > Department of Land, Air and Water Resources > > University of California, Davis > > One Shields Avenue > > Davis, CA 95616 > > Phone: 415-763-5476 > > AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307 > > > > __ > > 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. > [[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 do R calculates the number of intervals between tick-marks
Hi,
With axTicks I did what I wanted. Thank you all for the attention.
Below is the code I did to have a plot with three y axes.
Suggestions to improve the routine are welcome.
All the best,
Antonio Olinto
Data in a spreadsheet
YEARL NFB NFT
200526,352158 150025,014
200632,514789 210035,024
200727,458126 180030,254
200828,568971 150043,254
200932,564789 300060,32
# comma is a decimal symbol
dat.lnbnt <- read.delim("clipboard",dec=",",header=T)
attach(dat.lnbnt)
summary(dat.lnbnt)
YEARL NFBNFT
Min. :2005 Min. :26.35 Min. :1500 Min. :25.01
1st Qu.:2006 1st Qu.:27.46 1st Qu.:1500 1st Qu.:30.25
Median :2007 Median :28.57 Median :1800 Median :35.02
Mean :2007 Mean :29.49 Mean :1980 Mean :38.77
3rd Qu.:2008 3rd Qu.:32.51 3rd Qu.:2100 3rd Qu.:43.25
Max. :2009 Max. :32.56 Max. :3000 Max. :60.32
y1min <- 26
y1max <- 33
y1range <- y1max-y1min
y2min <- 1500
y2max <- 3000
y2range <- y2max-y2min
y3min <- 25
y3max <- 61
y3range <- y3max-y3min
par(mar=c(6, 6, 2,12))
plot(L, ylim=c(y1min,y1max), ylab="landings (×1000 t)",
type="l",col="red",las=0, cex.axis=1.2,cex.lab=1.4,xaxt="n",xlab="")
points(y1range*((NFB-y2min)/y2range)+y1min,type="l",col="blue")
points(y1range*((NFT-y3min)/y3range)+y1min,type="l",col="darkgreen")
yticks <- axTicks(2)
interval <- yticks[2]-yticks[1]
divisions <- length(yticks)-1
axis(1,at=axTicks(1),labels=as.character(YEAR),las=2,cex.axis=1.2)
axis(4,at=seq(y1min,y1max,interval),labels=as.integer(seq(y2min,y2max,y2range/divisions)),las=0,cex.axis=1.2)
axis(4,at=seq(y1min,y1max,interval),labels=as.integer(seq(y3min,y3max,y3range/divisions)),cex.axis=1.2,
las=0,line=5)
mtext("nº of fishing boats",4,3,cex=1.4)
mtext("nº of fishing trips (×1000)",4,8, cex=1.4)
legend(4,33,c("L","Nº FB","Nº
FT"),bty="n",lty=c(1,1,1),col=c("red","blue","darkgreen"),cex=1.4)
Citando David Winsemius :
On Aug 24, 2010, at 7:05 PM, Antonio Olinto wrote:
Hello,
I want to know how do R calculates the number of intervals between
tick-marks in the y axis in a plot.
?axTicks # and then look at the other citations and the code as needed
I'm making a three y-axes plot and this information would help me a lot.
Thanks in advance.
Antonio Olinto
--
David Winsemius, MD
West Hartford, CT
Webmail - iBCMG Internet
http://www.ibcmg.com.br
__
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 does this warning message (from optim function) mean?
the deteminant is a nonpositive value. log(det(...)) produce NaNs... -- View this message in context: http://r.789695.n4.nabble.com/What-does-this-warning-message-from-optim-function-mean-tp2338689p2338719.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] Checking a point behind two polygons
Hi, I am working on a problem of checking whether a point is behind two convex hulls. By 'behind', I mean a point which is outside of two convex hulls can't reach either of these two hulls without reaching the other. My idea is to find the segments between the point and the vertices of a hull first. And then check whether there is any intersection between the segments and the edges of the other hull. If for all the edges to be checked, there is no intersection point found outside of the edges --> the point is said to be behind that hull. Basically, I am doing line-line intersection. I know if I use the line-line intersection checking, I will have to deal with lots of edges and segments. Besides, the line-line intersection checking may not be working...I wonder someone knows a simply way and more accurate way to check it, e.g.line-polygon intersection. Thank you. -- View this message in context: http://r.789695.n4.nabble.com/Checking-a-point-behind-two-polygons-tp2338823p2338823.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] What does this warning message (from optim function) mean?
Hi,
You did not give us any information about your likelihood function, f, nor
did you provide a reproducible example. So, I cannot tell for sure whether
the parameter estimates are reliable.
Ravi.
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Sally Luo
Sent: Wednesday, August 25, 2010 11:26 AM
To: r-help@r-project.org
Subject: [R] What does this warning message (from optim function) mean?
Hi R users,
I am trying to use the optim function to maximize a likelihood funciton, and
I got the following warning messages.
Could anyone explain to me what messege 31 means exactly? Is it a cause for
concern?
Since the value of convergence turns out to be zero, it means that the
converging is successful, right?
So can I assume that the parameter estimates generated thereafter are
reliable MLE estimates?
Thanks a lot for your help.
Maomao
> p<-optim(c(0,0,0), f, method ="BFGS", hessian =T, y=y,X=X,W=W)
There were 31 warnings (use warnings() to see them)
> warnings()
Warning messages:
1: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
2: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
3: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
4: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
5: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
6: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
7: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
8: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
9: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
10: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
11: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
12: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
13: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
14: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
15: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
16: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
17: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
18: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
19: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
20: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
21: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
22: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
23: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
24: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
25: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
26: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
27: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
28: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
29: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
30: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
31: In if (hessian) { ... :
the condition has length > 1 and only the first element will be used
> p$counts
function gradient
148 17
> p$convergence
[1] 0
[[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] Draw a perpendicular line?
Hi:
The function below plots the line segment between two points A and B as well
as the normal from C to AB, with a dot as the intersection point (D), which
is returned with the function call. The aspect ratio is kept at one so that
the orthogonality between the two lines is not distorted by the axis
scaling.
The input is a 3 x 2 matrix with point A in the first row, B in the second
and C in the third.
normalLine <- function(pts) {
# A function that takes points A, B, C (rows of matrix pts)
# as input
# find slope of line between A and B
# find equation of normal that passes through C
# find the intersection point D between AB and its normal
# through C
if(nrow(pts) != 3 || ncol(pts) != 2) stop('input matrix must be 3 x 2')
A <- pts[1, , drop = TRUE]
B <- pts[2, , drop = TRUE]
C <- pts[3, , drop = TRUE]
slopeAB <- (B[2] - A[2])/(B[1] - A[1])
slopeNorm <- -1/slopeAB
a <- A[2] - slopeAB * A[1]
b <- C[2] - slopeNorm * C[1]
xintersect <- (b - a)/(slopeAB - slopeNorm)
yintersect <- b + slopeNorm * xintersect
mxlim <- max(pts)
mnlim <- min(pts)
plot(pts[, 1], pts[, 2], pch = c('A', 'B', 'C'), xlab = "", ylab = "",
xlim = c(mnlim, mxlim), ylim = c(mnlim, mxlim))
segments(A[1], A[2], B[1], B[2])
segments(xintersect, yintersect, C[1], C[2])
points(xintersect, yintersect, pch = 16, cex = 1.2)
c(xintersect, yintersect)
}
## Test:
pts <- matrix(c(1, 2, 5, 9, 3, 4), ncol = 2, byrow = TRUE)
normalLine(pts)
HTH,
Dennis
On Wed, Aug 25, 2010 at 8:54 AM, William Revelle wrote:
> At 3:04 PM -0700 8/23/10, CZ wrote:
>
>> Hi,
>>
>> I am trying to draw a perpendicular line from a point to two points.
>> Mathematically I know how to do it, but to program it, I encounter some
>> problem and hope can get help. Thanks.
>>
>> I have points, A, B and C. I calculate the slope and intercept for line
>> drawn between A and B.
>> I am trying to check whether I can draw a perpendicular line from C to
>> line
>> AB and get the x,y value for the point D at the intersection.
>>
>> Assume I get the slope of the perpendicular line, I will have my point (D)
>> using variable x and y which is potentially on line AB. My idea was
>> using
>> |AC|*|AC| = |AD|*|AD|+ |CD|*|CD|. I don't know what function I may need
>> to
>> call to calculate the values for point D (uniroot?).
>>
>
>
> This is easier than you think.
> Think of the x,y coordinates of each point :
> Then, the slope is slope = rise/run = (By- Ay)/(Bx- Ax)
> The Dx coordinate = Cx and the Dy = (Dx - Ax) * slope
> Then, to draw the line segment from C to D
> lines(C,D)
>
> In R:
>
> A <- c(2,4)
> B <- c(4,1)
> C <- c(8,10)
> slope <-( C[2]- A[2])/(C[1]-A[1]) #rise/run
> D <- c(B[1],(B[1]-A[1])*slope + A[2]) # find D
>
> my.data <- rbind(A,B,C,D)
> colnames(my.data) <- c("X","Y")
> my.data#show it
> plot(my.data,type="n") #graph it without the points
> text(my.data,rownames(my.data)) #put the points in
> segments(A[1],A[2],C[1],C[2]) #draw the line segments
> segments(B[1],B[2],D[1],D[2])
>
> Bill
>
>
>
>
>> Thank you.
>>
>>
>>
>> --
>> View this message in context:
>> http://r.789695.n4.nabble.com/Draw-a-perpendicular-line-tp2335882p2335882.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.
>
[[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] Change value of a slot of an S4 object within a method.
Howdy, On Wed, Aug 25, 2010 at 1:21 PM, Joris Meys wrote: > Hi Steve, > > thanks for the tip. I'll definitely take a closer look at your > solution for implementation for future use. But right now I don't > have the time to start rewriting my class definitions. > > Luckily, I found where exactly things were going wrong. After reading > into the documentation about the evaluation in R, I figured out I have > to specify the environment where substitute should look explicitly as > parent.frame(1). I still don't understand completely why exactly, but > it does the job. > > Thus : > eval(eval(substitute(expression(obj...@extra[[name]] <<- value > > should become : > > eval( > eval( > substitute( > expression(obj...@extra[[name]] <<- value) > ,env=parent.frame(1) ) > ) > ) My eyes start to gloss over on their first encounter of `substitute` ... add enough `eval`s and (even) an `expression` (!) to that, and you'll probably see me running for the hills ... but hey, if it makes sense to you, more power to you ;-) Glad you found a fix that works, -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ 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] package MuMI
Marino Taussig De Bodonia, Agnese imperial.ac.uk> writes: > I am using the package "MuMI" to run all the possible combinations > of variables in my full model, and select > my best models. When I enter my variables in the original model I > write them like this > > lm(y~ a +b +c +a:b) > > However, "MuMI" will also use the variable b:a, which I do not want > in my model. > > How do I stop that from happening? > (1) I think you mean "MuMIn". (2) Please send a reproducible example! (Hint: see the posting guide that is referred to at the bottom of all R-help messages ...) (3) [now proceeding to try to guess what you mean ...] a:b and b:a have identical meanings in the R formula syntax. I tried to do something like what I thought you might have meant and got what seemed to be reasonable answers. > library(MuMIn) > dat <- data.frame(y=runif(100),a=factor(sample(1:5,replace=TRUE,size=100)), +b= factor(sample(1:5,replace=TRUE,size=100)), c=runif(100)) > m1 <- lm(y~a+b+c+a:b,data=dat) > dredge(m1) Global model: lm(formula = y ~ a + b + c + a:b, data = dat) --- Model selection table (Int) a b c a:b k R.sqAdj.R.sq RSS AIC AICc delta weight 4 0.6385 -0.1643 3 0.03019 0.020290 7.386 29.23 29.48 0. 0.420 1 0.5530 2 0.0 0.00 7.616 30.29 30.42 0.9392 0.262 7 0.6283 + -0.1461 7 0.09397 0.045770 6.900 30.42 31.64 2.1650 0.142 3 0.5533 + 6 0.07067 0.031540 7.077 30.96 31.87 2.3890 0.127 6 0.6134 + -0.2199 7 0.06700 0.017370 7.105 33.36 34.57 5.0980 0.033 2 0.5202 +6 0.01989 -0.021380 7.464 36.28 37.19 7.7100 0.009 8 0.6232 + + -0.1941 11 0.12050 0.032540 6.698 35.45 38.45 8.9760 0.005 5 0.5352 + + 10 0.08490 0.004455 6.969 37.42 39.89 10.4100 0.002 10 0.7433 + + -0.2504 + 27 0.27800 0.034140 5.498 47.71 68.71 39.2400 0.000 9 0.6193 + + + 26 0.23150 -0.014420 5.853 51.96 71.19 41.7200 0.000 __ 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 odfWeave questions
Ben,
> 1a. am I right in believing that odfWeave does not respect the
> 'keep.source' option? Am I missing something obvious?
I believe it does, since this gets passed directly to Sweave.
> 1b. is there a way to set global options analogous to \SweaveOpts{}
> directives in Sweave? (I looked at odfWeaveControl, it doesn't seem to
> do it.)
Yes. There are examples of this in the 'examples' package directory.
> 2. I tried to write a Makefile directive to process files from the
> command line:
>
> %.odt: %_in.odt
> $(RSCRIPT) -e "library(odfWeave);
> odfWeave(\"$*_in.odt\",\"$*.odt\");"
>
> This works, *but* the resulting output file gives a warning ("The file
> 'odftest2.odt' is corrupt and therefore cannot be opened.
> OpenOffice.org can try to repair the file ..."). Based on looking at
> the contents, it seems that a spurious/unnecessary 'Rplots.pdf' file is
> getting
> created and zipped in with the rest of the archive; when I unzip, delete
> the Rplots.pdf file and re-zip, the ODT file opens without a warning.
> Obviously I could post-process but it would be nicer to find a
> workaround within R ...
Get the latest version form R-Forge. I haven't gotten this fix onto
CRAN yet (I've been on a caret streak lately).
> 3. I find the requirement that all file paths be specified as absolute
> rather than relative paths somewhat annoying -- I understand the reason,
> but it goes against one practice that I try to encourage for
> reproducibility, which is *not* to use absolute file paths -- when
> moving a same set of data and analysis files across computers, it's hard
> to enforce them all ending up in the same absolute location, which then
> means that the recipient has to edit the ODT file. It would be nice if
> there were hooks for read.table() and load() as there are for plotting
> and package/namespace loading -- then one could just copy them into the
> working directory on the fly.
> has anyone experienced this/thought of any workarounds?
> (I guess one solution is to zip any necessary source files into the archive
> beforehand,
> as illustrated in the vignette.)
You can set the working directory with the (wait for it...) 'workDir'
argument. Using 'workDir = getwd()' will pack and unpack the files in
the current location and you wouldn't need to worry about setting the
path. I use the temp directory because I started over-wrting files.
Max
__
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] GLM outputs in condensed versus expanded table
On 08/25/2010 05:17 PM, francogrex wrote: > > Hi I'm having different outputs from GLM when using a "condensed" table > V1V2 V3 Present Absent > 0 0 0 3 12 > 0 0 1 0 0 > 0 1 0 0 0 > 0 1 1 1 0 > 1 0 0 7 20 > 1 0 1 0 0 > 1 1 0 3 0 > 1 1 1 6 0 > > > resp=cbind(Present, Absent) > glm(resp~V1+V2+V3+I(V1*V2*V3),family=binomial) >> Deviance Residuals: > [1] 0 0 0 0 0 0 0 0 > etc and also coefficients... > > And when using the same but "expanded" table > > V1 V2 V3 condition (1 present 0 abscent) > Id1 1 0 0 1 > id2 1 1 1 1 > ... etc > glm(condition~V1+V2+V3+I(V1*V2*V3),family=binomial) >> Deviance Residuals: > Min 1Q Median 3Q Max > -0.7747317 -0.7747317 -0.6680472 0.0001315 1.7941226 > and also coefficients are different from above. > > What could I be doing wrong? > > Not necessarily anything. Anything technical, that is. You have 3 uninformative combinations where the total is zero. The model has 5 parameters. This is quite likely to generate a perfect fit to the aggregated data. With the groups having zeros in the "absent" category, the fit probably diverged so some coefficients are numerically large. Refitting with individual data will likely give slightly different coefficients, since it sort of depends on how far you came on the way to infinity. With the aggregated data, a perfect fit gives residuals of zero, but with individual data, the 0's and 1's give negative and positive residuals. Try z <- rep(0:1,5) zz <- cbind(5,5) summary(glm(z~1, binomial)) summary(glm(zz~1, binomial)) -- 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.
Re: [R] lattice help required
.. thanks again, richard.
and you swiftly saw the next problem comming up - when using par.settings =
list(layout.widths = list(axis.panel = c(1, 0))) getting rid of the double
tick labeling would be natural -
but i'll leave it at that for today.
many thanks,
kay
Richard M. Heiberger wrote:
>
> The multiple y axes are protecting you in this situation.
>
>
> z <- cbind(rnorm(100,c(1,10),1), rnorm(100,c(20,30),1))
> dotplot(z[,1]+z[,2] ~ facs$Treatment|facs$Sites,
> outer=TRUE,
> scales = list(
> y = list(
> relation="free")),
> ylab=c("y1", "y2"),
> xlab=c("Site 1", "Site 2"),
> strip=FALSE,
> main="problem")
>
> dotplot(z[,1]+z[,2] ~ facs$Treatment|facs$Sites,
> outer=TRUE,
> scales = list(
> y = list(
> relation="free",
> limits=list(c(-5,13),c(-5,13),c(18,32),c(18,32,
> ylab=c("y1", "y2"),
> xlab=c("Site 1", "Site 2"),
> strip=FALSE, main="protected")
>
> For more control (such as suppressing the y-tick labels in the right-hand
> column,
> I recommend Deepayan Sarkar's book.
>
> Rich
>
> [[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.
>
>
-
Kay Cichini
Postgraduate student
Institute of Botany
Univ. of Innsbruck
--
View this message in context:
http://r.789695.n4.nabble.com/lattice-help-required-tp2338382p2338707.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] approxfun-problems (yleft and yright ignored)
The plots did not come through, see the posting guide for which attachments are
allowed. It will be easier for us to help if you can send reproducible code
(we can copy and paste to run, then examine, edit, etc.). Try finding a subset
of your data for which the problem still occurs, then send the data if
possible, or similar simulated data if you cannot send original data.
--
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 Samuel Wuest
> Sent: Wednesday, August 25, 2010 8:20 AM
> To: r-help@r-project.org
> Subject: [R] approxfun-problems (yleft and yright ignored)
>
> Dear all,
>
> I have run into a problem when running some code implemented in the
> Bioconductor panp-package (applied to my own expression data), whereby
> gene
> expression values of known true negative probesets (x) are interpolated
> onto
> present/absent p-values (y) between 0 and 1 using the *approxfun -
> function*{stats}; when I have used R version 2.8, everything had
> worked fine,
> however, after updating to R 2.11.1., I got unexpected output
> (explained
> below).
>
> Please correct me here, but as far as I understand, the yleft and
> yright
> arguments set the extreme values of the interpolated y-values in case
> the
> input x-values (on whose approxfun is applied) fall outside range(x).
> So if
> I run approxfun with yleft=1 and yright=0 with y-values between 0 and
> 1,
> then I should never get any values higher than 1. However, this is not
> the
> case, as this code-example illustrates:
>
> > ### define the x-values used to construct the approxfun, basically
> these
> are 2000 expression values ranging from ~ 3 to 7:
> > xNeg <- NegExprs[, 1]
> > xNeg <- sort(xNeg, decreasing = TRUE)
> >
> > ### generate 2000 y-values between 0 and 1:
> > yNeg <- seq(0, 1, 1/(length(xNeg) - 1))
> > ### define yleft and yright as well as the rule to clarify what
> should
> happen if input x-values lie outside range(x):
> > interp <- approxfun(xNeg, yNeg, yleft = 1, yright = 0, rule=2)
> Warning message:
> In approxfun(xNeg, yNeg, yleft = 1, yright = 0, rule = 2) :
> collapsing to unique 'x' values
> > ### apply the approxfun to expression data that range from ~2.9 to
> 13.9
> and can therefore lie outside range(xNeg):
> > PV <- sapply(AllExprs[, 1], interp)
> > range(PV)
> [1]0.000 6208.932
> > summary(PV)
> Min. 1st Qu.Median Mean 3rd Qu. Max.
> 0.000e+00 0.000e+00 2.774e-03 1.299e+00 3.164e-01 6.209e+03
>
> So the resulting output PV object contains data ranging from 0 to 6208,
> the
> latter of which lies outside yleft and is not anywhere close to extreme
> y-values that were used to set up the interp-function. This seems wrong
> to
> me, and from what I understand, yleft and yright are simply ignored?
>
> I have attached a few histograms that visualize the data distributions
> of
> the objects I xNeg, yNeg, AllExprs[,1] (== input x-values) and PV (the
> output), so that it is easier to make sense of the data structures...
>
> Does anyone have an explanation for this or can tell me how to fix the
> problem?
>
> Thanks a million for any help, best, Sam
>
> > sessionInfo()
> R version 2.11.1 (2010-05-31)
> x86_64-apple-darwin9.8.0
>
> locale:
> [1] en_IE.UTF-8/en_IE.UTF-8/C/C/en_IE.UTF-8/en_IE.UTF-8
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] panp_1.18.0 affy_1.26.1 Biobase_2.8.0
>
> loaded via a namespace (and not attached):
> [1] affyio_1.16.0 preprocessCore_1.10.0
>
>
> --
> -
> Samuel Wuest
> Smurfit Institute of Genetics
> Trinity College Dublin
> Dublin 2, Ireland
> Phone: +353-1-896 2444
> Web: http://www.tcd.ie/Genetics/wellmer-2/index.html
> Email: wue...@tcd.ie
> --
__
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 does this warning message (from optim function) mean?
You mean 'TRUE': 'T' is a variable in R, with initial value TRUE.
On Wed, 25 Aug 2010, Sally Luo wrote:
Hi R users,
I am trying to use the optim function to maximize a likelihood funciton, and
I got the following warning messages.
Could anyone explain to me what messege 31 means exactly? Is it a cause for
concern?
Since the value of convergence turns out to be zero, it means that the
converging is successful, right?
So can I assume that the parameter estimates generated thereafter are
reliable MLE estimates?
Thanks a lot for your help.
Maomao
p<-optim(c(0,0,0), f, method ="BFGS", hessian =T, y=y,X=X,W=W)
There were 31 warnings (use warnings() to see them)
warnings()
Warning messages:
1: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
2: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
3: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
4: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
5: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
6: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
7: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
8: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
9: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
10: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
11: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
12: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
13: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
14: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
15: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
16: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
17: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
18: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
19: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
20: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
21: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
22: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
23: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
24: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
25: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
26: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
27: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
28: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
29: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
30: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
31: In if (hessian) { ... :
the condition has length > 1 and only the first element will be used
p$counts
function gradient
148 17
p$convergence
[1] 0
[[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.
--
Brian D. Ripley, rip...@stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax: +44 1865 272595
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Removing inter-bar spaces in barchart
Using the barplot function in base graphics you just set space=0, but that function does not have a box.ratio argument which would imply that you are using something else. If you let us know which function (and which package it is in) then it is easier (possible) for us to help you, even better is to give the information asked for at the bottom of every post and the posting guide. -- 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 Jonathan Greenberg > Sent: Tuesday, August 24, 2010 8:21 PM > To: r-help > Subject: [R] Removing inter-bar spaces in barchart > > Rhelpers: > > I'm trying to make a barchart of a 2-group dataset > (barchart(x~y,data=data,groups=z,horizontal=FALSE)). My problem is > that I can't, for the life of me, seem to get rid of the inter-bar > space -- box.ratio set to 1 doesn't do much. Any ideas? I'd > ideally want zero space between the bars. Thanks! > > --j > > -- > Jonathan A. Greenberg, PhD > Assistant Project Scientist > Center for Spatial Technologies and Remote Sensing (CSTARS) > Department of Land, Air and Water Resources > University of California, Davis > One Shields Avenue > Davis, CA 95616 > Phone: 415-763-5476 > AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307 > > __ > 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] output values from within a function
Put the line:
cat(z1,'\n')
in your function. You may also want to put the flush.console() command right
after that.
Hope this helps,
--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111
> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> project.org] On Behalf Of Maas James Dr (MED)
> Sent: Friday, August 20, 2010 7:41 AM
> To: r-help@r-project.org
> Subject: [R] output values from within a function
>
> Is it possible to get R to output the value of an expression, that is
> being calculated within a function? I've attached a very simple
> example but for more complicated ones would like to be able to debug by
> seeing what the value of specific expressions are each time it cycles
> through a loop that executes the expression. I'm relatively new to
> this so there may be much simpler more elegant ways to do it.
>
> For example, is there a command I can put within the function "funct01"
> that will output the value of z1 to the screen?
>
> Thanks
>
> Jim
>
> ==
> ## Practice file to try out evaluations
>
> y <- 5
>
> (x <- y^2)
>
> funct01 <- function (x,y) {
>
> z1 <- x + y
> z2 <- x * y
> z3 <- x^y
>
> results <- data.frame (z1=z1, z2=z2, z3=z3)
> return(results)
>
> }
>
> funct01(7,9)
>
>
> ===
> Dr. Jim Maas
> University of East Anglia
> Norwich, UK
> NR4 7TJ
>
> [[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] SEM : Warning : Could not compute QR decomposition of Hessian
Dear Anne,
I started to diagram your model but stopped when I noticed some problems:
(1) Some variables, such as pays_alti, are clearly endogenous, since they
have arrows pointing to them, yet are declared as fixed exogenous variables;
that clearly doesn't make sense. (2) You've placed conflicting constraints
on factor loadings and the variances of latent variables, for example
setting both the path from the factor type_paysage to the indicator
pays_alti and the variance of type_paysage to 1; again, that doesn't make
sense.
Do you have a path diagram for the model that you're trying to fit? It would
in particular be helpful to have a diagram of the structural part of the
model.
Regards,
John
John Fox
Senator William McMaster
Professor of Social Statistics
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox
> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
On
> Behalf Of Anne Mimet
> Sent: August-25-10 10:27 AM
> To: r-help@r-project.org
> Subject: [R] SEM : Warning : Could not compute QR decomposition of Hessian
>
> Hi useRs,
>
> I'm trying for the first time to use a sem. The model finally runs,
> but gives a warning saying :
> "In sem.default(ram = ram, S = S, N = N, param.names = pars, var.names
> = vars, : Could not compute QR decomposition of Hessian.
> Optimization probably did not converge. "
>
> I found in R-help some posts on this warning, but my attemps to modify
> the code didn't change the warning message (i tried to give an error
> of 1 to the latente variables). I can't figure what the problem is.
> Here is the code :
>
>
> tab<-read.table("F:/Mes documents/stats/sem/donnees_corr.txt",
> header=T, sep="",na.strings = "NA")
>
> tab[,46]<-as.factor(tab[,46])
> tab[,24]<-as.factor(tab[,24])
> tab[,40]<-as.factor(tab[,40])
>
> fct_cor<-hetcor(tab, ML=T)
> cor_tab<- fct_cor$correlations
>
> moment_tab<-read.moments(diag=F, names=c('c1','c2', 'c3','c4','c5',
> 'c6','c7', 'c8', 'c9', 'ind_plando', 'long_sup15', 'long_inf15',
> 'pente', 'est', 'sud','ouest', 'nord' ,'reg_hydriq', 'prof_sol',
> 'pierro', 'efferv', 'struct','drainage','texture', 'route1_pond',
> 'route2_pond',
> 'pourcactif', 'tx_chomage', 'pourcagric', 'pourc_jeunes', 'pop99',
> 'rev_imp_foyer','eq_CONC', 'eq_sante', 'eq_edu', 'sold_nat',
> 'sold_mig', 'tps_dom_emp','TXEMPLOI','ORIECO','dist_paris','axe1',
> 'axe2', 'axe3', 'dist_protect','urbanisation','pays_incli','pays_alti'))
> # after comes the moment matrix (triangular)
>
> ram_tab<-specify.model()
> type_paysage->pays_alti,NA,1
> type_paysage->pays_incli, pays2, NA
> pedo->reg_hydriq, NA, 1
> pedo->prof_sol, ped8, NA
> pedo->pierro, ped9, NA
> pedo->efferv, ped10, NA
> pedo->struct, ped11, NA
> pedo->drainage, ped12, NA
> pedo->texture, ped13, NA
> adj_99->c1, NA,1
> adj_99->c2, adj2,NA
> adj_99->c3, adj3,NA
> adj_99->c4, adj4,NA
> adj_99->c5, adj5,NA
> adj_99->c6, adj6,NA
> adj_99->c7, adj7,NA
> adj_99->c8, adj8,NA
> adj_99->c9, adj9,NA
> etat_hexa->axe1, NA, 1
> etat_hexa->axe2, et2, NA
> etat_hexa->axe3, et3, NA
> socioBV->sold_mig, BV1, NA
> socioBV->sold_nat, BV2, NA
> socioBV->TXEMPLOI, BV3, NA
> socioBV->ORIECO, BV4, NA
> socioBV->tps_dom_emp, NA, 1
> eqBV->eq_CONC, NA, 1
> eqBV->eq_sante, eq2, NA
> eqBV->eq_edu, eq3, NA
> socio_com->pourcactif , NA, 1
> socio_com->tx_chomage, com2, NA
> socio_com->pourcagric, com3, NA
> socio_com->pourc_jeunes, com4, NA
> socio_com->pop99, com5, NA
> socio_com->rev_imp_foyer, com7, NA
> access_hexa->route1_pond, NA, 1
> access_hexa->route2_pond, acc2, NA
> hydro->ind_plando, NA, 1
> hydro->long_sup15, eau2, NA
> hydro->long_inf15, eau3, NA
> topog->pente, NA, 1
> topog->est, top2, NA
> topog->sud, top3, NA
> topog->nord, top4, NA
> topog->ouest, top5, NA
> dist_protect-> urbanisation, cor1,NA
> dist_protect-> adj_99, cor2, NA
> dist_protect-> etat_hexa, cor3, NA
> topog-> urbanisation, cor4, NA
> topog-> adj_99, cor5, NA
> topog-> etat_hexa, cor6, NA
> topog-> access_hexa, cor7, NA
> topog<->hydro, cor8, NA
> topog<->pedo, cor9, NA
> pedo-> urbanisation, cor10, NA
> pedo-> adj_99, cor11, NA
> pedo-> etat_hexa, cor12, NA
> pedo<->hydro, cor1, NA
> hydro-> urbanisation, cor13, NA
> hydro-> adj_99, cor14, NA
> hydro-> etat_hexa, cor15, NA
> access_hexa-> urbanisation, cor16, NA
> access_hexa-> etat_hexa, cor17, NA
> socio_com-> etat_hexa, cor18, NA
> socio_com-> adj_99, cor19, NA
> socio_com-> urbanisation, cor20, NA
> dist_paris-> socio_com, cor21, NA
> dist_paris-> access_hexa, cor22, NA
> dist_paris-> adj_99, cor23, NA
> dist_paris-> etat_hexa, cor24, NA
> dist_paris-> urbanisation, cor25, NA
> dist_paris-> socioBV, cor26, NA
> socioBV-> eqBV, cor27, NA
> socioBV-> urbanisation, cor28, NA
> socioBV-> adj_99, cor29, NA
> socioBV-> etat_hexa, cor30, NA
> eqBV-> etat_hexa, cor31, NA
> eqBV-> adj_99, cor32, NA
> eqBV-> urbanisation, cor33, NA
> etat_hexa-> urbanisation, cor34, NA
> eta
[R] [R-pkgs] stringr: version 0.4
Strings are not glamorous, high-profile components of R, but they do play a big role in many data cleaning and preparations tasks. R provides a solid set of string operations, but because they have grown organically over time, they can be inconsistent and a little hard to learn. Additionally, they lag behind the string operations in other programming languages, so that some things that are easy to do in languages like Ruby or Python are rather hard to do in R. The `stringr` package aims to remedy these problems by providing a clean, modern interface to common string operations. More concretely, `stringr`: * Processes factors and characters in the same way. * Gives functions consistent names and arguments. * Simplifies string operations by eliminating options that you don't need 95% of the time. * Produces outputs than can easily be used as inputs. This includes ensuring that missing inputs result in missing outputs, and zero length inputs result in zero length outputs. * Completes R's string handling functions with useful functions from other programming languages. New in stringr 0.4: * all functions now vectorised with respect to string, pattern (and where appropriate) replacement parameters * fixed() function now tells stringr functions to use fixed matching, rather than escaping the regular expression. Should improve performance for large vectors. * new ignore.case() modifier tells stringr functions to ignore case of pattern. * str_replace renamed to str_replace_all and new str_replace function added. This makes str_replace consistent with all functions. * new str_sub<- function (analogous to substring<-) for substring replacement * str_sub now understands negative positions as a position from the end of the string. -1 replaces Inf as indicator for string end. * str_pad side argument can be left, right, or both (instead of center) * str_trim gains side argument to better match str_pad * stringr now has a namespace and imports plyr (rather than requiring it) -- Assistant Professor / Dobelman Family Junior Chair Department of Statistics / Rice University http://had.co.nz/ ___ R-packages mailing list r-packa...@r-project.org https://stat.ethz.ch/mailman/listinfo/r-packages __ 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] multivariate analysis of variance
Hi, can anyone tell me how i may find out the between-group variance and within-group variance for multivariate case? I have 6 Groups and 73 Variables. (with MANOVA ? wie) dim(data) [1] 2034 76 Thanks Celal __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problem with clusterCall, "Error in checkForRemoteErrors(lapply(cl, recvResult)) : "
Oh, I forgot to say, ui.Next() is a function I defined to sample from the proposal distribution given the current state. -- View this message in context: http://r.789695.n4.nabble.com/Problem-with-clusterCall-Error-in-checkForRemoteErrors-lapply-cl-recvResult-tp2338375p2338378.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] pairwise correlation of many samples
Dear list, I have a csv file like below, and I have similar 150 files (responses from 150 samples) in a folder. Id like to get the mean score of pairwise correlation score among 150 respondents by hour(0~20hour). All I can think of is bring up two files and get the pairwise correlation score and repeating this (e.g. (p1,p2),(p1,P3) (P4,P5)..)and get the total mean of correlation score. However, since therere 150 samples I cant think of doing this manually.. Hour Response 0 -0.11703 0.016667 -0.10427 0.03 -0.091347 0.05 -0.078313 0.07 -0.065226 0.08 -0.052147 0.1 -0.039136 . . Can anyone please point out the best way to do? .. Its been only a month since I started to use R Thank you! [[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] What does this warning message (from optim function) mean?
Hi R users,
I am trying to use the optim function to maximize a likelihood funciton, and
I got the following warning messages.
Could anyone explain to me what messege 31 means exactly? Is it a cause for
concern?
Since the value of convergence turns out to be zero, it means that the
converging is successful, right?
So can I assume that the parameter estimates generated thereafter are
reliable MLE estimates?
Thanks a lot for your help.
Maomao
> p<-optim(c(0,0,0), f, method ="BFGS", hessian =T, y=y,X=X,W=W)
There were 31 warnings (use warnings() to see them)
> warnings()
Warning messages:
1: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
2: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
3: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
4: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
5: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
6: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
7: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
8: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
9: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
10: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
11: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
12: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
13: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
14: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
15: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
16: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
17: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
18: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
19: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
20: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
21: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
22: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
23: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
24: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
25: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
26: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
27: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
28: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
29: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
30: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced
31: In if (hessian) { ... :
the condition has length > 1 and only the first element will be used
> p$counts
function gradient
148 17
> p$convergence
[1] 0
[[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] Problem with clusterCall, "Error in checkForRemoteErrors(lapply(cl, recvResult)) : "
Hi all,
I am trying to use snow package to do a parallel MCMC. I have read a few
guides and articles, the following is that I came up with.
When I run it I got the error message:
Error in checkForRemoteErrors(lapply(cl, recvResult)) :
4 nodes produced errors; first error: could not find function "ui.Next"
The data is a longitudinal data with few repeated readings on a number of
individuals. However, the response is organised in a vector rather than a
matrix. E.g. (y_11 , y_12 , y_13 , y_14 , y_21 , y_22 , ... , y_n1 , ... ,
y_nTn )^T
X , Z are covariates.
beta is a matrix of coefficients associated with X
C is the latent class membership indicator
sigma.sq is the diagonal elements of the covariance matrix of u, which is a
random effect parameter and the one I am trying to sample with the function.
sd is the diagonal elements of the covariance matrix of the proposal
distribution.
The following is the particular function:
ui.Full.Sample.Cluster = function( levels , Y , X , beta , Z , C , sigma.sq
, sd , burnin , iteration )
{
cluster = makeCluster( 4 , type = "MPI");
clusterSetupRNG( cluster );
patients = unique( levels );
q = length( X[1,] );
u = ui.Ini( q , length( patients ) );
n = levels[ length(patients) ];
marker = levels == patients[1];
y = Y[ marker ];
x = X[ marker, ];
z = Z[ marker, ];
u[1,] = ui.1.Sample( u[1,] , y , x , beta , z , C[1] , sigma.sq , sd ,
burnin , iteration )$uFinal;
for( i in 2:n)
{
marker = levels == patients[i];
y = Y[ marker ];
x = X[ marker, ];
z = Z[ marker, ];
print( i );
u[i, ] = clusterCall( cluster , ui.1.Sample, u[i,] , y , x ,
beta , z ,
C[i] , sigma.sq , sd , burnin , iteration )$uFinal;
print( i );
}
stopCluster( cluster );
u;
}
If anyone could help that would be much appreciated! But big thanks for
taking your time to read the post!
TelM8
--
View this message in context:
http://r.789695.n4.nabble.com/Problem-with-clusterCall-Error-in-checkForRemoteErrors-lapply-cl-recvResult-tp2338375p2338375.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] How to remove all objects except a few specified objects?
Thanks to De-Jian and Peter. Peter's way is neat and cool! -- View this message in context: http://r.789695.n4.nabble.com/How-to-remove-all-objects-except-a-few-specified-objects-tp2335651p2338221.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] how to plot y-axis on the right of x-axis
Dear Elaine
I'm developing a code to make a 3 y-axes plot. It may can help you.
Also your question leads to mine, posted yesterday: How does R
calculate the interval between tick-marks.
Below follows the code I'm developing.
Data:
ANO CKG NUP NDE
200526352158150025014
200632514789210035024
200727458126180030254
200828568971150043254
200932564789300060320
Code:
# I use comma as decimal symbol so I use dec=","
dat.caupde <- read.delim("clipboard",dec=",",header=T)
attach(dat.caupde)
summary(dat.caupde)
ANOCKGNUPNDE
Min. :2005 Min. :26352158 Min. :1500 Min. :25014
1st Qu.:2006 1st Qu.:27458126 1st Qu.:1500 1st Qu.:30254
Median :2007 Median :28568971 Median :1800 Median :35024
Mean :2007 Mean :29491767 Mean :1980 Mean :38773
3rd Qu.:2008 3rd Qu.:32514789 3rd Qu.:2100 3rd Qu.:43254
Max. :2009 Max. :32564789 Max. :3000 Max. :60320
# here I indicate the limits of each axes and calculate their ranges
y1min <- 26
y1max <- 33
y1range <- y1max-y1min
y2min <- 1500
y2max <- 3000
y2range <- y2max-y2min
y3min <- 25
y3max <- 61
y3range <- y3max-y3min
# making the plot
# y1range*((NUP-y2min)/y2range)+y1min calculates the proportion to
between axes
par(mar=c(6, 6, 2,12))
plot(CKG/100, ylim=c(y1min,y1max), ylab="landings (1000 t)",
type="l",col="red",las=0, cex.axis=1.2,cex.lab=1.4,xaxt="n",xlab="")
points( y1range*((NUP-y2min)/y2range)+y1min,type="l",col="blue")
points( y1range*((NDE/1000-y3min)/y3range)+y1min,type="l",col="darkgreen")
# the number 1 in axis(4,at=seq(y1min,y1max,1 ... is the interval
between each tick-mark
axis(1,at=1:5,labels=as.character(ANO),las=2,cex.axis=1.2)
axis(4,at=seq(y1min,y1max,1),labels=as.integer(seq(y2min,y2max,y2range/y1range)),las=0,cex.axis=1.2)
axis(4,at=seq(y1min,y1max,1),
labels=as.integer(seq(y3min,y3max,y3range/y1range)),cex.axis=1.2,
las=0,line=5)
mtext("nº of fishing boats",4,3,cex=1.4)
mtext("nº of fishing trips (x1000)",4,8, cex=1.4)
legend(4,33,c("L","Nº FB","Nº
FT"),bty="n",lty=c(1,1,1),col=c("red","blue","darkgreen"),cex=1.4)
Now I want to replace the number "1" by the formula used to calculate
the interval between tick-marks. Why, with the range 26-33, R choose
unitary intervals for the y axis (26, 27, 28 ...)?
All the best,
Antonio Olinto
Citando elaine kuo :
Dear List,
I have a richness data distributing across 20 N to 20 S latitude. (120 E-140
E longitude).
I would like to draw the richness in the north hemisphere and a regression
line in the plot
(x-axis: latitude, y-axis: richness in the north hemisphere).
The above demand is done using plot.
Then, south hemisphere richness and regression are required to be generated
using
the same y-axis above but an x-axis on the left side of the y-axis.
(The higher latitude in the south hemisphere, the left it would move)
Please kindly share how to design the south plot and regression line for
richness.
Also, please advise if any more info is in need.
Elaine
[[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.
Webmail - iBCMG Internet
http://www.ibcmg.com.br
__
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] NewJerseyR Meeting
Mango Solutions are proud to announce the first NewJerseyR Meeting to be held
on the 16th September 2010.
Please see details below:
NewJerseyR
Date: Tuesday 16th September
Time: 6.30pm - 9.30pm
Venue: Renaissance Woodbridge Hotel - Iselin, New Jersey
The agenda has yet to be finalised, I shall send details out as soon as they
have been confirmed.
Please contact: newjers...@mango-solutions.com to be added to the NewJerseyR
mailing list to receive meeting updates and R news.
To register for this event please contact: newjers...@mango-solutions.com
We currently hold 2 very successful R meetings - LondonR and BaselR. NewJerseyR
will be our first R meeting in America. NewJerseyR will be a free informal
event, held quarterly and intended to serve as a platform for all local (and
regional) R users to present and exchange their experiences and ideas around
the usage of R. A typical meeting will consist of 3-4 talks of about 20-25 min
to give plenty of room for sharing your R experiences, discussions and exchange
of ideas.
For presentations from our other R meetings can be found at:
www.londonr.org and www.baselr.org
The NewJerseyR website will be up and running shortly.
For further information, please contact newjers...@mango-solutions.com
Sarah Lewis
Hadley Wickham, Creator of ggplot2 - first time teaching in the UK. 1st - 2nd
November 2010.
To book your seat please go to http://mango-solutions.com/news.html
T: +44 (0)1249 767700 Ext: 200
F: +44 (0)1249 767707
M: +44 (0)7746 224226
www.mango-solutions.com
Unit 2 Greenways Business Park
Bellinger Close
Chippenham
Wilts
SN15 1BN
UK
LEGAL NOTICE
This message is intended for the use o...{{dropped:9}}
__
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] GLM outputs in condensed versus expanded table
Hi I'm having different outputs from GLM when using a "condensed" table V1 V2 V3 Present Absent 0 0 0 3 12 0 0 1 0 0 0 1 0 0 0 0 1 1 1 0 1 0 0 7 20 1 0 1 0 0 1 1 0 3 0 1 1 1 6 0 resp=cbind(Present, Absent) glm(resp~V1+V2+V3+I(V1*V2*V3),family=binomial) >Deviance Residuals: [1] 0 0 0 0 0 0 0 0 etc and also coefficients... And when using the same but "expanded" table V1 V2 V3 condition (1 present 0 abscent) Id1 1 0 0 1 id2 1 1 1 1 ... etc glm(condition~V1+V2+V3+I(V1*V2*V3),family=binomial) >Deviance Residuals: Min 1Q Median 3Q Max -0.7747317 -0.7747317 -0.6680472 0.0001315 1.7941226 and also coefficients are different from above. What could I be doing wrong? -- View this message in context: http://r.789695.n4.nabble.com/GLM-outputs-in-condensed-versus-expanded-table-tp2338407p2338407.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] Plot bar lines like excel
Great .. thanks for the to much help and i too appreciate hel p and explanation Cheers -- View this message in context: http://r.789695.n4.nabble.com/Plot-bar-lines-like-excel-tp2337089p2338158.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] Showing a line function on the ploat area
I have a collection of results. I use R to get the linearization presented. how can I get R to show the equation and R^2 value on the plot area next to the graph? -- View this message in context: http://r.789695.n4.nabble.com/Showing-a-line-function-on-the-ploat-area-tp2338389p2338389.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] Merging two data set in R,
try the following "merge" command. merge(A,B, by = intersect(names(A), names(B)), all.x=FALSE, all.y=FALSE) or merge(A,B, by = "ID", all.x=FALSE, all.y=FALSE) Dannemora On Wed, Aug 25, 2010 at 5:35 AM, Mangalani Peter Makananisa < pmakanan...@sars.gov.za> wrote: > Dear R Gurus, > > > > I am currently working on the two dataset ( A and B), they both have the > same fields:ID , REGION, OFFICE, CSTART, CEND, NCYCLE, STATUS and > CB. > > I want to merge the two data set by ID. The problem I have is that the > in data A, the ID's are unique. However in the data set B, the ID's are > not unique, thus some repeat themselves. > > > > How do I the merge or retrieve the common ones? > > Please advise. > > > > Kind Regards > > > > Peter > > > > South Africa > > +27 12 422 7357 > > +27 82 456 4669 > > > > > > > > > > > > > > > > > Please Note: This email and its contents are subject to our email legal > notice which can be viewed at http://www.sars.gov.za/Email_Disclaimer.pdf > >[[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. > [[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] Secant Method Convergence (Method to replicate Excel XIRR/IRR)
Hi,
I am new to R, and as a first exercise, I decided to try to implement an XIRR
function using the secant method. I did a quick search and saw another posting
that used the Bisection method but wanted to see if it was possible using the
secant method.
I would input a Cash Flow and Date vector as well as an initial guess. I
hardcoded today's initial date so I could do checks in Excel. This code seems
to only converge when my initial guess is very close to the correct IRR.
Maybe I have some basic errors in my coding/logic? Any help would be greatly
appreciated.
The Wikipedia article to secant method and IRR:
http://en.wikipedia.org/wiki/Internal_rate_of_return#Numerical_solution
Thanks!
ANXIRR <- function (cashFlow, cfDate, guess){
cfDate<-as.Date(cfDate,format="%m/%d/%Y")
irrprev <- c(0); irr<- guess
pvPrev<- sum(cashFlow)
pv<-
sum(cashFlow/((1+irr)^(as.numeric(difftime(cfDate,"2010-08-24",units="days"))/360)))
print(pv)
print("Hi")
while (abs(pv) >= 0.001) {
t<-irrprev; irrprev<- irr;
irr<-irr-((irr-t)*pv/(pv-pvPrev));
pvPrev<-pv;
pv<-sum(cashFlow/((1+irr)^(as.numeric(difftime(cfDate,"2010-08-24",units="days"))/365)))
print(irr);print(pv)
}
}
Please consider the environment before printing this e-mail.
[[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] lattice help required
The multiple y axes are protecting you in this situation.
z <- cbind(rnorm(100,c(1,10),1), rnorm(100,c(20,30),1))
dotplot(z[,1]+z[,2] ~ facs$Treatment|facs$Sites,
outer=TRUE,
scales = list(
y = list(
relation="free")),
ylab=c("y1", "y2"),
xlab=c("Site 1", "Site 2"),
strip=FALSE,
main="problem")
dotplot(z[,1]+z[,2] ~ facs$Treatment|facs$Sites,
outer=TRUE,
scales = list(
y = list(
relation="free",
limits=list(c(-5,13),c(-5,13),c(18,32),c(18,32,
ylab=c("y1", "y2"),
xlab=c("Site 1", "Site 2"),
strip=FALSE, main="protected")
For more control (such as suppressing the y-tick labels in the right-hand
column,
I recommend Deepayan Sarkar's book.
Rich
[[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] lattice help required
... i added relation="free" to account for diffferent ranges of y1 and y2:
dotplot(y1+y2 ~ facs$Treatment|facs$Sites,
outer=TRUE,
scales = list(y = list(relation="free"), x = list(rot = 90,
tck=c(1,0))),
ylab=c("y1", "y2"),
xlab=c("Site 1", "Site 2"),
strip=FALSE)
but then i get four different y-lims.
more suited there should be only two - one for each response.
can someone tell me how i have to change the code?
best,
kay
> Kay,
>
> doe this do what you want?
>
>
> dotplot(y1+y2 ~ facs$Treatment|facs$Sites,
> outer=TRUE,
> scales = list(x = list(rot = 90, tck=c(1,0))),
> ylab=c("y1", "y2"),
> xlab=c("Site 1", "Site 2"),
> strip=FALSE)
>
>
> On Wed, Aug 25, 2010 at 11:04 AM, Kay Cichini
> wrote:
>
>>
>> hello,
>>
>> i want to stack two lattice plots beneath each other using one x-axis and
>> sharing the same text-panels,
>> like:
>> #
>> library(lattice)
>>
>> y1 <- rnorm(100,100,10)
>> y2 <- rnorm(100,10,1)
>> facs<-expand.grid(Sites=rep(c("Site I","Site
>> II"),25),Treatment=c("A","B"))
>> pl1<-dotplot(y1 ~ facs$Treatment|facs$Sites,
>> scales = list(x = list(rot = 90, tck=c(1,0
>> pl2<-dotplot(y2 ~ facs$Treatment|facs$Sites,
>> scales = list(x = list(rot = 90, tck=c(1,0
>>
>> print(pl1, split=c(1,2,1,2), more=TRUE)
>> print(pl2, split=c(1,1,1,2))
>> #
>>
>> but as said, ideally the plots should be stacked with only the lower plot
>> giving the x-axis annotation
>> and only the upper plot with text-panels.
>>
>> thanks a lot,
>> kay
>>
>> -
>>
>> Kay Cichini
>> Postgraduate student
>> Institute of Botany
>> Univ. of Innsbruck
>>
>>
>> --
>> View this message in context:
>> http://r.789695.n4.nabble.com/lattice-help-required-tp2338382p2338382.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.
-
Kay Cichini
Postgraduate student
Institute of Botany
Univ. of Innsbruck
--
View this message in context:
http://r.789695.n4.nabble.com/lattice-help-required-tp2338382p2338641.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] nls self starting function
I need the simple function for the following set of data. This is a toy version of my data, but the error is persistent in both. To compare with excel, I would just 1) format trendline, 2) display equation and R-squared on chart. I obviously tried to use a nls and wanted to use the self startup function SSasymp rm(list=ls()) t <- c( 0,.1,.2,.7,.4,.5,.6,1,.8,.4,1.2,1.3,1.4,1.5,1.6, 1.6,1.7,1.8,1.9,2.0, 3,6,11) y<-c(2,3.5,3.01,3.09,3,3.25,3.36,3.49,3.64, 3.81, 4.44, 4.69, 4.96, 5.25, 6, 6.1, 5.89, 6.24, 6.61, 7.00, 12.00, 39, 124.00) plot(y~t ,xlab="t", ylab="y") model<-nls(y~SSasymp(t,a,b)) summary(model) Thanks in advance, keith -- M. Keith Cox, Ph.D. Alaska NOAA Fisheries, National Marine Fisheries Service Auke Bay Laboratories 17109 Pt. Lena Loop Rd. Juneau, AK 99801 keith@noaa.gov marlink...@gmail.com U.S. (907) 789-6603 [[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 obtain seed after generating random number?
If you find yourself doing things like this often, but don't want to explicitly
set the seed, or save seeds before simulating, then you can run the following
code (or put it into .Rprofile or similar):
.Last.Random.seed <- .Random.seed
addTaskCallback( function(expr, val, ok, visible){
if(!isTRUE( all.equal(.Last.Random.seed, .Random.seed)) ) {
.Last.Random.seed <- .Random.seed
}
TRUE
})
Then the previous seed will be stored in .Last.Random.seed and you can restore
the seed to be able to rerun the same values, e.g.:
> rnorm(10)
[1] -0.28361138 0.86951931 -0.54435528 0.62880324 -1.42233446 -1.22751263
[7] -1.67410552 0.08439848 -0.20612566 1.44187164
> .Random.seed <- .Last.Random.seed
> rnorm(10)
[1] -0.28361138 0.86951931 -0.54435528 0.62880324 -1.42233446 -1.22751263
[7] -1.67410552 0.08439848 -0.20612566 1.44187164
> rnorm(10)
[1] -0.0417821 1.3537545 1.9452253 -0.4909382 0.3884391 -0.8448933
[7] 0.7379904 -1.0797603 -1.0264739 0.2887934
The above code only keeps the most recent seed, but the code could be modified
to store a longer history if that were desired.
If you want to set your own seeds (a little easier to save, pass to others),
but find integers to unimaginative, then look at the char2seed function in the
TeachingDemos package.
Hope this helps,
--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111
> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> project.org] On Behalf Of Bogaso Christofer
> Sent: Tuesday, August 24, 2010 11:12 AM
> To: r-help@r-project.org
> Subject: [R] How to obtain seed after generating random number?
>
> Dear all, I was doing an experiment to disprove some theory therefore
> performing lot of random simulation. Goal is to show the audience that
> although something has very rare chance to occur but it doesn't mean
> that
> event would be impossible.
>
>
>
> In this case after getting that rare event I need to show that same
> scenario
> for multiple times to explain other audience. Hence I need to somehow
> save
> that seed which generates that random numbers after doing the
> experiment.
> However as it is very rare event it is not very practical to start with
> a
> fixed seed and then generate random numbers. Hence I am looking for
> some way
> which will tell me about that corresponding seed which was responsible
> to
> generate that particular series of random numbers responsible for
> occurrence
> of that rare event.
>
>
>
> In short, I need to know the seed ***after*** generating the random
> numbers.
>
>
>
> Is there any possibility to know this?
>
>
>
> Thanks and regards,
>
>
> [[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] Change value of a slot of an S4 object within a method.
Hi Steve,
thanks for the tip. I'll definitely take a closer look at your
solution for implementation for future use. But right now I don't
have the time to start rewriting my class definitions.
Luckily, I found where exactly things were going wrong. After reading
into the documentation about the evaluation in R, I figured out I have
to specify the environment where substitute should look explicitly as
parent.frame(1). I still don't understand completely why exactly, but
it does the job.
Thus :
eval(eval(substitute(expression(obj...@extra[[name]] <<- value
should become :
eval(
eval(
substitute(
expression(obj...@extra[[name]] <<- value)
,env=parent.frame(1) )
)
)
Tried it out and it works.
Cheers
Joris
On Wed, Aug 25, 2010 at 6:21 PM, Steve Lianoglou
wrote:
> Hi Joris,
>
> On Wed, Aug 25, 2010 at 11:56 AM, Joris Meys wrote:
>> Dear all,
>>
>> I have an S4 class with a slot "extra" which is a list. I want to be
>> able to add an element called "name" to that list, containing the
>> object "value" (which can be a vector, a dataframe or another S4
>> object)
>>
>> Obviously
>>
>> setMethod("add.extra",signature=c("PM10Data","character","vector"),
>> function(object,name,value){
>> obj...@extra[[name]] <- value
>> }
>> )
>>
>> Contrary to what I would expect, the line :
>> eval(eval(substitute(expression(obj...@extra[[name]] <<- value
>>
>> gives the error :
>> Error in obj...@extra[[name]] <<- value : object 'object' not found
>>
>> Substitute apparently doesn't work any more in S4 methods...
>>
>> I found a work-around by calling the initializer every time, but this
>> influences the performance. Returning the object is also not an
>> option, as I'd have to remember to assign that one each time to the
>> original name.
>>
>> Basically I'm trying to do some call by reference with S4, but don't
>> see how I should do that. How would I be able to tackle this problem
>> in an efficient and elegant way?
>
> In lots of my own S4 classes I define a slot called ".cache" which is
> an environment for this exact purpose.
>
> Using this solution for your scenario might look something like this:
>
> setMethod("add.extra",signature=c("PM10Data","character","vector"),
> function(object, name, value) {
> obj...@.cache$extra[[name]] <- value
> })
>
> I'm not sure what your entire problem looks like, but to "get" your
> extra list, or a value form it, you could:
>
> setMethod("extra", signature="PM10Data",
> function(object, name=NULL) {
> if (!is.null(name)) {
> obj...@.cache$extra[[name]]
> } else {
> obj...@.cache$extra
> })
>
> ... or something like that.
>
> The last thing you have to be careful of is that you nave to make sure
> that each new("PM10Data") object you have initializes its *own* cache:
>
> setClass("PM10Data", representation(..., .cache='environment'))
> setMethod("initialize", "PM10Data",
> function(.Object, ..., .cache=new.env()) {
> callNextMethod(.Object, .cache=.cache, ...)
> })
>
> Does that make sense?
>
> --
> Steve Lianoglou
> Graduate Student: Computational Systems Biology
> | Memorial Sloan-Kettering Cancer Center
> | Weill Medical College of Cornell University
> Contact Info: http://cbio.mskcc.org/~lianos/contact
>
__
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] frequency, count rows, data for heat map
Hi:
Here are a couple of ways to render a basic 2D table. Let's call your input
data frame dat:
> names(dat) <- c('samp', 'sequen')
> ssTab <- as.data.frame(with(dat, table(samp, sequen)))
> ssTab # data frame version
samp sequen Freq
1 111abc1
2 1079abc1
3 5576abc1
4 111sdf1
5 1079sdf0
6 5576sdf2
7 111xyz1
8 1079xyz2
9 5576xyz0
> with(dat, table(samp, sequen)) # table version
sequen
samp abc sdf xyz
1111 1 1
1079 1 0 2
5576 1 2 0
HTH,
Dennis
On Wed, Aug 25, 2010 at 7:53 AM, rtsweeney wrote:
>
> Hi all,
> I have read posts of heat map creation but I am one step prior --
> Here is what I am trying to do and wonder if you have any tips?
> We are trying to map sequence reads from tumors to viral genomes.
>
> Example input file :
> 111 abc
> 111 sdf
> 111 xyz
> 1079 abc
> 1079 xyz
> 1079 xyz
> 5576 abc
> 5576 sdf
> 5576 sdf
>
> How may xyz's are there for 1079 and 111? How many abc's, etc?
> How many times did reads from sample (1079) align to virus xyz.
> In some cases there are thousands per virus in a give sample, sometimes
> one.
> The original file (two columns by tens of thousands of rows; 20 MB) is
> text file (tab delimited).
>
> Output file:
> abc sdf xyz
> 111 1 1 1
> 1079 1 0 2
> 5576 1 2 0
>
> Or, other ways to generate this data so I can then use it for heat map
> creation?
>
> Thanks for any help you may have,
>
> rtsweeney
> palo alto, ca
> --
> View this message in context:
> http://r.789695.n4.nabble.com/frequency-count-rows-data-for-heat-map-tp2338363p2338363.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.
[R] several odfWeave questions
[Sending both to the maintainer and to R-help, in case anyone else has
answers ...]
I've looked in odfWeave documentation, vignette, and poked around on
the web some, and haven't found answers yet.
1a. am I right in believing that odfWeave does not respect the
'keep.source' option? Am I missing something obvious?
1b. is there a way to set global options analogous to \SweaveOpts{}
directives in Sweave? (I looked at odfWeaveControl, it doesn't seem to
do it.)
2. I tried to write a Makefile directive to process files from the
command line:
%.odt: %_in.odt
$(RSCRIPT) -e "library(odfWeave);
odfWeave(\"$*_in.odt\",\"$*.odt\");"
This works, *but* the resulting output file gives a warning ("The file
'odftest2.odt' is corrupt and therefore cannot be opened.
OpenOffice.org can try to repair the file ..."). Based on looking at
the contents, it seems that a spurious/unnecessary 'Rplots.pdf' file is getting
created and zipped in with the rest of the archive; when I unzip, delete
the Rplots.pdf file and re-zip, the ODT file opens without a warning.
Obviously I could post-process but it would be nicer to find a
workaround within R ...
3. I find the requirement that all file paths be specified as absolute
rather than relative paths somewhat annoying -- I understand the reason,
but it goes against one practice that I try to encourage for
reproducibility, which is *not* to use absolute file paths -- when
moving a same set of data and analysis files across computers, it's hard
to enforce them all ending up in the same absolute location, which then
means that the recipient has to edit the ODT file. It would be nice if
there were hooks for read.table() and load() as there are for plotting
and package/namespace loading -- then one could just copy them into the
working directory on the fly.
has anyone experienced this/thought of any workarounds?
(I guess one solution is to zip any necessary source files into the archive
beforehand,
as illustrated in the vignette.)
My test files are posted at
http://www.math.mcmaster.ca/~bolker/junk/odftest_in.odt and
http://www.math.mcmaster.ca/~bolker/junk/odftest.odt
thanks,
Ben Bolker
> > sessionInfo()
>
R version 2.11.1 (2010-05-31)
i486-pc-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] odfWeave_0.7.16 XML_3.1-1 lattice_0.19-9
loaded via a namespace (and not attached):
[1] grid_2.11.1
__
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] Repeat the first day data through all the day. Zoo
On Wed, Aug 25, 2010 at 11:18 AM, skan wrote:
>
> # duplicated / na.locf doesn't work
> it says Error in fix.by(by.x, x) : 'by' must specify valid column(s)
>
> if I use ifelse instead of ifelse.zoo it works but it gives me a non zoo
> vector.
> Myabe is because my zoo version is older.
>
They all work:
> library(zoo)
> library(chron)
> z <- zoo(1:10, chron(0:9/5))
>
> # aggregate / na.locf
> z.ag <- aggregate(z, as.Date, head, 1)
> na.locf(z.ag, xout = time(z))
(01/01/70 00:00:00) (01/01/70 04:48:00) (01/01/70 09:36:00) (01/01/70 14:24:00)
1 1 1 1
(01/01/70 19:12:00) (01/02/70 00:00:00) (01/02/70 04:48:00) (01/02/70 09:36:00)
1 6 6 6
(01/02/70 14:24:00) (01/02/70 19:12:00)
6 6
>
> # duplicated / na.locf
> z.na <- ifelse.zoo(!duplicated(as.Date(time(z))), z, NA)
> na.locf(z.na)
(01/01/70 00:00:00) (01/01/70 04:48:00) (01/01/70 09:36:00) (01/01/70 14:24:00)
1 1 1 1
(01/01/70 19:12:00) (01/02/70 00:00:00) (01/02/70 04:48:00) (01/02/70 09:36:00)
1 6 6 6
(01/02/70 14:24:00) (01/02/70 19:12:00)
6 6
>
> # ave - as before
> zz <- z
> zz[] <- ave(coredata(z), as.Date(time(z)), FUN = function(x) head(x, 1))
> zz
(01/01/70 00:00:00) (01/01/70 04:48:00) (01/01/70 09:36:00) (01/01/70 14:24:00)
1 1 1 1
(01/01/70 19:12:00) (01/02/70 00:00:00) (01/02/70 04:48:00) (01/02/70 09:36:00)
1 6 6 6
(01/02/70 14:24:00) (01/02/70 19:12:00)
6 6
> packageDescription("zoo")$Version
[1] "1.6-4"
> R.version.string
[1] "R version 2.11.1 Patched (2010-05-31 r52167)"
__
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] lattice help required
exactly -
thanks a lot, richard!
kay
Zitat von "RICHARD M. HEIBERGER" :
Kay,
doe this do what you want?
dotplot(y1+y2 ~ facs$Treatment|facs$Sites,
outer=TRUE,
scales = list(x = list(rot = 90, tck=c(1,0))),
ylab=c("y1", "y2"),
xlab=c("Site 1", "Site 2"),
strip=FALSE)
On Wed, Aug 25, 2010 at 11:04 AM, Kay Cichini wrote:
hello,
i want to stack two lattice plots beneath each other using one x-axis and
sharing the same text-panels,
like:
#
library(lattice)
y1 <- rnorm(100,100,10)
y2 <- rnorm(100,10,1)
facs<-expand.grid(Sites=rep(c("Site I","Site II"),25),Treatment=c("A","B"))
pl1<-dotplot(y1 ~ facs$Treatment|facs$Sites,
scales = list(x = list(rot = 90, tck=c(1,0
pl2<-dotplot(y2 ~ facs$Treatment|facs$Sites,
scales = list(x = list(rot = 90, tck=c(1,0
print(pl1, split=c(1,2,1,2), more=TRUE)
print(pl2, split=c(1,1,1,2))
#
but as said, ideally the plots should be stacked with only the lower plot
giving the x-axis annotation
and only the upper plot with text-panels.
thanks a lot,
kay
-
Kay Cichini
Postgraduate student
Institute of Botany
Univ. of Innsbruck
--
View this message in context:
http://r.789695.n4.nabble.com/lattice-help-required-tp2338382p2338382.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.
[R] package MuMI
Hello, I am using the package "MuMI" to run all the possible combinations of variables in my full model, and select my best models. When I enter my variables in the original model I write them like this lm(y~ a +b +c +a:b) However, "MuMI" will also use the variable b:a, which I do not want in my model. How do I stop that from happening? Thank you, Agnese __ 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] Odp: Finding pairs
Hi: I'm just ideating here (think IBM commercial...) but perhaps a graphical model approach might be worth looking into. It seems to me that Mr. Rhodes is looking for clusters of banks that are under the same ownership umbrella. That information is not directly available in a single variable, but can evidently be inferred from the matches between the two variables: B[i] controls A[i] if B[i] is nonempty. In the bank 144 -> 147 -> 149 example, 149 controls 147 and 147 controls 144, so it appears that some transitive relation holds among the set of matches as well. (Why is PacMan going through my head? :) I know next to nothing about graphical models, but I'm thinking about igraph and some of the tools in the statnet bundle to tackle this problem. Does that make sense to anyone? Alternatives? FWIW, Dennis On Wed, Aug 25, 2010 at 2:24 AM, Mike Rhodes wrote: > Dear Mr Petr PIKAL > After reading the R code provided by you, I realized that I would have > never figured out how this could have been done. I am going to re-read again > and again your code to understand the logic and the commands you have > provided. > Thanks again from the heart for your kind advice. > Regards > Mike > > --- On Wed, 25/8/10, Petr PIKAL wrote: > > From: Petr PIKAL > Subject: Re: [R] Odp: Finding pairs > To: "Mike Rhodes" > Cc: r-help@r-project.org > Date: Wednesday, 25 August, 2010, 9:01 > > Hm > > r-help-boun...@r-project.org napsal dne 25.08.2010 09:43:26: > > > Dear Mr Petr Pikal > > > > I am extremely sorry for the manner I have raised the query. Actually > that was > > my first post to this R forum and in fact even I was also bit confused > while > > drafting the query, for which I really owe sorry to all for consuming > the > > precious time. Perhaps I will try to redraft my query in a better way as > follows. > > > > I have two datasets "A" and "B" containing the names of branch offices > of a > > particular bank say XYZ plc bank. The XYZ bank has number of main branch > > > offices (say Parent) and some small branch offices falling under the > purview > > of these main branch offices (say Child). > > > > The datalist "A" and "B" consists of these main branch office names as > well as > > small branch office names. B is subset of A and these branch names are > coded. > > Thus we have two datasets A and B as (again I am using only a > > portion of a large database just to have some idea) > > > > > > A B > > 144 >what is here in B? Empty space?, > > 145 > > 146 > > 147 144 > > How do you know that 144 from B relates to 147 in A? Is it according to > its positions? I.e. 4th item in B belongs to 4.th item in A? > > > 148 145 > > > > 149 147 > > 151 148 > > > > > > > > Now the branch 144 appears in A as well as in B and in B it is mapped > with > > 147. This means branch 147 comes under the purview of main branch 144. > Again > > 147 is controlling the branch 149 (since 147 also has appeared in B and > is > > mapped with 149 of A). > > > > Similarly, branch 145 is controlling branch 148 which further controls > > operations of bank branch 151 and like wise. > > Well as you did not say anything about structure of your data > A<-144:151 > B<-144:148 > data.frame(A,B) > A B > 1 144 NA > 2 145 NA > 3 146 NA > 4 147 144 > 5 148 145 > 6 149 146 > 7 150 147 > 8 151 148 > DF<-data.frame(A,B) > main<-DF$A[is.na(DF$B)] > branch1<-DF[!is.na(DF$B),] > selected.branch1<-branch1$A[branch1$B%in%main] > branch2<-branch1[!branch1$B%in%main,] > selected.branch2<-branch2$A[branch2$B%in%selected.branch1] > > and for cbinding your data which has uneven number of values see Jim > Holtman's answer to this > > How to cbind DF:s with differing number of rows? > > Regards > Petr > > > > > > So in the end I need an output something like - > > > > Main Branch Branch office1 Branch > > office2 > > 144 147 149 > > 145 148 151 > > > 146 NA > > NA > > > > ... > > > > .. > > > > > > I understand again I am not able to put forward my query properly. But I > must > > thank all of you for giving a patient reading to my query and for > reverting > > back earlier. Thanks once again. > > > > With warmest regards > > > > Mike > > > > > > --- On Wed, 25/8/10, Petr PIKAL wrote: > > > > From: Petr PIKAL > > Subject: Odp: [R] Finding > > pairs > > To: "Mike Rhodes" > > Cc: r-help@r-project.org > > Date: Wednesday, 25 August, 2010, 6:39 > > > > Hi > > > > without other details it is probably impossible to give you any > reasonable > > advice. Do you have your data already in R? Wha
Re: [R] Change value of a slot of an S4 object within a method.
Hi Joris,
On Wed, Aug 25, 2010 at 11:56 AM, Joris Meys wrote:
> Dear all,
>
> I have an S4 class with a slot "extra" which is a list. I want to be
> able to add an element called "name" to that list, containing the
> object "value" (which can be a vector, a dataframe or another S4
> object)
>
> Obviously
>
> setMethod("add.extra",signature=c("PM10Data","character","vector"),
> function(object,name,value){
> obj...@extra[[name]] <- value
> }
> )
>
> Contrary to what I would expect, the line :
> eval(eval(substitute(expression(obj...@extra[[name]] <<- value
>
> gives the error :
> Error in obj...@extra[[name]] <<- value : object 'object' not found
>
> Substitute apparently doesn't work any more in S4 methods...
>
> I found a work-around by calling the initializer every time, but this
> influences the performance. Returning the object is also not an
> option, as I'd have to remember to assign that one each time to the
> original name.
>
> Basically I'm trying to do some call by reference with S4, but don't
> see how I should do that. How would I be able to tackle this problem
> in an efficient and elegant way?
In lots of my own S4 classes I define a slot called ".cache" which is
an environment for this exact purpose.
Using this solution for your scenario might look something like this:
setMethod("add.extra",signature=c("PM10Data","character","vector"),
function(object, name, value) {
obj...@.cache$extra[[name]] <- value
})
I'm not sure what your entire problem looks like, but to "get" your
extra list, or a value form it, you could:
setMethod("extra", signature="PM10Data",
function(object, name=NULL) {
if (!is.null(name)) {
obj...@.cache$extra[[name]]
} else {
obj...@.cache$extra
})
... or something like that.
The last thing you have to be careful of is that you nave to make sure
that each new("PM10Data") object you have initializes its *own* cache:
setClass("PM10Data", representation(..., .cache='environment'))
setMethod("initialize", "PM10Data",
function(.Object, ..., .cache=new.env()) {
callNextMethod(.Object, .cache=.cache, ...)
})
Does that make sense?
--
Steve Lianoglou
Graduate Student: Computational Systems Biology
| Memorial Sloan-Kettering Cancer Center
| Weill Medical College of Cornell University
Contact Info: http://cbio.mskcc.org/~lianos/contact
__
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] Compiling Fortran for R : .so: mach-o, but wrong architecture
Hi Marie, this link may be helpful if you want to build it for both architectures, 32-bit and 64-bit. I struggled with the same problem not too long ago. http://cran.rakanu.com/bin/macosx/RMacOSX-FAQ.html#Building-universal-package I ended up writing a very simple Makefile to accompany [my version of] the bar.f code. There is almost certainly a better way to do it. Maybe with Makevars? David -- View this message in context: http://r.789695.n4.nabble.com/Compiling-Fortran-for-R-so-mach-o-but-wrong-architecture-tp2337198p2338507.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] multiple x factors using the sciplot package
Hi: You can probably do what you want in either ggplot2 or lattice, but I would recommend at least a couple different approaches: (1) Plot individual bar charts by combinations of year and period. This is easy to do in both ggplot2 and lattice: in ggplot2, one would use geom_bar(x) + facet_grid(year ~ period), where x is the count variable. The same thing in lattice would be barchart( ~ x | year * period, data = mydata). Suitable options can be added to either plot for better presentation. All the plots are on the same horizontal and vertical scales so that visual comparisons across the rows and columns are easy to make. (2) Use a Cleveland dot chart. This is much less wasteful of ink than bar charts, and multiple groups can be compared without too much difficulty. See dotplot() in lattice for details. The VADeaths example is likely the best one to emulate for your problem. Both the input data structure and plot code matter. (3) If your purpose is to compare the two groups with respect to a count variable, a third option is to consider a mosaic plot, which can be thought of as a 'visual chi-square test'. These can be found in package vcd, which contains a nice 40+ page vignette to show you in detail how to create a suitable mosaic plot (and how to interpret it). I prefer to see R used to promote good graphics and sound statistical principles. In that spirit, performing a two-factor comparison with count data is easier, and quite possibly more instructive, with a 2D visual metaphor (which all three of the suggestions above have in common) than with the 1D metaphor of multiple dodged/side-by-side bar charts. Of course, this may be against the standard practice in your field, in which case you have a decision to make. HTH, Dennis On Wed, Aug 25, 2010 at 6:16 AM, Jordan Ouellette-Plante < jordan_opla...@hotmail.com> wrote: > > Dear R community, > > > > I am a beginner using the sciplot package to graph barplots. I would like > to be able to graph two x factors (Sampling.year and Period). > > Sampling.year: 2006, 2007, 2008, 2009 > > Period: First, Second, Total > > > > The parameter "group" is the different species I looked at. They can be > seen in the legend. > > The parameter "response" is the Average percentage cover category of these > species. > > > > I would like the graph so you see on the x axis the years (in a bigger > font) and the period (in smaller font). It would look a bit like the barplot > that can be seen at > http://onertipaday.blogspot.com/2007/05/make-many-barplot-into-one-plot.html, > but with the advantage of the sciplot package. > > > > > > Thanks a lot for your help! > > > > Jordan > >[[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. > [[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] Change value of a slot of an S4 object within a method.
Dear all,
I have an S4 class with a slot "extra" which is a list. I want to be
able to add an element called "name" to that list, containing the
object "value" (which can be a vector, a dataframe or another S4
object)
Obviously
setMethod("add.extra",signature=c("PM10Data","character","vector"),
function(object,name,value){
obj...@extra[[name]] <- value
}
)
Contrary to what I would expect, the line :
eval(eval(substitute(expression(obj...@extra[[name]] <<- value
gives the error :
Error in obj...@extra[[name]] <<- value : object 'object' not found
Substitute apparently doesn't work any more in S4 methods...
I found a work-around by calling the initializer every time, but this
influences the performance. Returning the object is also not an
option, as I'd have to remember to assign that one each time to the
original name.
Basically I'm trying to do some call by reference with S4, but don't
see how I should do that. How would I be able to tackle this problem
in an efficient and elegant way?
Thank you in advance
Cheers
Joris
--
Joris Meys
Statistical consultant
Ghent University
Faculty of Bioscience Engineering
Department of Applied mathematics, biometrics and process control
tel : +32 9 264 59 87
joris.m...@ugent.be
---
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php
__
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] Draw a perpendicular line?
At 3:04 PM -0700 8/23/10, CZ wrote:
Hi,
I am trying to draw a perpendicular line from a point to two points.
Mathematically I know how to do it, but to program it, I encounter some
problem and hope can get help. Thanks.
I have points, A, B and C. I calculate the slope and intercept for line
drawn between A and B.
I am trying to check whether I can draw a perpendicular line from C to line
AB and get the x,y value for the point D at the intersection.
Assume I get the slope of the perpendicular line, I will have my point (D)
using variable x and y which is potentially on line AB. My idea was using
|AC|*|AC| = |AD|*|AD|+ |CD|*|CD|. I don't know what function I may need to
call to calculate the values for point D (uniroot?).
This is easier than you think.
Think of the x,y coordinates of each point :
Then, the slope is slope = rise/run = (By- Ay)/(Bx- Ax)
The Dx coordinate = Cx and the Dy = (Dx - Ax) * slope
Then, to draw the line segment from C to D
lines(C,D)
In R:
A <- c(2,4)
B <- c(4,1)
C <- c(8,10)
slope <-( C[2]- A[2])/(C[1]-A[1]) #rise/run
D <- c(B[1],(B[1]-A[1])*slope + A[2]) # find D
my.data <- rbind(A,B,C,D)
colnames(my.data) <- c("X","Y")
my.data#show it
plot(my.data,type="n") #graph it without the points
text(my.data,rownames(my.data)) #put the points in
segments(A[1],A[2],C[1],C[2]) #draw the line segments
segments(B[1],B[2],D[1],D[2])
Bill
Thank you.
--
View this message in context:
http://r.789695.n4.nabble.com/Draw-a-perpendicular-line-tp2335882p2335882.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] lattice help required
Kay,
doe this do what you want?
dotplot(y1+y2 ~ facs$Treatment|facs$Sites,
outer=TRUE,
scales = list(x = list(rot = 90, tck=c(1,0))),
ylab=c("y1", "y2"),
xlab=c("Site 1", "Site 2"),
strip=FALSE)
On Wed, Aug 25, 2010 at 11:04 AM, Kay Cichini wrote:
>
> hello,
>
> i want to stack two lattice plots beneath each other using one x-axis and
> sharing the same text-panels,
> like:
> #
> library(lattice)
>
> y1 <- rnorm(100,100,10)
> y2 <- rnorm(100,10,1)
> facs<-expand.grid(Sites=rep(c("Site I","Site II"),25),Treatment=c("A","B"))
> pl1<-dotplot(y1 ~ facs$Treatment|facs$Sites,
> scales = list(x = list(rot = 90, tck=c(1,0
> pl2<-dotplot(y2 ~ facs$Treatment|facs$Sites,
> scales = list(x = list(rot = 90, tck=c(1,0
>
> print(pl1, split=c(1,2,1,2), more=TRUE)
> print(pl2, split=c(1,1,1,2))
> #
>
> but as said, ideally the plots should be stacked with only the lower plot
> giving the x-axis annotation
> and only the upper plot with text-panels.
>
> thanks a lot,
> kay
>
> -
>
> Kay Cichini
> Postgraduate student
> Institute of Botany
> Univ. of Innsbruck
>
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/lattice-help-required-tp2338382p2338382.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.
Re: [R] Repeat the first day data through all the day. Zoo
# duplicated / na.locf doesn't work it says Error in fix.by(by.x, x) : 'by' must specify valid column(s) if I use ifelse instead of ifelse.zoo it works but it gives me a non zoo vector. Myabe is because my zoo version is older. cheers -- View this message in context: http://r.789695.n4.nabble.com/Repeat-the-first-day-data-through-all-the-day-Zoo-tp2338069p2338409.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] Estimate average standard deviation of mean of two dependent groups
Hi Jokel,
If I remember correctly
1) The variance of the sum of two variables X and Y is:
var(X) + var(Y) + (2 * cov(X, Y))
If you create some sample data, you can verify this by:
var((X + Y))
2) The variance of a random variable (X) multiplied by some constant (K)
is equal to the the variance of X multiplied by the square of K. So:
var(X * K) = var(X) * K^2
I have never seen these combined, but I imagine they could be.
Let Z = X + Y
The mean of X and Y is (X + Y)/2 = Z/2
>From the formula above (1), the variance of Z is:
var(X) + var(Y) + (2 * cov(X, Y))
Because Z/2 = Z * 0.5, the variance of Z * 0.5 is given by:
var(Z) * (0.5^2)
and substituting, the variance of the mean of X and Y is:
(var(X) + var(Y) + (2 * cov(X, Y))) * (0.5^2)
This held up under several empirical tests I did where I had actual
data for X and Y.
***Discalimer
I am piecing this together, and I am far from an expert on the
subject. I do not know if it is actually true, and there may be
additional assumptions.
I tested this using:
myfun <- function() {
x <- rnorm(100)
y <- rnorm(100)
var.x <- var(x)
var.y <- var(y)
cov.xy <- cov(x, y)
calc.var.xplusy <- var.x + var.y + 2*cov.xy
var.meanxy <- calc.var.xplusy * (.5^2)
empirical.var.meanxy <- var( (x + y)/2 )
output <- all.equal(var.meanxy, empirical.var.meanxy)
return(output)
}
temp <- vector("logical", 1000)
for(i in 1:1000) {temp[i] <- myfun()}
all(temp)
HTH,
Josh
On Wed, Aug 25, 2010 at 5:29 AM, Jokel Meyer wrote:
> Dear R-experts!
>
> I am currently running a meta-analysis with the help of the great metafor
> package. However I have some difficulties setting up my raw data to enter it
> to the meta-analysis models.
> I have a group of subjects that have been measured in two continuous
> variables (A & B). I have the information about the mean of the two
> variables, the group size (n) and the standard deviations of the two
> variables.
> Now I would like to average both variables (A & B) and get the mean and
> standard deviation of the merged variable (C).
> As for the mean this would be quiet easy: I would just take the mean of mean
> A and mean B to get the mean of C.
> However for the standard deviation this seems more tricky as it is to assume
> that standard deviations in A & B correlate. I assume (based on further
> analysis) a correlation of r =0.5.
> I found the formula to get the standard deviation of the SUM (not the mean)
> of two variables:
> SD=SQRT(SD_A^2 + SD_B^2 + 2*r*SD_A*SD_B)
>
> with SD_B and SD_B being the standard deviation of A and B. And r*SD_A*SD_B
> being the covariance of A and B.
>
> Would this formula also be valid if I want to average (and not sum) my two
> variables?
>
> Many thanks for any help & best wishes,
> Jokel
>
> [[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] frequency, count rows, data for heat map
Your problem is not completely clear to me, but perhaps something like
data <- data.frame(
a = rep(c(1,2), each=10),
b = rep(c('a', 'b', 'c', 'd'), 5))
library(plyr)
daply(data, a ~ b, nrow)
does what you need.
Regards,
Jan
On Wed, Aug 25, 2010 at 4:53 PM, rtsweeney wrote:
>
> Hi all,
> I have read posts of heat map creation but I am one step prior --
> Here is what I am trying to do and wonder if you have any tips?
> We are trying to map sequence reads from tumors to viral genomes.
>
> Example input file :
> 111 abc
> 111 sdf
> 111 xyz
> 1079 abc
> 1079 xyz
> 1079 xyz
> 5576 abc
> 5576 sdf
> 5576 sdf
>
> How may xyz's are there for 1079 and 111? How many abc's, etc?
> How many times did reads from sample (1079) align to virus xyz.
> In some cases there are thousands per virus in a give sample, sometimes one.
> The original file (two columns by tens of thousands of rows; 20 MB) is
> text file (tab delimited).
>
> Output file:
> abc sdf xyz
> 111 1 1 1
> 1079 1 0 2
> 5576 1 2 0
>
> Or, other ways to generate this data so I can then use it for heat map
> creation?
>
> Thanks for any help you may have,
>
> rtsweeney
> palo alto, ca
> --
> View this message in context:
> http://r.789695.n4.nabble.com/frequency-count-rows-data-for-heat-map-tp2338363p2338363.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] lmer() causes segfault
Dennis, just wow. Thank you so much. I knew it was something trivial - in this case the variable type of the of the grouping variables. However, something as trivial as this should not throw a segfault IMHO. I tried subscribing to R-sig-mixed this morning, but the corresponding mail server at the ETH's stats department seems to be down. And thank you so much for changing the model, that is a great new starting point. Can you recommend a good book that deals with multilevel models in lmer() that include longitudinal data? I was not aware of the difference between scalar random effects and random slopes and would like to read up on that. Again, thanks a lot. Regards, Bertolt Am 25.08.2010 um 13:47 schrieb Dennis Murphy: Hi: Let's start with the data: > str(test.data) 'data.frame': 100 obs. of 4 variables: $ StudentID: num 17370 17370 17370 17370 17379 ... $ GroupID : num 1 1 1 1 1 1 1 1 1 1 ... $ Time : num 1 2 3 4 1 2 3 4 1 2 ... $ Score: num 76.8 81.8 89.8 92.8 75.9 ... Both StudentID and GroupID are numeric; in the model, they would be treated as continuous covariates rather than factors, so we need to convert: test.data$StudentID <- factor(test.data$StudentID) test.data$GroupID <- factor(test.data$GroupID) Secondly, I believe there are some flaws in your model. After converting your variables to factors, I ran library(lme4) mlmoded1.lmer <- lmer(Score ~ Time + (Time | GroupID/StudentID), data = test.data) You have two groups, so they should be treated as a fixed effect - more specifically, as a fixed blocking factor. The StudentIDs are certainly nested within GroupID, and Time is measured on each StudentID, so it is a repeated measures factor. The output of this model is > mlmoded1.lmer Linear mixed model fit by REML Formula: Score ~ Time + (Time | GroupID/StudentID) Data: test.data AIC BIC logLik deviance REMLdev 393.1 416.5 -187.5376.9 375.1 Random effects: GroupsNameVariance Std.Dev. Corr StudentID:GroupID (Intercept) 0.504131 0.71002 Time 0.083406 0.28880 1.000 GroupID (Intercept) 12.809567 3.57905 Time 3.897041 1.97409 -1.000 Residual 1.444532 1.20189 Number of obs: 100, groups: StudentID:GroupID, 25; GroupID, 2 Fixed effects: Estimate Std. Error t value (Intercept) 72.803 2.552 28.530 Time 4.474 1.401 3.193 Correlation of Fixed Effects: (Intr) Time -0.994 The high correlations among the random effects and then among the fixed effects suggests that the model specification may be a bit off. The above model fits random slopes to GroupIDs and StudentIDs, along with random intercepts, but GroupID is a between-subject effect and should be at the top level. Time is a within-subject effect and StudentIDs are the observational units. I modified the model to provide fixed effects for GroupIDs, scalar random effects for StudentIDs and random slopes for StudentIDs. > mod3 <- lmer(Score ~ 1 + GroupID + Time + (1 | StudentID) + + (0 + Time | StudentID), data = test.data) > mod3 Linear mixed model fit by REML Formula: Score ~ 1 + GroupID + Time + (1 | StudentID) + (0 + Time | StudentID) Data: test.data AIC BIC logLik deviance REMLdev 430.9 446.5 -209.4418.4 418.9 Random effects: GroupsNameVariance Std.Dev. StudentID (Intercept) 4.2186e-13 6.4951e-07 StudentID Time1.8380e+00 1.3557e+00 Residual 1.6301e+00 1.2768e+00 Number of obs: 100, groups: StudentID, 25 Fixed effects: Estimate Std. Error t value (Intercept) 70.7705 0.4204 168.33 GroupID2 4.0248 0.58546.88 Time 4.5292 0.2942 15.39 Correlation of Fixed Effects: (Intr) GrpID2 GroupID2 -0.668 Time -0.264 0.000 I didn't check the quality of the fit, but on the surface it seems to be more stable, FWIW. Perhaps one could also add a term (GroupID | StudentID), but I don't know offhand if that would make any sense. Another issue to consider is whether to fit by REML or ML, but that is secondary to getting the form of the model equation right. I don't claim this as a final model, but rather a 're-starting point'. It may well be in need of improvement, so comments are welcome. The confusion between subjects nested in time or vice versa has occurred several times this week with respect to repeated measures/ longitudinal models using lmer(), so perhaps it merits a comment: subjects/experimental units are NOT nested in time. Measurements taken on an individual at several time points *entails* that time be nested within subject. Just saying... This discussion may be better continued on the R-sig-mixed list, so I've cc-ed to that group as well. HTH, Dennis On Wed, Aug 25, 2010 at 1:27 AM, Bertolt Meyer wrote: Ben Bolker gmail.com> writes: Bertolt Meyer
[R] lattice help required
hello,
i want to stack two lattice plots beneath each other using one x-axis and
sharing the same text-panels,
like:
#
library(lattice)
y1 <- rnorm(100,100,10)
y2 <- rnorm(100,10,1)
facs<-expand.grid(Sites=rep(c("Site I","Site II"),25),Treatment=c("A","B"))
pl1<-dotplot(y1 ~ facs$Treatment|facs$Sites,
scales = list(x = list(rot = 90, tck=c(1,0
pl2<-dotplot(y2 ~ facs$Treatment|facs$Sites,
scales = list(x = list(rot = 90, tck=c(1,0
print(pl1, split=c(1,2,1,2), more=TRUE)
print(pl2, split=c(1,1,1,2))
#
but as said, ideally the plots should be stacked with only the lower plot
giving the x-axis annotation
and only the upper plot with text-panels.
thanks a lot,
kay
-
Kay Cichini
Postgraduate student
Institute of Botany
Univ. of Innsbruck
--
View this message in context:
http://r.789695.n4.nabble.com/lattice-help-required-tp2338382p2338382.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] frequency, count rows, data for heat map
Hi all, I have read posts of heat map creation but I am one step prior -- Here is what I am trying to do and wonder if you have any tips? We are trying to map sequence reads from tumors to viral genomes. Example input file : 111 abc 111 sdf 111 xyz 1079 abc 1079 xyz 1079 xyz 5576 abc 5576 sdf 5576 sdf How may xyz's are there for 1079 and 111? How many abc's, etc? How many times did reads from sample (1079) align to virus xyz. In some cases there are thousands per virus in a give sample, sometimes one. The original file (two columns by tens of thousands of rows; 20 MB) is text file (tab delimited). Output file: abc sdf xyz 111 1 1 1 1079 1 0 2 5576 1 2 0 Or, other ways to generate this data so I can then use it for heat map creation? Thanks for any help you may have, rtsweeney palo alto, ca -- View this message in context: http://r.789695.n4.nabble.com/frequency-count-rows-data-for-heat-map-tp2338363p2338363.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] SEM : Warning : Could not compute QR decomposition of Hessian
Hi useRs,
I'm trying for the first time to use a sem. The model finally runs,
but gives a warning saying :
"In sem.default(ram = ram, S = S, N = N, param.names = pars, var.names
= vars, : Could not compute QR decomposition of Hessian.
Optimization probably did not converge. "
I found in R-help some posts on this warning, but my attemps to modify
the code didn't change the warning message (i tried to give an error
of 1 to the latente variables). I can't figure what the problem is.
Here is the code :
tab<-read.table("F:/Mes documents/stats/sem/donnees_corr.txt",
header=T, sep="",na.strings = "NA")
tab[,46]<-as.factor(tab[,46])
tab[,24]<-as.factor(tab[,24])
tab[,40]<-as.factor(tab[,40])
fct_cor<-hetcor(tab, ML=T)
cor_tab<- fct_cor$correlations
moment_tab<-read.moments(diag=F, names=c('c1','c2', 'c3','c4','c5',
'c6','c7', 'c8', 'c9', 'ind_plando', 'long_sup15', 'long_inf15',
'pente', 'est', 'sud','ouest', 'nord' ,'reg_hydriq', 'prof_sol',
'pierro', 'efferv', 'struct','drainage','texture', 'route1_pond',
'route2_pond',
'pourcactif', 'tx_chomage', 'pourcagric', 'pourc_jeunes', 'pop99',
'rev_imp_foyer','eq_CONC', 'eq_sante', 'eq_edu', 'sold_nat',
'sold_mig', 'tps_dom_emp','TXEMPLOI','ORIECO','dist_paris','axe1',
'axe2', 'axe3', 'dist_protect','urbanisation','pays_incli','pays_alti'))
# after comes the moment matrix (triangular)
ram_tab<-specify.model()
type_paysage->pays_alti,NA,1
type_paysage->pays_incli, pays2, NA
pedo->reg_hydriq, NA, 1
pedo->prof_sol, ped8, NA
pedo->pierro, ped9, NA
pedo->efferv, ped10, NA
pedo->struct, ped11, NA
pedo->drainage, ped12, NA
pedo->texture, ped13, NA
adj_99->c1, NA,1
adj_99->c2, adj2,NA
adj_99->c3, adj3,NA
adj_99->c4, adj4,NA
adj_99->c5, adj5,NA
adj_99->c6, adj6,NA
adj_99->c7, adj7,NA
adj_99->c8, adj8,NA
adj_99->c9, adj9,NA
etat_hexa->axe1, NA, 1
etat_hexa->axe2, et2, NA
etat_hexa->axe3, et3, NA
socioBV->sold_mig, BV1, NA
socioBV->sold_nat, BV2, NA
socioBV->TXEMPLOI, BV3, NA
socioBV->ORIECO, BV4, NA
socioBV->tps_dom_emp, NA, 1
eqBV->eq_CONC, NA, 1
eqBV->eq_sante, eq2, NA
eqBV->eq_edu, eq3, NA
socio_com->pourcactif , NA, 1
socio_com->tx_chomage, com2, NA
socio_com->pourcagric, com3, NA
socio_com->pourc_jeunes, com4, NA
socio_com->pop99, com5, NA
socio_com->rev_imp_foyer, com7, NA
access_hexa->route1_pond, NA, 1
access_hexa->route2_pond, acc2, NA
hydro->ind_plando, NA, 1
hydro->long_sup15, eau2, NA
hydro->long_inf15, eau3, NA
topog->pente, NA, 1
topog->est, top2, NA
topog->sud, top3, NA
topog->nord, top4, NA
topog->ouest, top5, NA
dist_protect-> urbanisation, cor1,NA
dist_protect-> adj_99, cor2, NA
dist_protect-> etat_hexa, cor3, NA
topog-> urbanisation, cor4, NA
topog-> adj_99, cor5, NA
topog-> etat_hexa, cor6, NA
topog-> access_hexa, cor7, NA
topog<->hydro, cor8, NA
topog<->pedo, cor9, NA
pedo-> urbanisation, cor10, NA
pedo-> adj_99, cor11, NA
pedo-> etat_hexa, cor12, NA
pedo<->hydro, cor1, NA
hydro-> urbanisation, cor13, NA
hydro-> adj_99, cor14, NA
hydro-> etat_hexa, cor15, NA
access_hexa-> urbanisation, cor16, NA
access_hexa-> etat_hexa, cor17, NA
socio_com-> etat_hexa, cor18, NA
socio_com-> adj_99, cor19, NA
socio_com-> urbanisation, cor20, NA
dist_paris-> socio_com, cor21, NA
dist_paris-> access_hexa, cor22, NA
dist_paris-> adj_99, cor23, NA
dist_paris-> etat_hexa, cor24, NA
dist_paris-> urbanisation, cor25, NA
dist_paris-> socioBV, cor26, NA
socioBV-> eqBV, cor27, NA
socioBV-> urbanisation, cor28, NA
socioBV-> adj_99, cor29, NA
socioBV-> etat_hexa, cor30, NA
eqBV-> etat_hexa, cor31, NA
eqBV-> adj_99, cor32, NA
eqBV-> urbanisation, cor33, NA
etat_hexa-> urbanisation, cor34, NA
etat_hexa<-> adj_99, cor35, NA
adj_99-> urbanisation, cor36, NA
type_paysage-> urbanisation, cor37, NA
type_paysage-> adj_99, cor38, NA
type_paysage-> etat_hexa, cor39, NA
dist_paris<->dist_paris, auto1, NA
dist_protect<->dist_protect, auto2, NA
c1 <-> c1, auto4, NA
c2 <-> c2 , auto5, NA
c3 <-> c3 , auto6, NA
c4 <-> c4 , auto7, NA
c5 <-> c5 , auto8, NA
c6 <-> c6 , auto9, NA
c7 <-> c7 , auto10, NA
c8 <-> c8 , auto11, NA
c9 <-> c9 , auto12, NA
ind_plando <-> ind_plando, auto13, NA
long_sup15 <-> long_sup15 , auto14, NA
long_inf15 <-> long_inf15 , auto15, NA
pente <-> pente , auto16, NA
est<-> est , auto17, NA
sud <-> sud , auto18, NA
ouest<-> ouest , auto19, NA
nord <-> nord , auto20, NA
reg_hydriq <-> reg_hydriq , auto21, NA
prof_sol<-> prof_sol , auto22, NA
pierro <-> pierro , auto23, NA
efferv <-> efferv , auto24, NA
struct <-> struct , auto25, NA
drainage <-> drainage, auto26, NA
texture <-> texture , auto27, NA
route1_pond <->route1_pond , auto30, NA
route2_pond <-> route2_pond , auto31, NA
pourcactif <-> pourcactif , auto32, NA
tx_chomage <-> tx_chomage, auto33, NA
pourcagric <-> pourcagric , auto34, NA
pourc_jeunes<-> pourc_jeunes , auto36, NA
pop99<-> pop99 , auto36, NA
rev_imp_foyer <-> rev_imp_foyer , auto38, NA
eq_CONC <-> eq_CONC, auto39, NA
eq_sante <->eq_sante , auto40, NA
eq_edu<->eq_edu , auto41, NA
sold_nat<-> sold_nat , aut
Re: [R] Comparing samples with widely different uncertainties
On Aug 25, 2010, at 8:57 AM, Sandy Small wrote: > Hi > This is probably more of a statistics question than a specific R > question, although I will be using R and need to know how to solve the > problem in R. > > I have several sets of data (ejection fraction measurements) taken in > various ways from the same set of (~400) patients (so it is paired data). > For each individual measurement I can make an estimate of the percentage > uncertainty in the measurement. > Generally the measurements in data set A are higher but they have a > large uncertainty (~20%) while the measurements in data set Bare lower > but have a small uncertainty (~4%). > I believe, from the physiology, that the true value is likely to be > nearer the value of A than of B. > I need to show that, despite the uncertainties in the measurements > (which are not themselves normally distributed), there is (or is not) a > difference between the two groups, (a straight Wilcoxon signed ranks > test shows a difference but it cannot include that uncertainty data). > > Can anybody suggest what I should be looking at? Is there a language > here that I don't know? How do I do it in R? > Many thanks for your help > Sandy The first place that I would start is at Martin Bland's page pertaining to the design and analysis of measurement studies: http://www-users.york.ac.uk/~mb55/meas/meas.htm The papers he co-authored with Doug Altman are the "go to" resource for this domain. HTH, Marc Schwartz __ 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] Comparing samples with widely different uncertainties
On Aug 25, 2010, at 3:57 PM, Sandy Small wrote: > Hi > This is probably more of a statistics question than a specific R > question, although I will be using R and need to know how to solve the > problem in R. > > I have several sets of data (ejection fraction measurements) taken in > various ways from the same set of (~400) patients (so it is paired data). > For each individual measurement I can make an estimate of the percentage > uncertainty in the measurement. > Generally the measurements in data set A are higher but they have a > large uncertainty (~20%) while the measurements in data set Bare lower > but have a small uncertainty (~4%). > I believe, from the physiology, that the true value is likely to be > nearer the value of A than of B. > I need to show that, despite the uncertainties in the measurements > (which are not themselves normally distributed), there is (or is not) a > difference between the two groups, (a straight Wilcoxon signed ranks > test shows a difference but it cannot include that uncertainty data). > > Can anybody suggest what I should be looking at? Is there a language > here that I don't know? How do I do it in R? > Many thanks for your help > Sandy Hm, well... I don't think the issue is entirely well-defined, but let me try and give you some pointers: For bivariate normal data (X,Y), the situation is that even if V(X) != V(Y) you still end up looking at X-Y if the question is whether the means are the same. It's sort of the only thing you _can_ do... For non-normal data, it is not clear what the null hypothesis really is. The signed-rank test assumes that X-Y has a symmetric distribution, which is dubious if X is not symmetric and its variation dominates that of Y. You could also do a sign test and see if the differences has a median of zero (bear in mind that the median of a difference is different from the difference of the medians, but it could actually suffice.) I'd probably start off with a simple plot of Y vs X and look for fan-out effects indicating that the variance depends on the mean. If it does, perhaps take logs or square roots and see if it makes the variances appear more stable and perhaps improve on the normality. Then maybe just do a paired t-test. -- Peter Dalgaard Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark 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] approxfun-problems (yleft and yright ignored)
Dear all,
I have run into a problem when running some code implemented in the
Bioconductor panp-package (applied to my own expression data), whereby gene
expression values of known true negative probesets (x) are interpolated onto
present/absent p-values (y) between 0 and 1 using the *approxfun -
function*{stats}; when I have used R version 2.8, everything had
worked fine,
however, after updating to R 2.11.1., I got unexpected output (explained
below).
Please correct me here, but as far as I understand, the yleft and yright
arguments set the extreme values of the interpolated y-values in case the
input x-values (on whose approxfun is applied) fall outside range(x). So if
I run approxfun with yleft=1 and yright=0 with y-values between 0 and 1,
then I should never get any values higher than 1. However, this is not the
case, as this code-example illustrates:
> ### define the x-values used to construct the approxfun, basically these
are 2000 expression values ranging from ~ 3 to 7:
> xNeg <- NegExprs[, 1]
> xNeg <- sort(xNeg, decreasing = TRUE)
>
> ### generate 2000 y-values between 0 and 1:
> yNeg <- seq(0, 1, 1/(length(xNeg) - 1))
> ### define yleft and yright as well as the rule to clarify what should
happen if input x-values lie outside range(x):
> interp <- approxfun(xNeg, yNeg, yleft = 1, yright = 0, rule=2)
Warning message:
In approxfun(xNeg, yNeg, yleft = 1, yright = 0, rule = 2) :
collapsing to unique 'x' values
> ### apply the approxfun to expression data that range from ~2.9 to 13.9
and can therefore lie outside range(xNeg):
> PV <- sapply(AllExprs[, 1], interp)
> range(PV)
[1]0.000 6208.932
> summary(PV)
Min. 1st Qu.Median Mean 3rd Qu. Max.
0.000e+00 0.000e+00 2.774e-03 1.299e+00 3.164e-01 6.209e+03
So the resulting output PV object contains data ranging from 0 to 6208, the
latter of which lies outside yleft and is not anywhere close to extreme
y-values that were used to set up the interp-function. This seems wrong to
me, and from what I understand, yleft and yright are simply ignored?
I have attached a few histograms that visualize the data distributions of
the objects I xNeg, yNeg, AllExprs[,1] (== input x-values) and PV (the
output), so that it is easier to make sense of the data structures...
Does anyone have an explanation for this or can tell me how to fix the
problem?
Thanks a million for any help, best, Sam
> sessionInfo()
R version 2.11.1 (2010-05-31)
x86_64-apple-darwin9.8.0
locale:
[1] en_IE.UTF-8/en_IE.UTF-8/C/C/en_IE.UTF-8/en_IE.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] panp_1.18.0 affy_1.26.1 Biobase_2.8.0
loaded via a namespace (and not attached):
[1] affyio_1.16.0 preprocessCore_1.10.0
--
-
Samuel Wuest
Smurfit Institute of Genetics
Trinity College Dublin
Dublin 2, Ireland
Phone: +353-1-896 2444
Web: http://www.tcd.ie/Genetics/wellmer-2/index.html
Email: wue...@tcd.ie
--
__
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] Odp: Finding pairs
Hi well, I will add some explanation r-help-boun...@r-project.org napsal dne 25.08.2010 11:24:38: > Dear Mr Petr PIKAL > After reading the R code provided by you, I realized that I would have never > figured out how this could have been done. I am going to re-read again and > again your code to understand the logic and the commands you have provided. > Thanks again from the heart for your kind advice. > Regards > Mike > > --- On Wed, 25/8/10, Petr PIKAL wrote: > > From: Petr PIKAL > Subject: Re: [R] Odp: Finding pairs > To: "Mike Rhodes" > Cc: r-help@r-project.org > Date: Wednesday, 25 August, 2010, 9:01 > > Hm > > r-help-boun...@r-project.org napsal dne 25.08.2010 09:43:26: > > > Dear Mr Petr Pikal > > > > I am extremely sorry for the manner I have raised the query. Actually > that was > > my first post to this R forum and in fact even I was also bit confused > while > > drafting the query, for which I really owe sorry to all for consuming > the > > precious time. Perhaps I will try to redraft my query in a better way as > follows. > > > > I have two datasets "A" and "B" containing the names of branch offices > of a > > particular bank say XYZ plc bank. The XYZ bank has number of main branch > > > offices (say Parent) and some small branch offices falling under the > purview > > of these main branch offices (say Child). > > > > The datalist "A" and "B" consists of these main branch office names as > well as > > small branch office names. B is subset of A and these branch names are > coded. > > Thus we have two datasets A and B as (again I am using only a > > portion of a large database just to have some idea) > > > > > > A B > > 144 >what is here in B? Empty space?, > > 145 > > 146 > > 147 144 > > How do you know that 144 from B relates to 147 in A? Is it according to > its positions? I.e. 4th item in B belongs to 4.th item in A? > > > 148 145 > > > > 149 147 > > 151 148 > > > > > > > > Now the branch 144 appears in A as well as in B and in B it is mapped > with > > 147. This means branch 147 comes under the purview of main branch 144. > Again > > 147 is controlling the branch 149 (since 147 also has appeared in B and > is > > mapped with 149 of A). > > > > Similarly, branch 145 is controlling branch 148 which further controls > > operations of bank branch 151 and like wise. > > Well as you did not say anything about structure of your data > A<-144:151 > B<-144:148 > data.frame(A,B) > A B > 1 144 NA > 2 145 NA > 3 146 NA > 4 147 144 > 5 148 145 > 6 149 146 > 7 150 147 > 8 151 148 > DF<-data.frame(A,B) This was just making a data frame with 2 columns to have some data to play with > main<-DF$A[is.na(DF$B)] Above are codes from A which are NA in B > branch1<-DF[!is.na(DF$B),] Above is data frame of remaining codes (other than main) > selected.branch1<-branch1$A[branch1$B%in%main] Above is codes from column A for which B column and main are the same > branch2<-branch1[!branch1$B%in%main,] This is the rest of yet not selected rows > selected.branch2<-branch2$A[branch2$B%in%selected.branch1] and this is selection of values from column A for which B column and selected.branch1 values are same. But it works for this particular data, I am not sure how it behaves with duplicates and further issues. It also depends on how your data is organised. And if you are in reading you could also go through setdiff, merge and maybe sqldf package and Rdata Import/export manual Regards Petr > > and for cbinding your data which has uneven number of values see Jim > Holtman's answer to this > > How to cbind DF:s with differing number of rows? > > Regards > Petr > > > > > > So in the end I need an output something like - > > > > Main Branch Branch office1 Branch > > office2 > > 144 147 149 > > 145 148 151 > > > 146 NA > > NA > > > ... > > > .. > > > > > > I understand again I am not able to put forward my query properly. But I > must > > thank all of you for giving a patient reading to my query and for > reverting > > back earlier. Thanks once again. > > > > With warmest regards > > > > Mike > > > > > > --- On Wed, 25/8/10, Petr PIKAL wrote: > > > > From: Petr PIKAL > > Subject: Odp: [R] Finding > > pairs > > To: "Mike Rhodes" > > Cc: r-help@r-project.org > > Date: Wednesday, 25 August, 2010, 6:39 > > > > Hi > > > > wit
Re: [R] how do R calculates the number of intervals between tick-marks
Take a look at pretty function. On Tue, Aug 24, 2010 at 8:05 PM, Antonio Olinto wrote: > Hello, > > I want to know how do R calculates the number of intervals between > tick-marks in the y axis in a plot. > > I'm making a three y-axes plot and this information would help me a lot. > > Thanks in advance. > > Antonio Olinto > > > > > Webmail - iBCMG Internet > http://www.ibcmg.com.br > > __ > 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] Repeat the first day data through all the day. Zoo
On Wed, Aug 25, 2010 at 9:48 AM, skan wrote: > > thanks > I'll try them, > > Why do you use the brackets in zz[] ? So it stays a zoo object with the same index. We are only replacing the data part. __ 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] Comparing samples with widely different uncertainties
Hi
This is probably more of a statistics question than a specific R
question, although I will be using R and need to know how to solve the
problem in R.
I have several sets of data (ejection fraction measurements) taken in
various ways from the same set of (~400) patients (so it is paired data).
For each individual measurement I can make an estimate of the percentage
uncertainty in the measurement.
Generally the measurements in data set A are higher but they have a
large uncertainty (~20%) while the measurements in data set Bare lower
but have a small uncertainty (~4%).
I believe, from the physiology, that the true value is likely to be
nearer the value of A than of B.
I need to show that, despite the uncertainties in the measurements
(which are not themselves normally distributed), there is (or is not) a
difference between the two groups, (a straight Wilcoxon signed ranks
test shows a difference but it cannot include that uncertainty data).
Can anybody suggest what I should be looking at? Is there a language
here that I don't know? How do I do it in R?
Many thanks for your help
Sandy
--
Sandy Small
Clinical Physicist
NHS Greater Glasgow and Clyde
and
NHS Forth Valley
Phone: 01412114592
E-mail: sandy.sm...@nhs.net
This message may contain confidential information. If yo...{{dropped:21}}
__
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 do R calculates the number of intervals between tick-marks
On Aug 24, 2010, at 7:05 PM, Antonio Olinto wrote: Hello, I want to know how do R calculates the number of intervals between tick-marks in the y axis in a plot. ?axTicks # and then look at the other citations and the code as needed I'm making a three y-axes plot and this information would help me a lot. Thanks in advance. Antonio Olinto -- David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Repeat the first day data through all the day. Zoo
thanks I'll try them, Why do you use the brackets in zz[] ? -- View this message in context: http://r.789695.n4.nabble.com/Repeat-the-first-day-data-through-all-the-day-Zoo-tp2338069p2338266.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] percentage sign in expression
On Aug 25, 2010, at 4:32 AM, e-letter wrote:
On 24/08/2010, David Winsemius wrote:
On Aug 24, 2010, at 9:37 AM, e-letter wrote:
Readers,
According to the documentation for the function 'plotmath' there
is no
apparent possibility to add the percent sign (%) to a plot function,
Where did you see an assertion made???
Within R I entered the command:
?plotmath
Also accessed using:
help.start(browser="opera")
Navigating the web browser page:
packages
packages in /usr/lib/R/library
grdevices
plotmath
In the list headed 'syntax' and 'meaning' within the section
'details'.
e.g.
plot(a[,1]~b[,2],ylab=expression(x~%),xlab=expression(z))
How to achieve this please?
Read the plotmath helo page more carefully. The section immediatedly
below the plotmath expressions points you to the use of the symbol()
expression-function and to the points help page where generation of
the available glyphs proceeds according to the advice on
help(plotmath):
In my system the paragraph immediately after the list of features
(i.e. 'syntax','meaning') describes a note to TeX users. I cannot see
reference to 'symbol()'.
It's possible that my help page is different than yours. Right after
the syntax/meaning description on mine (which is a Mac OSX system) is
a paragraph:
"The symbol font uses Adobe Symbol encoding so, for example, a lower
case mu can be obtained either by the special symbol mu or by
symbol("m"). This provides access to symbols that have no special
symbol name, for example, the universal, or forall, symbol is
symbol("\042"). To see what symbols are available in this way
useTestChars(font=5) as given in the examples for points: some are
only available on some devices."
(In this case I would be surprised if the help pages were different
because this makes a cross-reference to the examples in points. I am
not surprised about cross-platform differences in descriptions of
graphical devices and would have included a caveat if I were
corresponding on rhelp about such. I suppose the font issues could be
platform specific so if you want to correct me on this point, I will
try to file it away. I did, however, give you the code needed to to
display Symbols and it sounds further on that it succeeded)
TestChars <- function(sign=1, font=1, ...)
+ {
+if(font == 5) { sign <- 1; r <- c(32:126, 160:254)
+} else if (l10n_info()$MBCS) r <- 32:126 else r <- 32:255
+if (sign == -1) r <- c(32:126, 160:255)
+par(pty="s")
+plot(c(-1,16), c(-1,16), type="n", xlab="", ylab="",
+ xaxs="i", yaxs="i")
+grid(17, 17, lty=1)
+for(i in r) try(points(i%%16, i%/%16, pch=sign*i,
font=font,...))
+ }
TestChars(font=5)
Notice that the "%" sign is three characters to the right (i.e.
higher) of the "forall" symbol that is produced by the example code
I can't see 'forall' in the code above.
Gavin has already explained why you did not. The upside-down A (==
"universal" or "forall") was a useful reference point in the indexing,
since it is only 3 glyphs away for the "%" symbol.
they offer. (The numbering proceeds from bottom up which confused me
at first.)
What numbering?
Actually I did not see any numbering either, which was why I remained
confused about the location of the "%" symbol for several minutes.
Perhaps I should have used the term "indexing".
The documentation makes reference to the command:
demo(plotmath)
I applied this command and could not see an instruction to produce the
percent (%) symbol.
I don't think I suggested it would.
--
David Winsemius, MD
West Hartford, CT
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] multiple x factors using the sciplot package
Dear R community, I am a beginner using the sciplot package to graph barplots. I would like to be able to graph two x factors (Sampling.year and Period). Sampling.year: 2006, 2007, 2008, 2009 Period: First, Second, Total The parameter "group" is the different species I looked at. They can be seen in the legend. The parameter "response" is the Average percentage cover category of these species. I would like the graph so you see on the x axis the years (in a bigger font) and the period (in smaller font). It would look a bit like the barplot that can be seen at http://onertipaday.blogspot.com/2007/05/make-many-barplot-into-one-plot.html, but with the advantage of the sciplot package. Thanks a lot for your help! Jordan [[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] multi day R run crashed - is there a log
Hi , thanks, the lower.tri idea is I guess the best way. Will try that. Cheers Martin On 8/25/2010 2:21 PM, Sarah Goslee wrote: Hi, On Wed, Aug 25, 2010 at 8:11 AM, Martin Tomko wrote: Hi Sarah, thank you very much for your answer. I have been spitting things out on screen, but unfortunately I have not run it as a batch log, but in the interactive window, so when the GUI crashed I was left without trace... I guess I should google how to run it from a batch. I have no idea how to do that in Windows, but I'm sure it's possible. :) That way the things written to the screen will be saved in a file instead. I should also explore using RData, I have ony been using csv files and tables so far, it seems that can also bring some added performance. If you need to save *everything*, RData is what you get when you use save(), or when you close a session and choose y to saving the data. It can be read in using load(). That's one way to be able to pick up where you left off. As the main output of my process is a matrix, I would really need to append to a matrix after each iteration. I have identified the write.table append parameter-based solution, but that would only append rows. Is there a way to slowly "grow" a matrix in both directions, meaning append columns as well (it is a big distance matrix). AFAIK, you can't append columns that way because of the way text files are written to disk. You'd need to rewrite the whole thing, or possibly write it out in lower triangular format with NA values as padding (assuming it's a symmetric distance). Or for that matter, you could just write it out as a really long vector, and turn it back into a matrix later if you need to read the saved file in after a crash. I'd recommend saving whatever variables are needed so that you can pick up exactly where you left off, if possible. Much nicer to pick up 12 hours in than to start over from the beginning. Not R, but I just finished a 5-week batch job. You can bet that I put a lot of thought into incremental save points and how to resume after an unexpected halt! Sarah -- Martin Tomko Postdoctoral Research Assistant Geographic Information Systems Division Department of Geography University of Zurich - Irchel Winterthurerstr. 190 CH-8057 Zurich, Switzerland email: martin.to...@geo.uzh.ch site: http://www.geo.uzh.ch/~mtomko mob:+41-788 629 558 tel:+41-44-6355256 fax:+41-44-6356848 __ 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.
