Re: [R] COVID-19 datasets...

[Thomas Petzoldt]

On 07.05.2020 at 13:12 Deepayan Sarkar wrote:

On Thu, May 7, 2020 at 4:16 PM Thomas Petzoldt  wrote:

On 07.05.2020 at 11:19 Deepayan Sarkar wrote:

On Thu, May 7, 2020 at 12:58 AM Thomas Petzoldt  wrote:

Sorry if I'm joining a little bit late.

I've put some related links and scripts together a few weeks ago. Then I
stopped with this, because there is so much.

The data format employed by John Hopkins CSSE was sort of a big surprise
to me.

Why? I find it quite convenient to drop the first few columns and
extract the data as a matrix (using data.matrix()).

-Deepayan

Many thanks for the hint to use data.matrix

My aim was not to say that it is difficult, especially as R has all the
tools for data mangling.

My surprise was that "wide tables" and non-ISO dates as column names are
not the "data base way" that we in general teach to our students

Well, I am all for long format data when it makes sense, but I would
disagree that that is always the "right approach". In the case of
regular multiple time series, as in this context, a matrix-like
structure seems much more natural (and nicely handled by ts() in R),
and I wouldn't even bother reshaping the data in the first place.

See, for example,

https://github.com/deepayan/deepayan.github.io/blob/master/covid-19/deaths.rmd

and

https://deepayan.github.io/covid-19/deaths.html

-Deepayan


Great, thank you for the link with the comprehensive lattice graphs and 
the explanations. I like your package very much and use it often, since 
it appeared on CRAN (3 of my CRAN packages depend on it). As "dynamic 
modeller", I consider time always as the first column, but I agree on 
the other hand, that long tables are often, but not always the right 
approach, let's think about gridded multi-dimensional netcdf data.


Many thanks for sharing your analysis publicly, I'll add your repo to my 
link list.


Thomas


With reshape2::melt or tidyr::gather resp. pivot_longer, conversion is
quite easy, regardless if one wants to use tidyverse or not, see example
below.

Again, thanks, Thomas


library("dplyr")
library("readr")
library("tidyr")

file <-
"https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv;

dat <- read_delim(file, delim=",")
names(dat)[1:2] <- c("Province_State", "Country_Region")
dat2 <-
dat %>%
## summarize Country/Region duplicates
group_by(Country_Region) %>% summarise_at(vars(-(1:4)), sum) %>%
## make it a long table
pivot_longer(cols = -Country_Region, names_to = "time") %>%
## convert to ISO 8601 date
mutate(time = as.POSIXct(time, format="%m/%e/%y"))




An opposite approach was taken in Germany, that organized it as a
big JSON trees.

Fortunately, both can be "tidied" with R, and represent good didactic
examples for our students.

Here yet another repo linking to the data:

https://github.com/tpetzoldt/covid


Thomas


On 04.05.2020 at 20:48 James Spottiswoode wrote:

Sure. COVID-19 Data Repository by the Center for Systems Science and 
Engineering (CSSE) at Johns Hopkins University is available here:

https://github.com/CSSEGISandData/COVID-19

All in csv fiormat.



On May 4, 2020, at 11:31 AM, Bernard McGarvey  
wrote:

Just curious does anyone know of a website that has data available in a format 
that R can download and analyze?

Thanks


Bernard McGarvey


Director, Fort Myers Beach Lions Foundation, Inc.


Retired (Lilly Engineering Fellow).

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


James Spottiswoode
Applied Mathematics & Statistics
(310) 270 6220
jamesspottiswoode Skype
ja...@jsasoc.com


--
Dr. Thomas Petzoldt
senior scientist

Technische Universitaet Dresden
Faculty of Environmental Sciences
Institute of Hydrobiology
01062 Dresden, Germany

https://tu-dresden.de/Members/thomas.petzoldt


__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] COVID-19 datasets...

[Thomas Petzoldt]

On 07.05.2020 at 11:19 Deepayan Sarkar wrote:

On Thu, May 7, 2020 at 12:58 AM Thomas Petzoldt  wrote:


Sorry if I'm joining a little bit late.

I've put some related links and scripts together a few weeks ago. Then I
stopped with this, because there is so much.

The data format employed by John Hopkins CSSE was sort of a big surprise
to me.


Why? I find it quite convenient to drop the first few columns and
extract the data as a matrix (using data.matrix()).

-Deepayan


Many thanks for the hint to use data.matrix

My aim was not to say that it is difficult, especially as R has all the 
tools for data mangling.


My surprise was that "wide tables" and non-ISO dates as column names are 
not the "data base way" that we in general teach to our students


With reshape2::melt or tidyr::gather resp. pivot_longer, conversion is 
quite easy, regardless if one wants to use tidyverse or not, see example 
below.


Again, thanks, Thomas


library("dplyr")
library("readr")
library("tidyr")

file <- 
"https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv;


dat <- read_delim(file, delim=",")
names(dat)[1:2] <- c("Province_State", "Country_Region")
dat2 <-
  dat %>%
  ## summarize Country/Region duplicates
  group_by(Country_Region) %>% summarise_at(vars(-(1:4)), sum) %>%
  ## make it a long table
  pivot_longer(cols = -Country_Region, names_to = "time") %>%
  ## convert to ISO 8601 date
  mutate(time = as.POSIXct(time, format="%m/%e/%y"))






An opposite approach was taken in Germany, that organized it as a
big JSON trees.

Fortunately, both can be "tidied" with R, and represent good didactic
examples for our students.

Here yet another repo linking to the data:

https://github.com/tpetzoldt/covid


Thomas


On 04.05.2020 at 20:48 James Spottiswoode wrote:

Sure. COVID-19 Data Repository by the Center for Systems Science and 
Engineering (CSSE) at Johns Hopkins University is available here:

https://github.com/CSSEGISandData/COVID-19

All in csv fiormat.



On May 4, 2020, at 11:31 AM, Bernard McGarvey  
wrote:

Just curious does anyone know of a website that has data available in a format 
that R can download and analyze?

Thanks


Bernard McGarvey


Director, Fort Myers Beach Lions Foundation, Inc.


Retired (Lilly Engineering Fellow).

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



James Spottiswoode
Applied Mathematics & Statistics
(310) 270 6220
jamesspottiswoode Skype
ja...@jsasoc.com



--
Dr. Thomas Petzoldt
senior scientist

Technische Universitaet Dresden
Faculty of Environmental Sciences
Institute of Hydrobiology
01062 Dresden, Germany

https://tu-dresden.de/Members/thomas.petzoldt

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] COVID-19 datasets...

[Thomas Petzoldt]

Sorry if I'm joining a little bit late.

I've put some related links and scripts together a few weeks ago. Then I 
stopped with this, because there is so much.


The data format employed by John Hopkins CSSE was sort of a big surprise 
to me. An opposite approach was taken in Germany, that organized it as a 
big JSON trees.


Fortunately, both can be "tidied" with R, and represent good didactic 
examples for our students.


Here yet another repo linking to the data:

https://github.com/tpetzoldt/covid


Thomas


On 04.05.2020 at 20:48 James Spottiswoode wrote:

Sure. COVID-19 Data Repository by the Center for Systems Science and 
Engineering (CSSE) at Johns Hopkins University is available here:

https://github.com/CSSEGISandData/COVID-19

All in csv fiormat.



On May 4, 2020, at 11:31 AM, Bernard McGarvey  
wrote:

Just curious does anyone know of a website that has data available in a format 
that R can download and analyze?
  
Thanks



Bernard McGarvey


Director, Fort Myers Beach Lions Foundation, Inc.


Retired (Lilly Engineering Fellow).

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



James Spottiswoode
Applied Mathematics & Statistics
(310) 270 6220
jamesspottiswoode Skype
ja...@jsasoc.com




__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] Partial differential equation

[Thomas Petzoldt]

On 19.12.2017 04:22, Timothy Axberg wrote:

Hello, I am having troubles with heat conduction problem. Below is the
given information. Should I move forward with this problem like any other
1-D PDE?


Yes. You may use your favorite search engine and search for:

"heat transfer" "deSolve" R

or have a look at:

http://desolve.r-forge.r-project.org



¶T/¶t = a* (¶^2T/¶x^2)



I.C. For t = 0 and 0 £ x £ 10, T = 0 °C

B.C. For x = 0 cm and all t , T = 100°C

For x = 10 cm and all t , T = 0 °C

The thermal diffusivity is a = 2.0 cm^2 /s.


I also added what I have for code in R. Any help will be appreciated. Thanks


There was no code appended.


- Timothy


Hope it helps,

Thomas

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] Date operation Question in R

[Thomas Petzoldt]

On 30.03.2017 23:34, Paul Bernal wrote:

Hello everyone,

Is there a way to use the function seq to generate a date sequence in
this kind of format: jan-2007?


format(seq(ISOdate(2017,1,1), ISOdate(2017,12,31), "months"), "%b-%Y")



Also, is there a way to change the Sys.Date() format to the one
mentioned above (jan-2007)?


format(Sys.Date(), "%b-%Y")

see ?strptime for details.

Thomas



Thanks in advance for your valuable help,

Best regards,

Paul


__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] non-linear optimisation ODE models

[Thomas Petzoldt]

Hi,

fitting ODE models may also be done with package FME, see:

Soetaert K, Petzoldt T. Inverse modelling, sensitivity and Monte Carlo 
analysis in R using package FME. Journal of Statistical Software. 
2010(33): 1–28. http://dx.doi.org/10.18637/jss.v033.i03


or the (interactive) poster at:

http://desolve.r-forge.r-project.org/user2014/examples/FME/fit_twocomp.svg

Regards,

Thomas P.


Am 16.02.2017 um 19:41 schrieb David Winsemius:



On Feb 15, 2017, at 1:43 PM, Jim Lemon  wrote:

Hi Malgorzata,
The function "rxnrate" seems to want three values in a list with the
names "k1", "k2" and "k3". If you are passing something with different
names, it is probably going to complain, so the names "A", "B" and "C"
may be your problem. I can't run the example, so this is a guess.


There's a more readable version at:
http://stackoverflow.com/questions/42256509/how-to-feed-data-into-ide-while-doing-optimisation

It can be run, but does not produce the errors offered when I do so.



__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] Differential equations

[Thomas Petzoldt]

Hi,

yes, the suggestion of Peter Dalgaard is correct, and it can also be 
done using the "event" mechanism of deSolve, see:


library("deSolve")
?events

... or the slides from our talk given at useR-2011

http://desolve.r-forge.r-project.org/slides/petz_soet2011.pdf

... or the tutorial from L.A.

http://desolve.r-forge.r-project.org/user2014/tutorial.pdf

Thomas Petzoldt

See also:

http://desolve.r-forge.r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models


Am 08.02.2017 um 16:19 schrieb peter dalgaard:

It's been a while, but I think I have gotten through this sort of situation by 
splitting the integration into intervals, i.e., you run from t=0 to t=20 witn 
initial condition c(100,0), yielding a value c(y1,y2), then you run from 20 to 
40 with initial condition c(y1+100, y2), etc.

-pd

On 08 Feb 2017, at 11:10 , Fanny Gallais <gallais.fa...@gmail.com> wrote:


Hi,

I'm working on a system of 2 differential equations. My initial condition
(t=0) is c(100,0) and i'm using lsoda function (from package deSolve) to
solve it.
My system reprensents the evoution of drug concentration in two
compartments throug time. Problem is I would like to model a repeated drug
administration. That is to say, not only 100 at t=0 but also at t=20,40,...
I can't find a solution to do so. I tried adding "100" to the first
differential equation at the times of interest but it doesn't work. Do you
have any idea?

Thank you
F.G.

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 x and y axis labels with the function plot.sensFun (package FME)

[Thomas Petzoldt]

Hi,

I reproduced your example and it worked for me with "Time" and 
"Population size" as axis labels. Are you using the most recent version?


Thomas



Am 25.01.2017 um 02:32 schrieb Marine Regis:

Hello,


How can I add x and y axis labels for a plot that is built from the function plot.sensFun (package 
FME) ? I tested xlab ="Time" and ylab="Population size" but this doesn't work 
for me.


Here is a code to generate the plot:


pars <- list(gmax = 0.5, eff = 0.5,
 ks = 0.5, rB = 0.01, dB = 0.01)

solveBact <- function(pars) {
  derivs <- function(t, state, pars) { # returns rate of change
with (as.list(c(state, pars)), {
  dBact <-  gmax * eff * Sub/(Sub + ks) * Bact - dB * Bact - rB * Bact
  dSub  <- -gmax   * Sub/(Sub + ks) * Bact + dB * Bact
  return(list(c(dBact, dSub)))
})
  }
  state   <- c(Bact = 0.1, Sub = 100)
  tout<- seq(0, 50, by = 0.5)
  ## ode solves the model by integration ...
  return(as.data.frame(ode(y = state, times = tout, func = derivs,
   parms = pars)))
}

## sensitivity functions
SF <- sensFun(func = solveBact, parms = pars,
  sensvar = c("Bact", "Sub"), varscale = 1)
par(mfrow = c(1,2))
plot(SF, which = c("Sub", "Bact"), mfrow = NULL, xlab="Time", ylab="Population 
size")

Thanks a lot for your time.
Marine


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] separate mixture of gamma and normal (with mixtools or ??)

[Thomas Petzoldt]

Dear Bert,

thank you very much for your suggestion. You are right, ill-conditioning 
was sometimes a problem for 3 components, but not in the two-component 
case. The modes are well separated, and the sample size is high.


My main problem is (1) the shape of the distributions and (2) the 
diversity of available packages and approaches to this topic.


In the mean time I made some progress in this matter by treating the 
data as a mixture of gamma distributions (package mixdist, see below), 
so what I want to know is the purely R technical question ;-)


Has someone else has ever stumbled across something like this and can 
make a suggestion which package to use?


Thanks, Thomas


## Approximate an Exponential+Gaussian mixture with
##   a mixture of Gammas

library("mixdist")
set.seed(123)
lambda <- c(0.25, 0.75)
N <- 2000
x <- c(rexp(lambda[1]*N, 1), rnorm(lambda[2]*N, 20, 4))

xx <- mixgroup(x, breaks=0:40)
pp <- mixparam(mu=c(1, 8), sigma=c(1, 3), pi=c(0.2, 0.5))
mix <-  mix(xx, pp, dist="gamma", emsteps=10)

summary(mix)
p <- coef(mix)
beta <- with(p, sigma^2/mu)
alpha <- with(p, mu /beta)
lambda <- p$pi

plot(mix, xlim=c(0, 35))
x1 <- seq(0, 35, 0.1)
lines(x1, lambda[1]*dgamma(x1, alpha[1], 1/beta[1]),
  col="orange", lwd=2)
lines(x1, lambda[2]*dgamma(x1, alpha[2], 1/beta[2]),
  col="magenta", lwd=2)



Am 24.01.2017 um 00:34 schrieb Bert Gunter:

Fitting multicomponent mixtures distributions -- and 3 is already a
lot of components -- is inherently ill-conditioned. You may need to
reassess your strategy. You might wish to post on stackexchange
instead to discuss such statistical issues.

Cheers,
Bert
Bert Gunter

"The trouble with having an open mind is that people keep coming along
and sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )


On Mon, Jan 23, 2017 at 2:32 PM, Thomas Petzoldt <t...@simecol.de> wrote:

Dear friends,

I am trying to separate bi- (and sometimes tri-) modal univariate mixtures
of biological data, where the first component is left bounded (e.g.
exponential or gamma) and the other(s) approximately Gaussian.

After checking several packages, I'm not really clear what to do. Here is an
example with "mixtools" that already works quite good, however, the left
component is not Gaussian (and not symmetric).

Any idea about a more adequate function or package for this problem?

Thanks a lot!

Thomas



library(mixtools)
set.seed(123)

lambda <- c(0.25, 0.75)
N  <- 200

## dist1 ~ gamma (or exponential as a special case)
#dist1 <- rexp(lambda[1]*N, 1)
dist1 <- rgamma(lambda[1]*N, 1, 1)

## dist2 ~ normal
dist2 <- rnorm(lambda[2]*N, 12, 2)

## mixture
x <- c(dist1, dist2)

mix <-  spEMsymloc(x, mu0=2, eps=1e-3, verbose=TRUE)
plot(mix, xlim=c(0, 25))
summary(mix)


--
Thomas Petzoldt
TU Dresden, Institute of Hydrobiology
http://www.tu-dresden.de/Members/thomas.petzoldt

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] separate mixture of gamma and normal (with mixtools or ??)

[Thomas Petzoldt]

Dear friends,

I am trying to separate bi- (and sometimes tri-) modal univariate 
mixtures of biological data, where the first component is left bounded 
(e.g. exponential or gamma) and the other(s) approximately Gaussian.


After checking several packages, I'm not really clear what to do. Here 
is an example with "mixtools" that already works quite good, however, 
the left component is not Gaussian (and not symmetric).


Any idea about a more adequate function or package for this problem?

Thanks a lot!

Thomas



library(mixtools)
set.seed(123)

lambda <- c(0.25, 0.75)
N  <- 200

## dist1 ~ gamma (or exponential as a special case)
#dist1 <- rexp(lambda[1]*N, 1)
dist1 <- rgamma(lambda[1]*N, 1, 1)

## dist2 ~ normal
dist2 <- rnorm(lambda[2]*N, 12, 2)

## mixture
x <- c(dist1, dist2)

mix <-  spEMsymloc(x, mu0=2, eps=1e-3, verbose=TRUE)
plot(mix, xlim=c(0, 25))
summary(mix)


--
Thomas Petzoldt
TU Dresden, Institute of Hydrobiology
http://www.tu-dresden.de/Members/thomas.petzoldt

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] Euler & Runge-Kutta

[Thomas Petzoldt]

Am 07.11.2016 um 21:52 schrieb Franklin Bretschneider:

Hello Tom Mosca,

Re:


Can someone help me with R code to perform approximations to second order 
differential equations and systems of first order differential equations using 
Euler's method and Runge-Kutta?  I am not a student and this is not for a test 
or graded assignment.

Examples (unrelated to each other):

h = 0.1

1.  3(t^2)y'' - 5ty' + 5y = 0
y(1) = 0, y'(1) = 2/3

2.  Lotka-Volterra
x' = x(3-y)
y' = y(x-3)




For solving differential equations, the famous "lsoda" solver is best, far better 
than the simple Euler method, and also better than Runge & Kutta.
Angels, or very nice people, implemented solving methods using the lsoda into R. Download 
the package "deSolve" and the world of numerical solutions is at your feet.


deSolve contains all of this: Euler, Runge-Kutta, lsoda and more ... and 
you will find many examples, papers, tutorials and books, and a special 
mailing list (special interest group) at:


http://desolve.r-forge.r-project.org

Have fun!

Thomas Petzoldt



Succes, Franklin
-


Thanks, Franklin :)






Franklin Bretschneider
Dept of Biology
Utrecht University
brets...@xs4all.nl

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] PDF form Rstudio

[Thomas Petzoldt]

Yes, you are right. Sweave depends on Latex too, so its no workaround in 
this case.


Hope it helps, thpe

Am 26.02.2016 um 06:54 schrieb Ulrik Stervbo:

My understanding is that Sweave also depends on LaTeX to generate pdfs, so
I am not sure Sweave is the solution.



Just follow the advice given in the error message:


No TeX installation detected (TeX is required to create PDF output). You
should install a recommended TeX distribution for your platform:

   Windows: MiKTeX (Complete) - http://miktex.org/2.9/setup
   (NOTE: Be sure to download the Complete rather than Basic

installation)


   Mac OS X: TexLive 2013 (Full) - http://tug.org/mactex/
   (NOTE: Download with Safari rather than Chrome _strongly_ recommended)

   Linux: Use system package manager



__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] Problems with the deSolve package

[Thomas Petzoldt]

Hi,

to diagnose your problem, I tried the model with your original 
parameters, using a Runge-Kutta fixed step integrator and smaller time 
steps. You may try to make them even smaller.


rk4 has no warranted accuracy, so it is less reliable than "lsoda" etc. 
However, it can be useful for debugging, because it runs constantly 
through, even if NA or NaN's occur. We see that the states grow 
drastically to extremely high (or even negative) values.


This points indeed to unrealistic parameter values, or to a 
model-misspecification, which means that necessary feedback mechanisms 
are missing, so that i12 can grow, even if nothing is left (i1, i2) from 
which it can grow from.


Regards, Thomas



init <- c(i1=10, i2=10, i12=0)
parameters <- c(alpha1=0.7, alpha2=0.5, beta1=0.5, beta2=0.3,
  gamma1=0.5, gamma2=0.5, delta=0.5, N=100)
times <- seq(0, 10, by=.01)

simul <- as.data.frame(ode(y = init, times = times, func = system,
  parms = parameters, method="rk4"))

head(simul, 200)
timei1i2  i12
1   0.00  1.00e+01  1.00e+01 0.00e+00
2   0.01  1.685572e+01  1.433028e+01 6.451953e-01
3   0.02  2.567008e+01  1.875164e+01 4.259949e+00
4   0.03  3.251911e+01  2.091655e+01 4.215409e+01
5   0.04 -1.270280e+00 -1.988876e+00 1.581789e+02
6   0.05 -8.108876e+02 -4.992572e+02 2.631133e+04
7   0.06 -4.166359e+73 -2.975971e+73 6.427467e+98
8   0.07   NaN   NaN  NaN
9   0.08   NaN   NaN      NaN

...



--
Dr. Thomas Petzoldt
Technische Universitaet Dresden
Faculty of Environmental Sciences
Institute of Hydrobiology
01062 Dresden, Germany

E-Mail: thomas.petzo...@tu-dresden.de
http://tu-dresden.de/Members/thomas.petzoldt

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] R-help mailing list activity

[Thomas Petzoldt]

Hi,

from my perspective as R user and package maintainer I would consider 
the normalization of the r-help mailing list a good sign. r-help is 
still a good place for general questions, while more specific 
discussions moved to the r-sig-... mailing lists.


Maybe a slight reduction can also be a motivation for more people to 
step in again answering questions.


Thomas

Am 23.01.2016 um 13:28 schrieb Jean-Luc Dupouey:

Dear members,

Not a technical question:

The number of threads in this mailing list, following a long period of
increase, has been regularly and strongly decreasing since 2010, passing
from more than 40K threads to less than 11K threads last year. The trend
is similar for most of the "ancient" mailing lists of the R-project. I
cannot imagine the total number of R-related inquiries on the Internet
decreased. It means that contributors have gone elsewhere. Indeed, in
the meantime, the number of R posts on stackoverflow passed from 2K to
100K between 2009 and 2015. Thus my question: what are the
specificities, the plus and minus of the R-project mailing lists, in
comparison with other lists, and especially in comparison with
stackoverflow? A lot of threads are duplicated on both lists, which
seems to me a little bit counterproductive.

I hope it is the wright place to ask this question. Thanks in advance,

Jean-Luc Dupouey

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] deSolve package

[Thomas Petzoldt]

Dear Andre,

some people would be happy to assist you, but your code fragment is 
incomplete, so diagnosing your problem was impossible. Please read the 
posting guide and provide a minimum reproducible example.


Furthermore I wonder why you use dede and not ode, and why sum() with 
empty parenthesis (that returns always zero).


Thomas Petzoldt

Am 10/15/2015 um 8:17 PM schrieb Andre Jackson:

I have the following differential equations and return list:

dCgd.dt = -kad*y[1]-kgd*y[1] # PK model equation gut d

   dCld.dt= kad*y[1]-rhyd-rmetd   #pk model equation liver d

   dCgl.dt = -kal*y2[1]-kgl*y2[1] # PK model equation gut l

   dCll.dt= kal*y2[1]-rhyl-rmetl   #pk model equation liver l



 # PK model equation

   return(list(c(dCgd.dt,dCld.dt,dCgl.dt,dCll.dt),c(massbalance=sum(

}



For the DDE solve I have :

out = dede(y=y0,times=times,func=model.LIDR,parms=parms)



This gives me the following error:

Error in func(time, state, parms, ...) : object 'p' not found



Can someone explain this error and its remedy?

Andre Jackson

jacksonan1...@gmail.com

"You must be the change you wish the world to be"




[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] fitting grid-based models

[Thomas Petzoldt]

Hi,

I have no example at hand, but the usual way could be indeed to
calculate a time series of an adequate statistic (e.g. spatial
statistic) from both, the observed and the simulated data and then to
apply standard model fitting.

Thomas

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] fitting grid-based models

[Thomas Petzoldt]

On 17.06.2014 10:37, Javier Rodríguez Pérez wrote:

Hi Thomas, Thank you very much for your help. I will try to calculate
some summary statistics and fit observed and simulated data using a
time series (as your examples). With observed data, I would use two
values (e.g. 0 and 10), because I do not have intermediate values.
Thanks! Javier


Fitting models is often more an art than a science, but there are also
approaches to overcome the dilemma, cf.:

Piou, C. Berger, U., Gririmm, V. (2009) Proposing an information
criterion for individual-based models developed in a pattern-oriented
modelling framework. Ecological modelling 220(17) 1957-1967. DOI:
10.1016/j.ecolmodel.2009.05.003

... and related work about individual-based models and pattern oriented
modelling, e.g. the book of Grimm and Railsback.

Thomas

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Difficulty coding time-forced functions in deSolve

[Thomas Petzoldt]

Hi,

does the following help you?

Thomas


require(deSolve)

SEIR - function(t, x, p) {

  if (t  7)
xlag - x
  else
xlag - lagvalue(t - 7)

  V - ifelse(xlag[3]/sum(xlag)  0.01, 0.25, 0)

  ## uncomment the following for printing some internal information
  #cat(t=, t,  -- )
  #cat(xlag,  -- , V, \n)

  N - sum(x)
  with(as.list(c(x,p)),{
dS-  b*N - d*S - beta*S*I/N - V*S
dE- -d*E+beta*S*I/N - epsilon*E - V*E
dI- -d*I + epsilon*E - gamma*I - mu*I - V*I
dR- -d*R + gamma*I + V*S + V*E + V*I
list(c(dS, dE, dI, dR), N = unname(N), V = unname(V))
  })
}

num_years  - 1
time_limit - num_years*365.00

b   - 1/(10.0*365)
d   - b
beta- 0.48
epsilon - 1/4
gamma   - 1/4
mu  - -log(1-0.25)*gamma
parms   - c(b=b, d=d, beta=beta, epsilon=epsilon, gamma=gamma, mu=mu)

xstart -c(S=999, E=0, I=1, R=0)
times   - seq(0.0, time_limit, 1.0)

out -  dede(xstart, times, SEIR, parms, rtol=1e-8, atol=1e-8)

plot(out)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Difficulty coding time-forced functions in deSolve

[Thomas Petzoldt]

On 4/2/2014 5:51 AM, Aimee Kopolow wrote:


Any pointers as to how I can code a function that relies on solutions
from previous time steps?


Such a system would be called a delay differential equation (DDE). It 
can be solved with the dede function, see ?dede for details.


However if you want to model something like this:

 Explicitly:
 I want to introduce vaccination 7 days after the proportion of I2/N2
 reaches 0.01.


Than this is called root finding, that can be combined with events, 
see example EVENTS triggered by a root function in ?events.


More can be found in the papers listed at:

http://desolve.r-forge.r-project.org

... or you may consider to ask the R-sig-dynamic models mailing list:

https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

Hope it helps

Thomas Petzoldt

Dr. Thomas Petzoldt
Technische Universitaet Dresden
Faculty of Environmental Sciences
Institute of Hydrobiology
01062 Dresden, Germany

http://tu-dresden.de/Members/thomas.petzoldt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] deSolve, unresolved namespace error -- solved

[Thomas Petzoldt]

We have been able to reproduce the reported issue on another Linux
system: Fedora 19, and the solution was quite simple:

The deSolve package must always to be loaded *before* loading the shared
library of the compiled model.

Thomas

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] deSolve, unresolved namespace error

[Thomas Petzoldt]

On 11/6/2013 6:50 PM, Adam Clark wrote: Addendum:

unloading and reloading deSolve.so does indeed fix the problem:

library.dynam.unload(deSolve, libpath=paste(.libPaths()[1],
//deSolve, sep=))
library.dynam(deSolve, package=deSolve, lib.loc=.libPaths()[1])

However, this is a little clunky, and seems like overkill. Does
anybody have an idea for a more elegant workaround?



Adam,

the solver lock is used for the ODEPACK solvers to prevent simultaneous
(i.e. nested) calls of the solver within R session, because the ODEPACK
algorithms use some global variables. The RK solvers do not use
global variables, so they have no lock.

If you use the ode solvers in the intended way, an internal call of:

 on.exit(.C(unlock_solver))

should always unlock the solver, even if the function is exited due to
an error. However, your example may point to another problem in the
code, for which we need a minimal reproducible example.


As a first workaround, you may try to call:

.C(unlock_solver)

... but calling an internal .C package function outside a package is
generally not recommended, so you should definitely try to find the real
cause of your problem.

Thomas



PS: please move this thread to either R-Devel or (preferred) to the
specialised mailing list:

mailto:r-sig-dynamic-mod...@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] force pdf()

[Thomas Petzoldt]

Hi,

this is not a problem of R it is a problem of the pdf viewer.

Solution: just use an alternative pdf viewer like gsview or (even 
simpler) SumatraPDF.


Hope it helps

Thomas P.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] ODE solver

[Thomas Petzoldt]

Hi,

just press ESC and the solver should stop.

BTW, warnings like this are usually a sign of unrealistic parameters,
problematic model equations or inappropriate tolerance settings.

See also:

http://cran.r-project.org/web/packages/deSolve/vignettes/deSolve.pdf#46

ThPe


On 5/2/2013 10:36 AM, Tjun Kiat Teo wrote:

I am trying to use the  package ode and periodically it will come up with
this error message

  Warning..Internal T (=R1) and H (=R2) are
  such that in the machine, T + H = T on the next step
  (H = step size). Solver will continue anyway.

And then the program just take very long to run. Is there anyway to get the
program to terminate when this warning is issued instead of continuing to
run ?


Tjun Kiat


[...]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] lsoda question from deSolve package

[Thomas Petzoldt]

Hi Andreas,

adding a varying amount of something over time is even easier than 
implementing an immediate event. It can be implemented as a forcing, 
see ?forcings and example below.


Hope it helps,

Thomas


require(deSolve)

pkmod - function(t, y, p) {
  inp - forc(t)

  if (t  tin) R - (D/tin) else R - 0
  dy1 - R - (p[k12]+p[k])* y[1] + p[k21]*y[2]   + inp
  dy2 - p[k12]* y[1] - p[k21] *y[2]
  list(c(dy1, dy2), inp=inp)
}

## an example, not your data
intimes - c(0, 2, 5, 10, 24)
input   - c(0, 200, 100, 0, 0)
forc - approxfun(intimes, input, method=constant)

times - seq(0, 24, 0.1)
tin  - 0.5
D- 400
V- 26.3
yini - c(Central = 0, Peripheral = 0)
p- c(k=0.05, k12=0.197118, k21=0.022665, V=26.3)

result - lsoda(yini, times, pkmod, p, rtol = 1e-08, atol = 1e-08)
plot(result)


--
Thomas Petzoldt
Technische Universitaet Dresden
Faculty of Environmental Sciences
Institute of Hydrobiology
01062 Dresden, Germany
http://tu-dresden.de/Members/thomas.petzoldt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] using and event in deSolve

[Thomas Petzoldt]

Hi Jannetta,

as far as I can see, your implementation was almost ;-) correct, except 
that:


init[2] - init[1] + d

should be:

init[2] - init[2] + d



The root function with:

return(init[1]-30)

is correct, because this triggers an event when init[1]-30 crosses the 
zero line, i.e. when init[1] == 30. This can also be seen at the figure 
if we decrease the external step size, see the corrected example below.


Note also, that ; at the end of line is not required (and should be 
avoided), because R is not C.


You can also send such questions to the special interest group:

r-sig-dynamic-mod...@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models


Thanks for using deSolve

Thomas



library(deSolve)

Izhikevich - function(time, init, parms) {
  with(as.list(c(init, parms)), {
dv - 0.04 * v^2 + 5 * v + 140 - u + I
du - a * (b * v - u)
list(c(dv, du))
  })
}

parms - c(a = 0.02, b = 0.2, c = -65, d = 2, I = 10)
#times - seq(from = 1, to = 1000, by = 0.1)
times - seq(from = 1, to = 100, by = 0.005)
init  - c(v = -65, u = 0.2)

root - function(time, init, parms) {
  return(v = init[1] - 30)
}

event - function(time, init, parms) {
  with(as.list(parms), {
init[2] - init[2] + d
init[1] - c
return(init)
  })
}

out - ode(y = init, times = times,
  func = Izhikevich, parms = parms,
  events = list(func = event, root = TRUE),
  rootfun = root)

plot(out)__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] R2OpenBUGS question with differential equations

[Thomas Petzoldt]

Dear Andreas,

I don't understand your question, because I can't see any fundamental 
difference between differential equation models and any other R 
function. So why not rewriting your diffeq. model in the same way like 
the other example?


And, please use package deSolve instead of odesolve because (1) the 
latter is deprecated and (2) the successor package deSolve is much more 
flexible and powerful, and it has much better documentation and 
examples, cf. http://desolve.r-forge.r-project.org


Thomas P.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] stiff delay differential equations

[Thomas Petzoldt]

On 12/5/2012 9:15 AM, Suzen, Mehmet wrote:

Hello List,

Can you recommend me if odeSolve can handle stiff delay differential
equations with discontinuities? Or any other package? Best, -m


Package deSolve (the successor of odesolve) can candle:

- stiff differential equations: yes (solvers: lsoda, lsode, vode, vode,
radau)

- delay differential equations: yes (function dede, that calls lsoda,
lsode, vode, radau, ...)

- discontinuities: yes (events and root finding)

A combination is also possible, but please keep in mind that this can
become difficult by definition, so that it is good modelling practise
to start with a simple system first - and to introduce all the
complications only when necessary.

I am not aware that any other package contains all the requested
features, but if please let me know if there is something new and we'll
add it to the DifferentialEquations task view.

Thomas Petzoldt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Solving a multinomial gompertz partial differential equation in r

[Thomas Petzoldt]

On 12/3/2012 9:12 AM, Suzen, Mehmet wrote:

Hi Brandon,

You can try ReacTran package:

cran.r-project.org/web/packages/ReacTran/vignettes/PDE.pdf

Best,
-m


ReacTran is for managing the tr5ansport in reactive transport models, is 
relies on package deSolve that contains the ODE/PDE solvers, so I would 
recommend to start with this first.


Papers, slides and tutorials with full source code can be found at:

http://desolve.r-forge.r-project.org

... and an overview about packages dealing with differential equations 
can be found in the Taskview:


http://cran.r-project.org/web/views/DifferentialEquations.html


Hope it helps

Thomas Petzoldt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Defining a variable outside of optim or differential equation solver.

[Thomas Petzoldt]

Hi Tjun Kat,

you can define variables outside the ode function, but normally NOT 
state variables, because their values need to be updated by the solver 
during the simulation process.


But, if you want to block this for any debugging purposes and want to 
e.g. fix a derivative to a certain value, even this is possible. Note 
however that this is a very special case and I suspect that you don't 
want this.


Can you please tell, why you want to define states outside? I guess you 
want to emulate a feature that is already available in deSolve, e.g. 
forcings or events. In that case, please have a look into the 
documentation and one of the papers tutorial slides etc. that can be 
found on:


http://desolve.r-forge.r-project.org



Note also that your code contains 3 errors:

1) The call must be function(t, y, p), i.e. with p even if this is 
not required by the model, because ode needs this interface.


2) the closing parenthesis ) of list is missing.

3) dvdpol vs. vdpol

Hope it helps

Thomas Petzoldt


On 7/18/2012 3:59 AM, Tjun Kiat Teo wrote:
 This is applicable to either using optim  or the differential equation
 solver or any similar solver

 Suppose I want to use the differential equation solver and this is my 
code


 d-y[2]

 vdpol-function(t,y)
 {
 list(c(1,
 d,
 3,
 4
)
 }


 stiff-ode(y=rep(0,4),times=c(0,1),func=dvdpol,parms=1)


 The thing is I want d to be composed of one of state variables in the
 differential function vdopl. Can it be done ?

 tjun kiat

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html

 and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] How to Read source code of function as string vector or matrix

[Thomas Petzoldt]

On 7/12/2012 11:00 AM, purushothaman wrote:

Hi,

I need to read source code of function as a string vector or matrix.

if i am using get  method it return as closure type.

example

readcsv-function(filepath)
{
output-read.csv(filepath)
}

how to get as string vector like this
test[1]=readcsv-function(filepath)
test[2]={
test[3]=output-read.csv(filepath)
test[4]=}

Thanks
B.Purushothaman


?readLines

Thomas

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 find best parameter values using deSolve n optim() ?

[Thomas Petzoldt]

On 6/6/2012 3:50 PM, mhimanshu wrote:

Hello Thomas,

This code seems to be fine and its now working well.

I read the about the FME package, but I have one doubt, as in the data set
given in the paper, it showing a nice kinetics of the viral growth, so my
question is what if there is a sudden increase in viral growth after some
interval, say Bimodal growth curve?

How does it fits the bimodal growth curve? I tried with FME but I am not
getting the desired results. May be you can explain me a little, I would be
really grateful to you. :)

Thanks a lot,
Himanshu


Hi and thanks for the feedback,

regarding your problem with a bimodal growth curve I am not completely 
sure what you mean. However, I suspect that failing to fit a bimodal 
behaviour may not be caused by parameter fitting with FME, but instead 
would need extension of the underlying ODE model.


Thomas

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 find best parameter values using deSolve n optim() ?

[Thomas Petzoldt]

On 4/5/2012 6:32 PM, mhimanshu wrote:

Hi Thomas,

Thank you so much for your suggestion.
I tried your code and it is working fine. Now when I change the values of Y
in yobs I am getting so many warnings.

say,
yobs- data.frame(
  time = 0:7,
  Y = c(0.00, 3.40, 4.60 ,5.80, 5.80, 6.00, 6.00 ),
  Z = c(0.1, 0.11, 0.119, 0.128, 0.136, 0.145, 0.153, 0.16)
)


This data set is wrong because time and Z have 8 elements but Y has only 
7. Let's assume the missing Y is 6.0.



So when i fit the model with the same code that you have written, i got the
following warnings:
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
   such that in the machine, T + H = T on the next step
  (H = step size). Solver will continue anyway.
   In above,  R1 =  0.1484502806322D+01   R2 =  0.2264549048113D-16

and I have got so many such warnings.


This means that the step adjustment algorithm was unable to determine an 
appropriate step size (h), possibly because of wrongly negative 
parameter values.



Can you explain me why this is happening??  and Secondly,  I dont understand
why i am getting parameters values in negatives after fitting??


You can restrict the parameter ranges by using the optional arguments 
upper and lower, see ?modFit for details.



Can you
please help me out in this... :)


The reason for this can be manifold, e.g. wrong model specification, 
unrealistic data, inappropriate accuracy settings, inadequate start 
parameters or just unidentifiable parameters, e.g. if some of them 
depend strongly on each other.


The example below showed clearly that the parameters are strongly 
correlated (0.998 or even 1.000 !!!). This means that not all parameters 
can be identified simultaneously because of high interdependency (see my 
comment in my first post), so you should try to reduce the number of 
fitted parameters. Note however, that Ymax needs to be fitted (or 
adjusted) as well.


Fitting nonlinear models can be an art and requires mathematical 
understanding of the applied techniques (differential equations, 
identifiability analysis, numerical optimization), a good understanding 
of the properties of your (differential equation) model -- and some 
feeling and experience.



Thomas

##-

library(deSolve)
library(FME)

## function derivs
derivs - function(time, y, pars) {
  with (as.list(c(y, pars)), {
dy = a1 * Y * (1 - Y/Ymax) - a2 * (Y + 0.001)
dz = a3 * Y - a4 * Z;
return(list(c(dy, dz)))
  })
}

## parameters
pars - c(a1 = 0.9, a2 = 0.7, a3 = 0.06, a4 = 0.02, Ymax=0.8)
#Ymax - c(0.8)
## initial values
y - c(Y = 0.2, Z = 0.1)
## time steps
times - c(seq(0, 10, 0.1))
## numerical solution of ODE
out - ode(y = y, parms = pars, times = times, func = derivs)


yobs - data.frame(
 time = 0:7,
 Y = c(0.00, 3.40, 4.60 ,5.80, 5.80, 6.00, 6.00, 6.00),
 Z = c(0.1, 0.11, 0.119, 0.128, 0.136, 0.145, 0.153, 0.16)
)

modelCost - function(p) {
  out - ode(y = y, func = derivs, parms = p, times = yobs$time)
  return(modCost(out, yobs, weight=mean))
}

## start values for the parameters
startpars - c(a1 = 5, a2 = 0.002, a3 = 0.002, a4 = 0.02, Ymax=7)
plower - c(a1 = 1, a2 = 0.0001, a3 = 0.0001, a4 = 0.0001, Ymax=0.001)
pupper - c(a1 = 10, a2 = 2, a3 = 1, a4 = 1, Ymax=10)


## fit the model; nprint = 1 shows intermediate results
fit - modFit(f = modelCost, p = startpars,
  upper=pupper, lower=plower, control = list(nprint = 1))
summary(fit)

## graphical result
out2 - ode(y = y, parms = startpars, times = times, func = derivs)
out3 - ode(y = y, parms = fit$par, times = times, func = derivs)
plot(out, out2, out3, obs = yobs)

legend(topleft, legend=c(original, startpars, fitted),
  col = 1:3, lty = 1:3)

summary(fit)

#Parameters:
#  Estimate Std. Error t value Pr(|t|)
#a1   5.610e+00  2.913e+03   0.0020.998
#a2   2.000e+00  2.908e+03   0.0010.999
#a3   1.872e-03  3.162e-03   0.5920.566
#a4   1.188e-04  1.224e-01   0.0010.999
#Ymax 8.908e+00  4.615e+03   0.0020.998
#
#Residual standard error: 0.08001 on 11 degrees of freedom
#
#Parameter correlation:
#   a1   a2   a3  a4 Ymax
#a11.0  1.0 -0.09916 -0.1031  1.0
#a21.0  1.0 -0.09917 -0.1032  1.0
#a3   -0.09916 -0.09917  1.0  0.9981 -0.09917
#a4   -0.10315 -0.10315  0.99807  1. -0.10316
#Ymax  1.0  1.0 -0.09917 -0.1032  1.0

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 find best parameter values using deSolve n optim() ?

[Thomas Petzoldt]

Hi Himanshu,

the use of optim is described on its help page.

In addition to this package FME provides additional functionality for 
fitting ODE models. See FME help files, package vignettes and the 
example below for details.


Hope it helps

Thomas Petzoldt



library(deSolve)
library(FME)

## function derivs
derivs - function(time, y, pars) {
  with (as.list(c(y, pars)), {
dy = a1 * Y * (1 - Y/Ymax) - a2 * (Y + 0.001)
dz = a3 * Y - a4 * Z;
return(list(c(dy, dz)))
  })
}

## parameters
pars - c(a1 = 0.9, a2 = 0.7, a3 = 0.06, a4 = 0.02)
Ymax - c(0.8)
## initial values
y - c(Y = 0.2, Z = 0.1)
## time steps
times - c(seq(0, 10, 0.1))
## numerical solution of ODE
out - ode(y = y, parms = pars, times = times, func = derivs)

## example observation data
yobs - data.frame(
  time = 0:7,
  Y = c(0.2, 0.195, 0.19, 0.187, 0.185, 0.183, 0.182, 0.18),
  Z = c(0.1, 0.11, 0.119, 0.128, 0.136, 0.145, 0.153, 0.16)
)

## model cost function, see help file and vignette for details
modelCost - function(p) {
  out - ode(y = y, func = derivs, parms = p, times = yobs$time)
  return(modCost(out, yobs))
}

## start values for the parameters
startpars - c(a1 = 1, a2 = 0.6, a3 = 0.1, a4 = 0.1)

## fit the model; nprint = 1 shows intermediate results
fit - modFit(f = modelCost, p = startpars, control = list(nprint = 1))

## Note the high correlation between a1 and a2, resp a3 and a4
## that can make parameter identification difficult
summary(fit)

## graphical result
out2 - ode(y = y, parms = startpars, times = times, func = derivs)
out3 - ode(y = y, parms = fit$par, times = times, func = derivs)
plot(out, out2, out3, obs = yobs)

legend(topleft, legend=c(original, startpars, fitted),
  col = 1:3, lty = 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.

Re: [R] Logistic population growth and deSolve

[Thomas Petzoldt]

Hi Thomas,

I've been out of office for a time, but maybe you are still waiting ...

As far as I see your model is a correctly implemented ODE (!!!) system 
and I don't understand what you mean with unexpected results.


Would it be possible that your original intention was not a *continuous* 
ODE but a discrete system instead?


In this case, just use the euler method:


## continuous system (ordinary differential equation, ODE)
output = ode(y = state, times = times, func = logGrowth, parms = parameters)
plot(output)

## discrete system (difference equation)
output = ode(y = state, times = times, func = logGrowth, parms = 
parameters, method=euler)

plot(output)


Note also that you can write your problem more elegantly using vector 
and matrix notation.


Thomas Petzoldt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] DLSODA error

[Thomas Petzoldt]

Kristian,

your example is rather complex and somewhat time consuming even on fast 
computers, so it is not easy to track your problem down to its original 
reason.


However, what is clear to me is that this is definitely *not*  not a 
problem of deSolve and most likely even not of optim, it is just a 
problem of your model specification. Your error messages point to a 
numerical problem:


DLSODA-  Warning..Internal T (=R1) and H (=R2) are
  such that in the machine, T + H = T on the next step
 (H = step size). Solver will continue anyway.

In above,  R1 =  0.2630651367072D+01   R2 =  0.2055665829864D-15

This is an internal message of the original lsoda Fortran code telling 
us that it cannot add H=0.2e-15 to to T=0.2e+01. Note that the relative 
numerical resolution of double precision arithmetic of the computer is 
about 1e-16.


To understand what I mean just try other ode methods, e.g.:


outmat - ode(y = state, times = times,
  func = Kristian, parms = parameters, method=adams)

## or:

outmat - ode(y = state, times = times,
  func = Kristian, parms = parameters, method=ode45)


Thomas P.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] precompiled ode with spline input

[Thomas Petzoldt]

Hi Daniel,

thanks for your positive response about using compiled ODE functions. 
Did you use package deSolve?


Regarding the use of splines (as forcing functions ?) we don't have an 
out of the box method yet, but you may consider to approximate the 
splines by linear segments or contribute your own spline code to the 
deSolve project.


The example Forcing_lv.R / Forcing_lv.c shows how to use (linearly) 
interpolated forcing data within compiled code, and the interpolation 
code itself (Initdeforc, updateforc) is found in forcings.c. Depending 
on your particular problem there may be other possibilities too.


I don't want to go more into the details here on r-help, but you may 
consider to post additional questions to the mailing list of the 
respective special interest group


https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

This is also a good idea in general for such questions for getting more 
immediate responses.


Thomas Petzoldt



On 3/15/2011 1:00 PM, Daniel Kaschek wrote:

I tried to use lsode with precompiled C code instead of an R function
defining the derivatives. No problem so far.

However, now I need to implement an ODE that contains spline functions,
i.e. the derivatives at given time points depend on the value of a
spline function at this time point. What is an efficient way to
implement this in precompiled C code?




__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] deSolve: Problem solving ODE including modulo-operator

[Thomas Petzoldt]

Dear Albert2002,

there is no problem with deSolve and, of course, no problem with R's 
modulo operator, but there are at least two errors in your model 
formulation:


1.) The order of the returned derivatives must be exactly the same as 
specified in the state variables. This is documented in the help files 
and mentioned in section Troubleshooting (11.2) in the deSolve 
vignette. You use:


  State - c(Theta = 1 ,  P = 1)
  return(list(c(dP, dTheta)))

but you should use:

  State - c(P = 1, Theta = 1)
  return(list(c(dP, dTheta)))


2.) Your model is **not a differential equation** but a difference 
equation. It is (1) discrete (not continuous) and returns the new state 
(not the derivative). As the name of deSolve suggests, this package is 
primarily for differential equations. Nevertheless, it can be useful for 
difference equations too, if one respects the distinction.


Solution A:
===

Use the new development version of deSolve from

  http://deSolve.r-forge.r-project.org

that has a solver method iteration for this type of models:

out1 - ode(func = standardmap1, y = State,parms = Parameter,
   times = Time, method = iteration)

plot(out1)


Solution B:
===

Rewrite your model so that it returns the 'derivative' and not the new 
state and use method=euler. This works already with recent versions of 
deSolve.


iterations - 100
Parameter - c(k = 0.6)
State - c(P=1, Theta = 1 )
Time - 0:iterations

standardmap2 - function(Time, State, Parameter){
  with(as.list(c(State, Parameter)), {
  P1  - (P + k * sin(Theta)) %% (2 * pi)
  Theta1 - (P + Theta) %% (2 * pi)

return(list(c(P1-P, Theta1-Theta)))
  })
}


out2 - ode(func = standardmap2, y = State, parms = Parameter,
   times = Time, method = euler)

plot(out2)

##

You may also consider to rewrite your problem in matrix form to get the 
full map easily and without using loops.


If you have further questions, consider to subscribe to the dynamic 
models in R mailing list:


r-sig-dynamic-mod...@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

Hope it helps

Thomas Petzoldt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] non-linear fourth-order differential equations

[Thomas Petzoldt]

Hi Yanika,

more information about deSolve (books, papers, tutorials) can be found 
on the deSolve homepage:


http://desolve.r-forge.r-project.org

and, of course, in the package documentations. If you need further help 
from the list, please provide a short reproducible example.


Thomas Petzoldt


--
Dr. Thomas Petzoldt
Limnology and Ecological Modelling

Technische Universitaet Dresden
Fakulty of Forest, Geo and Hydro Sciences
Institute of Hydrobiology
01062 Dresden, Germany
E-Mail: thomas.petzo...@tu-dresden.de
http://tu-dresden.de/Members/thomas.petzoldt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] interpretation of stress in NMDS

[Thomas Petzoldt]

On 7/22/2010 2:21 PM, Graves, Gregory wrote:

Among those users of Primer, stress values greater than 0.3 are
interpreted as questionable.  Using both isoMDS and metaMDS (vegan
package), the stress values returned are much higher using my own data
and using examples provided in R Help.  For example Rstress = 8.3, and
the stressplot r2 = 0.99 indicating (to me) that the ordination is OK.
I am guessing that the stress values across packages are not the same,
and googling about has not returned a satisfactory answer ... thus this
posting.  My concern being that reporting a stress value of 8 for what I
consider a satisfactory result may raise a few Primer-user's eyebrows.


Dear Gregory,

as the help file of ?isoMDS tells us:

stress   The final stress achieved (in percent)

So 8.3 is in reality 0.083.

Thomas Petzoldt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Desolve package: How to pass thousand of parameters to C compiled code?

[Thomas Petzoldt]

On 6/8/2010 2:03 PM, Ben Bolker wrote:

  cmmu_mileat  yahoo.com  writes:



Hi,

I have used DeSolve package for my ODE problem regarding
infectious disease transmission and currently am
trying to pass lots (roughly a thousand) of model parameters
to the C compiled model (I have to use C
compiled code instead of R code purely because of the speed).

I can't go define it one by one as it gonna take ages to finish
and also quite difficult to revise. I have read
the instructions in Writing Code in Compiled Languages, which
is very well written, but still have no
clue on how to proceed.

Please can anyone suggest the best way to go about this, and also the places

where I can find examples of

DeSolve C compiled code.



   There's another example of compiled code in
http://www.jstatsoft.org/v33/i09/paper , but it probably
won't say anything the 'Writing Code' guide doesn't already
cover.  I found the guide pretty complete ...

   I would suggest that you pack your parameters into a single
numeric vector (or several, if that makes more sense) and pass
them that way.
   What kind of infectious disease data do you have that will
allow you to fit a model with thousands of parameters ... ??
(Just curious.)



I agree in using vectors (and matrices) like recommended by Ben. In 
addition, you may consider to use structures and unions instead of 
#define macros:


union my_parms_vec {
  struct {double foo1, foo2, bar[15];};
  double value[18];
};

static my_parms_vec p;

// ...

// This structure has then to be initialized
// as vector 'value' in initmod:

/* initializer  */
void initmod(void (* odeparms)(int *, double *))
{
int N=3;
odeparms(N, p.value);
}

// so that you can use constructs like

p.foo or p.bar[7] in your model (i.e. derivs) function.


With respect of the large number of parameters you may think whether 
they are really parameters and not something like external forcings, 
for which recent versions of deSolve have separate mechanisms to deal 
with. See ?forcings


Hope it helps


Thomas Petzoldt

PS: There is also a dedicated mailing list for discussions like this:
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] finite difference scheme for 2D differential equations

[Thomas Petzoldt]

Hello Subodh Acharya,

I've been away for a field trip for two weeks, and I guess you have 
already found a solution for your problem.


If your question is still open, you may consider to look into package 
deSolve that now provides functions (ode1D, ode2D and ode3D) for solving 
1-3 dimensional equations using the method of lines approach. The 
package comes with many examples and extensive documentation and there 
is also a JSS publication:


http://www.jstatsoft.org/v33/i09

If you need more help, please write to the dedicated mailing list for 
dynamic models:


https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

Hope it helps

Thomas Petzoldt

PS: You may also have a look into package ReacTran.


At 26.04.2010 03:09 Subodh Acharya wrote:

Hello everyone,
I am trying to solve 2D differential equations using finite difference
scheme in R. I have been able to work with the equations with only one
spatial dimensions but I want to extend it to the two dimensional problem.
For example i can simulate one dimensional diffusion using a code like the
following. But I want to write a similar code for,say, a two dimensional
diffusion equation. Any kind of help/advice is highly appreciated.

Here is how I did for the 1D equation (in this case for 1D diffusion)
Equation: dc/dt  = D*d^2c/dx^2


dx = 10
dt = 0.1
D = 15
n = 20
tstep = 200
C = 50
dt = 0.1
dx = 10
conc-matrix(0,tstep,n)
conc[1,1:10] = C
  for (i in 2:tstep) {
for (j in 1:n){
if (j==1){
conc[i,j] = conc[i-1,j] - dt*(D/dx^2)*(conc[i-1,j] - conc[i-1,j+1])
}
conc[i,j] = conc[i,j]
  if(j1  j  n){
conc[i,j] = conc[i-1,j] + dt*(D/dx^2)* (conc[i-1,(j-1)] - 2*conc[i-1,j] +
conc[i-1,(j+1)])
}
if (j==n){
conc[i,n] = conc[i-1,n] + dt*(D/dx^2)* ( conc[i-1,n-1] - conc[i-1,n])
}
conc[i,j] = conc[i,j]
}
}

Now in 2D the equation will be like this

dc/dt  = Dx*d^2c/dx^2 + Dy*d^C/dy^2

So that when I solve it I will get C(x, y, t) at each node of the grid.

Thanks for the help in advance.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 generate Mackey-Glass time series with ddesolve package?

[Thomas Petzoldt]

Hello Mike,

please provide a reproducible example, so that we can see why and how 
ddesolve hangs. In addition you may consider using function dede from 
package deSolve which allows to choose between different solver methods 
for dealing with delay differential equations.


Thomas Petzoldt


Am 03.05.2010 08:20, schrieb Mike Beddo:

I could use some help generating a time series for the Mackey-Glass
equation: dx/dt = 0.2 x(t-tau)/(1 + x(t-tau)^10) - 0.1 x(t) for the
case there tau = 17. I tried the ddesolve package but the dde(...)
function seems to hang, not producing anything. Can someone show me
the R script how to do this?


-  Mike Beddo

[[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] estimating the starting value within a ODE using nls and lsoda

[Thomas Petzoldt]

Hi Dave,

first of all, fitting starting values of a dynamic model the same way 
like its parameters is indeed the usual method. In that case parameters 
*and* some or all initial value(s) of the dynamic model are both in fact 
'parameters' for the statistical model fitting problem.


Fitting a nonlinear model can be easy or problematic or even impossible, 
depending on the data and the model structure. In such cases one speaks 
about identifiability and there are several methods that can help to 
find parameter combinations that can be identified simultaneously.


It is not important whether such a statistical parameter was
originally a 'parameter' or an 'initial value' of the dynamic model,
it simply depends on collinearity of the problem:

http://en.wikipedia.org/wiki/Multicollinearity


Several methods for identifiability analysis are provided in the CRAN 
package FME (Flexible Modelling Environment), which comes with 
extensive documentation (package vignettes as pdf files) and examples 
and there is a recent paper in the Journal of Statistical Software:


http://www.jstatsoft.org/v33/i03


Look for function 'collin' that implements a collinearity index.

In addition, function 'modFit', that is also in this package provides an 
interface to use different optimizers of R in a unique way, namely nls, 
nlminb, optim, nls.lm (from package minpack), and pseudoOptim, a 
pseudo-random search algorithm.


Another hint, because you are still using odesolve (from Woodrow 
Setzer). This package has now a direct and compatible successor 
deSolve (Soetaert, Petzoldt, Setzer), see:


http://www.jstatsoft.org/v33/i09

deSolve is even more stable than the former version and contains many 
many more solvers, not only lsoda and rk4, can handle other types of 
differential equations too and has much more documentation.


Hope it helps!

Thomas Petzoldt


PS: There is also a dedicated mailing list for such questions:
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models



--
Thomas Petzoldt
Technische Universitaet Dresden
Institut fuer Hydrobiologie
01062 Dresden
GERMANY

http://tu-dresden.de/Members/thomas.petzoldt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] External signal in ODE written in C (using deSolve and approx1?)

[Thomas Petzoldt]

Dear Glenn, dear list,

this is just a short notice, that a new version 1.5 of package deSolve 
was released yesterday. It now supports the feature requested below. 
Details are documented in the package vignette Writing Code in Compiled 
Languages that comes with the package and is also available online:


http://cran.r-project.org/web/packages/deSolve/vignettes/compiledCode.pdf

See section 6. Using forcing functions.

Thomas Petzoldt (co-author of deSolve)


Glenn Woodart wrote:

Dear list

The deSolve package allows you to specify the model code in C or Fortran.
Thanks to the excellent vignette this works fine. However I have not yet
managed to use forcing functions in C code.

In pure R code this works very well with approxfun() specified outside the
model:


[...]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] New Special Interest Group: R-SIG-DYNAMIC-MODELS

[Thomas Petzoldt]

Dear useRs and developeRs,

In the last few months we have had an increasing number of requests
about deSolve and dynamic modelling. We have also received several
requests to establish a SIG of our own about such things. We think
this is a good idea  -- here it is!

---
R-sig-dynamic-models --
Special Interest Group for Dynamic Simulation Models in R

https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

mailto:r-sig-dynamic-mod...@r-project.org

A mailing list for discussing the use of R for implementation,
simulation and analysis of dynamic simulation models. The list covers
differential equation models (ODE, DAE, PDE, DDE) as well as other
dynamic systems (e.g. individual-based or agent-based simulations) in
different application areas ranging from natural sciences to economy.
---

The list is intended as a forum for users and developers with
different skills and background and for everybody who uses R for
modelling dynamic systems from basic to advanced, for answering
technical questions, for announcing new books, courses and
conferences, for discussing examples and solutions or for promoting
the development of new packages and tools.

Please feel free to subscribe, to ask questions and to contribute!


Thomas Petzoldt
Woodrow Setzer
Karline Soetaert

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] simulating a model

[Thomas Petzoldt]

Dear Rafael,

first of all, your simulation works, at least in a technical sense, so I 
don't understand what you mean with can't simulate it properly.


Second, your SIR-based model is a quite different from the SIR models I 
know (e.g. http://en.wikipedia.org/wiki/SIR_Model).


The R code, however, seems to be technically correct, if compared with 
your system of equations. To help you solving your problem, we need more 
information, e.g. how the equations where derived, where the parameters 
come from, what is the process behind, and, most important, why do you 
think that the outcome is wrong.


In addition, I guess that your simulation time is too long, compared 
with the speed of the process. Try something like


times = c(from=0, to=1, by=0.01)


Thomas



Rafael Moral wrote:

Dear useRs,

I have written an ecological model, based on the epidemiology SIR model.
I've been trying to simulate it in R.
However, I can't simulate it properly.
Two guesses: my script isn't right; I'm not setting the parameters properly

I have uploaded an image to the model here:
http://img24.imageshack.us/img24/743/imagemutr.jpg

The script I am using is as it follows:

require(simecol)

mod1 - new(odeModel,
  main = function(time, init, parms) {
  x - init
  p - parms
dx1 - p[K] - p[alpha]*x[1]*x[2] - p[gamma]*x[1]
dx2 - x[1]*x[2]*(p[alpha] - p[beta])
dx3 - p[beta]*x[1]*x[2] + p[gamma]*x[1]
list(c(dx1, dx2, dx3))
},
times = c(from=0, to=100, by=0.1),
parms = c(K=100, alpha=0.3, gamma=0.5, beta=0.2),
init = c(S=500, V=100, R=0),
solver = lsoda
)

plot(sim(mod1))

Thanks in advance!
Rafael.


  

[[elided Yahoo spam]]

[[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] object .trPaths not found

[Thomas Petzoldt]

Hi,

using a variable named .trPaths is a feature and a frequent nuisance of 
recent Tinn-R and also explained in the Tinn-R docs.


Solution 1:

In Tinn-R select the menu:

 R -- Configure -- and then Temporarily or Permanent

Solution 2:

Manually edit file R-installation-directory/etc/Rconfig.site

The minimal line to add is:

.trPaths - paste(paste(Sys.getenv('APPDATA'), '\\Tinn-R\\tmp\\', 
sep=''), c('', 'search.txt', 'objects.txt', 'file.r', 'selection.r', 
'block.r', 'lines.r'), sep='')



Hope it helps

Thomas P.


Knut Krueger schrieb:

Uwe Ligges schrieb:




I have to admit that I have no idea what we are talking about here 
(yes, I tend to forget many things these days) - and you have not 
cited the original message, unfortunately (nor have you specifies R 
versions, Tinn-R versions and both OS versions, but just one) ...


Best wishes,
Uwe

Original question
rkevinbur...@charter.net wrote:
I am running an R script with Tinn-R (2.2.0.1) and I get the error 
message

 
Error in source(.trPaths[4], echo = TRUE, max.deparse.length = 150) : 
  object .trPaths not found


Any solutions?

System

R 2.8.1
Tinn-R 2.2.02
windows Xp professional Servicepack 3

This combination works on my Pc without problems
And the code is working on the other system with Tinn-R problems
with the R-Editor without problems
so it seems to be a Tinn_R problem

but maybe anybody has a solution for that.

Regards Knut

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] External signal in ODE written in C (using deSolve and approx1?)

[Thomas Petzoldt]

Glenn Woodart wrote:

Dear Thomas

Thank you for your excellent help. I agree method 3 would be the best one,
and this is the one I have been trying to get working (like you described).
So far my skills with pointers, what is returned etc in pure C code is
limited, so my attempts have failed so far. 


That's why I use R for most cases, even if it is somewhat slower than C 
-- is simply much more comfortable.



If it is possible to get a
working example in the vignette it would be great.


Yes, this will surely come ..., but for the moment, see the ad-hoc 
example on:


http://www.simecol.de/models/c_ode_signal.zip

You will see that in addition to linear interpolation one needs a method 
to transfer the data of the signal to the C side. There are several 
possibilities and I'm not yet sure which one will eventually find its 
way to the package.


Note also that the code is a quick hack without any warranty.


Again, thank you for your brilliant work on implementing solvers for
differential equations in R.

Best wishes
Glenn


Thanks for giving feedback.

Thomas Petzoldt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] coupled ODE population model

[Thomas Petzoldt]

Justin Frank wrote:


I'm fairly new to R, and I'm trying to write out a population model that
satisfies the following;

the system consists of s species, i= 1, 2,...,s
network of interactions between species is specified by a (s x s) real matrix,
C[i,j]

x[i] being the relative population of the ith species (0 = x[i] = 1,
sum(x[i]=1)

the evolution rule being considered is as follows;

xprime[i] = f[i] if x[i]  0 or f[i] = 0

xprime[i] = 0 if x[i] = 0 and f[i]  0

where f[i] = sum(C[i,j]*x[j]) - x[i]*sum(C[k,j]*x[j])


I have a bit of attempted code written out, but are there any tricks or tips
that would condense or make this mess look nicer?

-Justin


Hi Justin,

you may consider using package simecol (http://www.simecol.de). You find 
several examples on the help pages and also some PDF files (so called 
vignettes). If you have many species, you should consider using matrix 
formulations. The example below demonstrates a multi-species 
Lotka-Volterra model with an interaction matrix A.


While the example is not yet your model it may serve as a starting 
point. If you have any further questions, please let me know.



Thomas Petzoldt





library(simecol)
##
## Basic Multi Species Predator Prey Model
##  1) parameters for pairwise single species interaction
##

LVPP - new(odeModel,
  main = function(t, n, parms) {
with(parms, {
  dn - r * n  + n * (A %*% n)
  return(list(c(dn)))
})
  },
  parms = list(
r = c(r1 = 0.1, r2 = 0.1, r3 = -0.1, r4 = -0.1),
## only pairwise interactions:
A = matrix(c(0.0, 0.0, -0.2, 0.0,  # prey 1
 0.0, 0.0, 0.0, -0.1,  # prey 2
 0.2, 0.0, 0.0, 0.0,   # predator 1; eats prey 1
 0.0, 0.1, 0.0, 0.0),  # predator 2; eats prey 2
 nrow = 4, ncol = 4, byrow=TRUE)
  ),
  times = seq(from=0, to=500, by = 0.1),
  init  = c(prey1=1, prey2=1, pred1=2, pred2=2),
  solver = lsoda
)

plot(sim(LVPP))

##
## 2) multi species interactions
##

## make two clones of LVPP
LVPPweak - LVPPstrong - LVPP

## a helper function

## this copies the negative of the upper triangular to the lower
makeSym - function(A, f = -1) {
  ind - lower.tri(A)
  A[ind] - f * t(A)[ind]
  A
}

## weak coupling
A - matrix(c(0.0, 0.0, -0.1, -0.001,   # prey 1
  NA,  0.0, -0.002, -0.2,   # prey 2
  NA,  NA,   0.0,   0.0,# predator 1
  NA,  NA,   NA,0.0),   # predator 2
  nrow = 4, ncol = 4, byrow=TRUE)

parms(LVPPweak)$A - makeSym(A)

## stronger coupling
A - matrix(c(0.0, 0.0, -0.1, -0.05,# prey 1
  NA,  0.0, -0.02, -0.2,# prey 2
  NA,  NA,   0.0,   0.0,# predator 1
  NA,  NA,   NA,0.0),   # predator 2
  nrow = 4, ncol = 4, byrow=TRUE)

parms(LVPPstrong)$A - makeSym(A)


LVPPweak - sim(LVPPweak)
LVPPstrong - sim(LVPPstrong)


plot(LVPPweak)
plot(LVPPstrong)


o - out(LVPPweak)
par(mfrow=c(2,2))
plot(o$prey1, o$pred1)
plot(o$prey2, o$pred2)

plot(o$prey1, o$pred2)
plot(o$prey2, o$pred1)


o - out(LVPPstrong)
par(mfrow=c(2,2))
plot(o$prey1, o$pred1)
plot(o$prey2, o$pred2)

plot(o$prey1, o$pred2)
plot(o$prey2, o$pred1)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] External signal in ODE written in C (using deSolve and approx1?)

[Thomas Petzoldt]

Glenn Woodart wrote:

Dear list

The deSolve package allows you to specify the model code in C or Fortran.
Thanks to the excellent vignette this works fine. However I have not yet
managed to use forcing functions in C code.

In pure R code this works very well with approxfun() specified outside the
model:


[... example deleted, see original post ...]


I would like to do the same thing in C, and my guess is that the approx1
function has to be used in some way. So far did not figure out how. If
anyone has managed to do so, or know how to approach this problem, please
let me know.

Best wishes
Glenn


Hi Glen,

this problem is on our radar for a while, however, we have not found the 
time yet to implement this in a general way. There are, in principal, 
several different possibilities:


1) use a solver with a fixed internal step size, e.g. rk4 and supply the 
data in the appropriate discretization (note also the intermediate 
steps). The upcoming deSolve 1.3 on R-Forge is shortly before to be 
released and has now all solvers in C resp. Fortran.


2) make a back-call from your C model to the R function approx. This is 
the most flexible method, but makes your code considerably slower. See 
Writing R Extensions how to do this.


3) write your own linear interpolation function in C. This is quite 
simple as even R uses the bisection method (in approx1, as you correctly 
identified). So take approx1, make a copy and simplify it to your needs.


Solution (1) needs the least effort in C programming and, in my opinion, 
well behaving ODE systems that require considerable interpolation effort 
for forcing data are one of the very few cases where rk4 may have 
retained its ecological niche. But even in such cases there is still the 
problem of unknown numerical accuracy.


The best method is certainly (3) and we would be very glad to integrate 
such a function (or a documented example) into the deSolve or simecol 
package.


Thomas Petzoldt


--
Dr. Thomas Petzoldt
Technische Universitaet Dresden
Institut fuer Hydrobiologiethomas.petzo...@tu-dresden.de
01062 Dresden  http://tu-dresden.de/hydrobiologie/
GERMANY

http://tu-dresden.de/Members/thomas.petzoldt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] deSolve question

[Thomas Petzoldt]

insun nam wrote:

Dear All,

I like to simulate a physiologically based pharmacokinetics model using R
but am having a problem with the daspk routine.

The same problem has been implemented in Berkeley madonna and Winbugs so
that I know that it is working. However, with daspk it is not, and the
numbers are everywhere!

Please see the following and let me know if I am missing something...

Thanks a lot in advance,
In-Sun


[ ... long example deleted, see original post ...]

Hi In-Sun,

as far as I see the *first* of your problems is your high demand on 
precision. Double precision allows 16 digits, so 1e-10 may lead 
internally to values that are close to the maximum number of digits in 
double precision arithmetics, because you set relative *and* absolute 
tolerance of all states to such a small value.


Reduce atol and rtol to reasonable values (e.g. the default of 1e-6) and 
it should work (it does on my machine). Another possibility would be to 
increase maxsteps, but this will not help if the requested precision is 
higher than would be possible with double precision.


If you try the following:

ODE - as.data.frame(daspk(y = y, times = times, func = Fun_ODE,
parms = pars, atol = 1e-6, rtol = 1e-6))


... then the simulation goes through in a technical sense, however state 
variables reach very high values which, again, are beyond what is 
possible with double precision arithmetics -- is this what you mean with 
the numbers are everywhere?


Are you sure that your model formulation is correct? Please check your 
equations, especially parentheses and signs.


Thomas Petzoldt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] using lsoda() and nls() together

[Thomas Petzoldt]

Hi Benoit,

your problem is not really a problem of lsoda. The reason of the crash 
is a violation of the statistical assumptions of least squares 
regression due to dependency of residual variance on x. Due to this, K1 
is varied over a very large range of values until numeric overflow occurs.


Note that you have an exponentially growing state, so log transformation 
will help:



res - nls(log(foo) ~
  log(func(K1)),start=list(K1=1),data=data.frame(foo=y), trace=TRUE)

summary(res)


You may also consider using packages simecol (on CRAN) or FME (on 
R-Forge) that both support constrained optimization of ode systems.


Hope it helps

Thomas Petzoldt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] self organizing map advice for categorical data

[Thomas Petzoldt]

George Chen schrieb:

Hello,

Could anybody offer any advice about implementing a Kohonen self organizing map 
for categorical data?  Specifically I am wondering if there are any 
pre-existent packages that can deal with categorical data and/or how one would 
compare the input vector of categoricals with the self organizing map nodes.

Thanks in advance.

George Chen



Yes, there is a very nice one with excellent graphics, see article:

Ron Wehrens, Lutgarde M. C. Buydens (2007) Self- and Super-organizing 
Maps in R: The kohonen Package. Journal of Statistical Software	21(5). 
http://www.jstatsoft.org/v21/i05



Hope it helps

Thomas Petzoldt

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] extract the p value of F statistics from the lm class

[Thomas Petzoldt]

Hi,

what about the following:

## some test data
x - 1:10
y - x + rnorm(x)

## model and summary
m - lm(y~x)
sm - summary(m)
sm

# str(sm)
# sm$fstatistic

## and now: the manual case
1 - pf(sm$fstatistic[1], sm$fstatistic[2], sm$fstatistic[3])


Hope it helps, ThPe


tedzzx schrieb:

Dear R users

I have run an regression and want to extract the p value of the F
statistics, but I can find a way to do that.

x-summary(lm(log(RV2)~log(IV.m),data=b))

Call:
lm(formula = log(RV2) ~ log(IV.m), data = b[[11]])

Residuals:
 Min   1Q   Median   3Q  Max 
-0.26511 -0.09718 -0.01326  0.11095  0.29777 


Coefficients:
Estimate Std. Error t value Pr(|t|)
(Intercept)  -0.3059 0.1917  -1.5950.121
log(IV.m) 0.9038 0.1065   8.488 1.38e-09 ***

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 


Residual standard error: 0.1435 on 31 degrees of freedom
Multiple R-squared: 0.6991,	Adjusted R-squared: 0.6894 
F-statistic: 72.04 on 1 and 31 DF,  p-value: 1.379e-09 


names(x)
 [1] call  terms residuals
 [4] coefficients  aliased   sigma
 [7] dfr.squared adj.r.squared

[10] fstatisticcov.unscaled

x$fstatistic
   valuenumdfdendf 
72.04064  1.0 31.0 

But can not find the p value of F statistics. 


Thanks

Ted





__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] R Style Guides?

[Thomas Petzoldt]

Robert Biddle schrieb:

Hi:

I am looking for style guides for larger R programs,
and would like suggestions of web pages or documents to look at.

I am particularly interested in guidelines for how to structure programs,
and such issues as managing scope, passing data around, data hiding,
and so forth.
I do have some understanding of how scope works, and how
object-oriented programming does or might work.
But I am not looking for more explanations of the detail mechanisms.

What I want is documented advice on how to use these mechanisms
for laying out larger programs, with rationale, and preferably with 
examples.


Suggestions, please?

Thanks
Robert Biddle



Hi Robert,

the following book review may help you. Don't be confused by the last 
word in the book title.


http://www.jstatsoft.org/v29/b08

Thomas P.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] ttest in R

[Thomas Petzoldt]

1Rnwb schrieb:

Thanks for the help, but ANOVA will give me a single pvalue, then how i can
make sure which group is showing the significant differences.


Hi,

ANOVA is fine and please have a look on ?TukeyHSD and don't forget to 
consult your statistics textbook about post-hoc testing. If you insist 
in using t.test, don't forget Bonferroni or Holm-correction (in R: 
?p.adjust) or use ?pairwise.t.test !!!


ThPe

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] ttest in R

[Thomas Petzoldt]

1Rnwb wrote:

since the estimation is not done pairwise so i cannot use pairwise.t.test,
how do i apply tukeyHSD


Note correct capitalization: TukeyHSD and follow the examples on the 
help page:


?TukeyHSD

You may also contact a statistics textbook. In addition I send you some 
links off list.


ThPe

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Problems in Recommending R

[Thomas Petzoldt]

Hi,

you are probably right, though I must say that I like *spartanic and 
efficient* homepages and I don't think that the example given by the 
first mail is a good prototype for the R homepage. But, yes, occasional 
face lifting may be adequate.  Anti-aliasing is of course simple, but 
that's probably not the point. (And I know that there are graphics 
experts with a masters in psychology between us.)


So, why not a new Homepage Graphics Competition 2009? There is still 
some time until useR!2009 in Rennes:


http://www2.agrocampus-ouest.fr/math/useR-2009/


Thomas Petzoldt


PS to all useRs: Don't miss to visit the useR!2009 homepage and note the 
upcoming submission deadline 2009-02-27 and the exciting location of the 
conference!




Stavros Macrakis wrote:

A first step that would make the current Web page look much better
would be to anti-alias the demonstration graphic.  The current graphic
makes R graphics seem (falsely!) to be very primitive. I'm afraid I
don't know how to do the anti-aliasing myself.

Replacing the fixed-width, typewriter-style font with something a bit
more elegant might also be good

-s

On Sun, Feb 1, 2009 at 9:52 PM, Ajay ohri ohri2...@gmail.com wrote:

Dear List,
One persistent feedback I am getting to people who are newly introduced to R
( especially in this cost cutting recession)  is -

1) The website looks a bit old.

While the current website does have a lot of hard work behind it, should n't
a world class statistics package have a better website instead.

You can check out www.knime.org which is an open source software , and free,
and supports R---and notice the change in perception .

Best Regards,

Ajay Ohri


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] 1-Pearson's R Distance

[Thomas Petzoldt]

Rodrigo Aluizio wrote:

Hi again List,

Well this time I?m writing for a friend (really J). He needs to create a
distance matrix based on an abundance matrix using the 1-Pearson?s R index.
Well I told him to look at the proxy package, but there is only Pearson
Index. He needs it to perform a clustering. Well, as soon as he told me
there proxy only had the Pearson index I thought: ?He could just do
something like


NewObject-1-PearsonMatrixObject?


But I didn?t tell him that because I?m not sure it?s the same thing, and it
probably will generate a strange cluster with the braches ends distant from
the base?

He told me that Statistica 7.0 has this Index, but he doesn?t own it. So? Is
there a way on R to do this correctly?

 


Thanks for the attention.


The help file of dist has an example how to do this:

?dist


## Use correlations between variables as distance
dd - as.dist((1 - cor(USJudgeRatings))/2)


Note the division by 2! Correlations can be between -1 ... +1 and 
negative distances make no sense.


Another approach is r^2, but this has, of course, a different 
interpretation.


ThPe

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] LaTeX and R-scripts/-results

[Thomas Petzoldt]

Oliver Bandel wrote:

Hello,

at some places I read about good interaction of
LaTeX and R.

Can you give me a starting point, where I can find
information about it?

Are there special LaTeX-packages for the support,
or does R have packages for support of LaTeX?
Or will an external Code-Generator be used?

TIA,
   Oliver



Hi Oliver,

you are right, LaTeX and R are perfect companions. Look for Sweave(*). 
You find an introduction of Fritz Leisch in R-News 2002, Vol 2/3:


http://cran.r-project.org/doc/Rnews/Rnews_2002-3.pdf

and an entire homepage about it:

http://www.statistik.lmu.de/~leisch/Sweave/

HTH Thomas P.

(*) Many thanks to Friedrich Leisch for this great peace of software!

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Count days of current year

[Thomas Petzoldt]

Felipe Carrillo wrote:

Hi:
Is there a function that counts the number of days of any given or current year based on the date?  


See ?Date and ?strptime and try the following:

y - as.numeric(format(as.Date(2008-11-25), %Y))

as.numeric(format(as.Date(paste(y, 12, 31, sep=-)), %j))


HTH Thomas P.



Felipe D. Carrillo  
Supervisory Fishery Biologist  
Department of the Interior  
US Fish  Wildlife Service  
California, USA


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Rendering Dendrograms

[Thomas Petzoldt]

Hi Mark,

similar questions were at least two times during the last weeks, see

http://tolstoy.newcastle.edu.au/R/e5/help/08/11/6722.html

http://tolstoy.newcastle.edu.au/R/e5/help/08/11/6736.html

or

http://tolstoy.newcastle.edu.au/R/e5/help/08/11/7790.html

or search the archives yourself:

http://tolstoy.newcastle.edu.au/R/

search for cluster graph or plot.hclust


Thomas P.


Mark Cook wrote:




Hello all
I've been using the hclust and as.dendrogram objects for hierarchical
clustering. The problem I have is that my sample set is now so large (circa
500 points) that it isn't possible to view the leaf nodes.

I'd like to be able to zoom in on specific areas of the graph by 
selecting a

region with the mouse.  I've deduced from trial and error that the xlim and
ylim parameters don't work with these plots and the only method I've had 
any

success with is to cut the tree at a specific point and display the lower
region.  However, the cut function returns a list of all subtrees below the
specified height and I have no way of selecting a particular subtree (that
is, cutting on the x-axis).

The identify and rect.hclust functions seem to do this OK  and  in fact 
these

partially solve my problem but I want to use getGraphicsEvent() if possible
since this also allows interaction with the keyboard (so that I can use the
backspace key for example to return to the previous plot). I'm fairly 
new to

R so may be wrong but it appears that the plots are single threaded and can
respond to either a set of graphics events or an identify loop but not 
both)


It has been suggested that I look at  Graphviz since this has more
comprehensive mapping facilities. Is this likely to help? Is there for
example a method for converting an R hclust object to an Graphviz ICLUST or
some alternative method for describing an hclust in the Graphviz .DOT
notation?

Or can anybody suggest an alternative approach ??

Regards

Mark








__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] lsoda warning too much accuracy requested

[Thomas Petzoldt]

A short addendum, resulting from an off-list discussion:

The reason why Colleen's code failed was raising a negative base to a 
fractional exponent in the third state equation for certain sets of 
parameters, esp. fractional values of beta.


Old versions of odesolve broke down, and recent versions of deSolve (or 
the deprecated odesolve) simply returned NaN for the third state. In 
order to track such problems down, it is helpful to call the system 
separately, e.g.:


 model(0, start, parms)
[[1]]
HBA N
0.3702103 1.2229436   NaN

and then isolate the problem in the respective state equation.

Thomas P.



Thomas Petzoldt wrote:

Hi Colleen,

this error was not uncommon and is usually a sign of a numerically 
problematic or wrongly implemented model. Please use package deSolve, 
the successor of odesolve, that is more robust and has also a bunch 
alternative solvers for difficult cases.


I tested your code with deSolve (on R 2.8.0, Windows) and it runs 
without problems.


Thomas Petzoldt


BTW: your system worked also with odesolve, so my question: which 
versions (R, odesolve) and operating system are you using?


Colleen Carlson wrote:

Dear list -

 

Does anyone have any ideas / comments about why I am receiving the 
following

warning when I run lsoda:

 

1:  lsoda--  at t (=r1), too much accuracy requested in: lsoda(start, 
times,

model, parms)
2:   for precision of machine..  see tolsf (=r2) in: lsoda(start, 
times,

model, parms)

 


I have tried changing both rtol and atol but without success.  I saw the
thread in the R-archive of 11 June 2004 but this has not helped me.

 


I have built the model in stages and the problem only occurs when the
exponent beta in the third DE is anything other than 0 or 1.  If beta 
= 0 or

1 then the solver gives me perfectly justifiable results.  Just changing
beta to 0.9 or similar causes the problem.

 


I am still new to R so I am unsure if it is my programming or my
understanding of the way lsoda works.

 


Any comments or input would be welcome.

Many thanks

Colleen

___

 


My code is:

 


library(odesolve)

SI - 80

model - function(t, x, parms) {

H  - x[1]

BA - x[2]

N  - x[3]

with(as.list(parms), {

dHdt - (b/c)*(((a**c)*((H)**(1-c))-H))

dBAdt - -(BA*b)*(c0+(c1*SI)-log(BA))/(log(1-((H/a)**c)))

dNdt -  N*alpha*(((log(1-((H/a)**c)))/b)**beta) - (gamma*BA)

list(c(dHdt, dBAdt, dNdt))

})

}

times - seq(0, 40, 1)

 


parms - c(a=(SI*1.258621)-1.32759, b=0.1, c=0.4, c0=4.6012, c1=0.013597,
alpha=0.0005, beta=0.5,  gamma=0.01)

start - c(H=0.1, BA=0.1, N=600)

 


out - as.data.frame(lsoda(start, times, model, parms))


[[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] lsoda warning too much accuracy requested

[Thomas Petzoldt]

Hi Colleen,

this error was not uncommon and is usually a sign of a numerically 
problematic or wrongly implemented model. Please use package deSolve, 
the successor of odesolve, that is more robust and has also a bunch 
alternative solvers for difficult cases.


I tested your code with deSolve (on R 2.8.0, Windows) and it runs 
without problems.


Thomas Petzoldt


BTW: your system worked also with odesolve, so my question: which 
versions (R, odesolve) and operating system are you using?


Colleen Carlson wrote:

Dear list -

 


Does anyone have any ideas / comments about why I am receiving the following
warning when I run lsoda:

 


1:  lsoda--  at t (=r1), too much accuracy requested in: lsoda(start, times,
model, parms) 


2:   for precision of machine..  see tolsf (=r2) in: lsoda(start, times,
model, parms)

 


I have tried changing both rtol and atol but without success.  I saw the
thread in the R-archive of 11 June 2004 but this has not helped me.

 


I have built the model in stages and the problem only occurs when the
exponent beta in the third DE is anything other than 0 or 1.  If beta = 0 or
1 then the solver gives me perfectly justifiable results.  Just changing
beta to 0.9 or similar causes the problem.

 


I am still new to R so I am unsure if it is my programming or my
understanding of the way lsoda works.

 


Any comments or input would be welcome.

Many thanks

Colleen

___

 


My code is:

 


library(odesolve)

SI - 80

model - function(t, x, parms) {

H  - x[1]

BA - x[2]

N  - x[3]

with(as.list(parms), {

dHdt - (b/c)*(((a**c)*((H)**(1-c))-H))

dBAdt - -(BA*b)*(c0+(c1*SI)-log(BA))/(log(1-((H/a)**c)))

dNdt -  N*alpha*(((log(1-((H/a)**c)))/b)**beta) - (gamma*BA)

list(c(dHdt, dBAdt, dNdt))

})

}

times - seq(0, 40, 1)

 


parms - c(a=(SI*1.258621)-1.32759, b=0.1, c=0.4, c0=4.6012, c1=0.013597,
alpha=0.0005, beta=0.5,  gamma=0.01)

start - c(H=0.1, BA=0.1, N=600)

 


out - as.data.frame(lsoda(start, times, model, parms))


[[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.



--
Thomas Petzoldt
Technische Universitaet Dresden
Institut fuer Hydrobiologie[EMAIL PROTECTED]
01062 Dresden  http://tu-dresden.de/hydrobiologie/
GERMANY

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Fitting a sine wave using solver

[Thomas Petzoldt]

Ben Zuckerberg schrieb:

Greetings,

I have several sets of oscillation data and would like to estimate the 
parameters of a sine function to each set (and hopefully automate 
this).  A colleague provided an excel sheet that uses solver to minimize 
the RSS after fitting the sine function to each data set, but this 
cumbersome and difficult to automate.  Is there a method in R for 
fitting a given sine function to a supplied data using maximum 
likelihood estimation (or minimizing the RSS).  Thanks in advance.




Hi Ben,

why not using FFT? See one of the examples below.

HTH Thomas P.


### Some periodic data 
n   - 360
ord - 180
x   - 0:(n-1)
y - sin((50 + x) * 2 * pi / n) + rnorm(n) * 0.1
plot(x, y,xlim=c(-10,370))

### example 1 -
## Fast Fourier Transform
pf - fft(y)

## Plot highest possible order
pf[(n/2):n] - 0
yy - 2*Re(fft(pf, inverse=TRUE)/n) - Re(pf[1])/n
lines(x, yy, col=gray, lwd=1)

### example 2 -
## Fast Fourier Transform
pf - fft(y)

## another, lower order
pf[4:n] - 0
yy - 2*Re(fft(pf, inverse=TRUE)/n) - Re(pf[1])/n
lines(x, yy, col=blue, lwd=2)


### example 3 -
## Fast Fourier Transform
pf - fft(y)

## synthesize it the traditional way
n  - length(y)
a0 - Re(pf[1])/n
pf - pf[-1]
a  -  2*Re(pf)/n
b  - -2*Im(pf)/n

harmonic - function(x, a0, a, b, n, ord) {
  k - (2 * pi * x/n) %*% t(1:ord)
  y -a0 + cos(k) %*% a[1:ord] +
   sin(k) %*% b[1:ord]
  y
}

lines(x, harmonic(x, a0, a, b, n, ord=2), col=red, lty=2, lwd=2)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] skeleton for package containing C code?

[Thomas Petzoldt]

Rainer M Krug wrote:

Hi

I looked in the writing R extensions and also into package.skeleton,
but could not find anything:

I am looking for a skeleton or simple example on how to create a
package which contains C++ code.

I only want to implement one function, so nothing complicated. The
idea is to use the package to make it easier to install the function
on several Linux computers - no windows required at the moment.

Thanks,

Rainer


Hi Rainer,

you can download a minimal R package with C source from here:

http://www.simecol.de/pkg/cpackage_0.1.tar.gz

You can also include C++ files in src/ if you provide plain C interfaces 
as described in Writing R Extensions.


Hope that helps ...


Thomas P.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Cluster Graph

[Thomas Petzoldt]

Jared Chapman schrieb:

I am an uber nubie with r but I have a question that I hope someone
can help me with.

I have a data set with 1200+ data points that I want to put into a
cluster graph. the problem is that when the cluster graph is generated
there are too many data points to view the labels, they are just a
jumbled mess at the bottom.

Is there a way to stretch the graph out so that the data points are
visible or view the data groupings in a way that make them readable?

I looked in the R introduction and I could not find a forum to search
for this answer. I am sure that it is something simple.


you are right. A similar question was answered a few days before. See:

http://tolstoy.newcastle.edu.au/R/e5/help/08/11/6722.html

HTH

Thomas P.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Changing the position of the origin

[Thomas Petzoldt]

Plantky wrote:

Hi all,

Can anyone tell me how I can make 0,0 start at the top left hand
corner of a graph, instead of the typical lower left hand corner? I've
tried to plot with axes=F and then putting on the axes later, but I
want the points to correspond to the axes.

Thanks,
Kang Min


Does this what you want?

x -  runif(10)
y - -runif(10)

plot(x, y, xlim=c(0, max(x)), ylim=c(min(y), 0), axes=FALSE)
axis(2)
axis(3)
box()

Thomas P.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Arms Race

[Thomas Petzoldt]

dennis campos wrote:
 hey can anybody help me? i have to simulate the richardson Arms race 
model on R.. for my simulation class...


Hi Dennis,

this list is not intended for solving your homework, however, the 
following may help you to go one step ahead.


You should have a look into the help pages of packages deSolve and/or 
simecol. Please find below the (IMHO rather elegant) simecol 
implementation but beware, your professor might expect the basic 
deSolve style, so be prepared about discussing both. And, don't forget 
to read the background.


## deSolve #
library(deSolve)
? ode

## simecol #

library(simecol)

## open the directory with R sourcecode examples
browseURL(paste(system.file(package=simecol), /examples, sep=))

 and modify the example lv.R


Thomas P.

##
library(simecol)

armsrace - new(odeModel,
  main = function (time, init, parms) {
x - init[1]
y - init[2]
with(as.list(parms), {
  dx - a * y - m * x + r
  dy - b * x - n * y + s
  list(c(dx, dy))
})
  },
  parms  = c(a=2, m=5, r=5, b=2, n=3,  s=5),
  times  = c(from=0, to=10, by=0.5),
  init   = c(x=0.5, y=1),
  solver = lsoda
)

ar1 - ar2 - ar3 - ar4 - armsrace

parms(ar2) - c(a=2, m=1, r=3, b=2, n=2, s=3)
parms(ar3) - c(a=1, m=4, r=-1, b=1, n=1, s=-2)
parms(ar4) - c(a=2, m=1, r=-3, b=2, n=1, s=-3)

ar1 - sim(ar1); plot(ar1)
ar2 - sim(ar2); plot(ar2)
ar3 - sim(ar3); plot(ar3)
ar4 - sim(ar4); plot(ar4)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Agent-based social simulation and R

[Thomas Petzoldt]

Hi Tom,

you may have a look at the CRAN package simecol, that has some examples
how to implement different types of dynamic models in R (differential
equations, grid models, individual based models).

Individual-based models (IBMs) are a model family used in ecology, which
are in its essence almost the same as ABMs in other areas.

See http://simecol.r-forge.r-project.org for the package, examples,
pdf's, and in particular the useR!2008 slides.

Thomas Petzoldt


Tom Backer Johnsen wrote:
Do anyone know anything about the use of R for agent-based social 
simulation?  It should be possible, and would be convenient for the 
simple reason that there are several nice packages containing useful 
stuff for SNA (Social Network Analysis).  Information about packages, 
web sites, experienced persons in the field, etc. would be very welcome.


Tom
++
| Tom Backer Johnsen, Psychometrics Unit,  Faculty of Psychology |
| University of Bergen, Christies gt. 12, N-5015 Bergen,  NORWAY |
| Tel : +47-5558-9185Fax : +47-5558-9879 |
| Email : [EMAIL PROTECTED]URL : http://www.galton.uib.no/ |
++

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] plot.hclust with lots of objects

[Thomas Petzoldt]

Hi Paul,

one possibility: write the tree to a wide pdf file and then zoom and 
scroll using Adobe Acrobat or another PDF reader. Here is a tree with 
1000 objects:


x - matrix(rnorm(3000), nrow=1000)
hc - hclust(dist(x))
pdf(tree.pdf, width=150)
plot(hc)
dev.off()


Thomas P.


paul murima wrote:
 Dear all,

 The default plotting method for hclust trees looks just fine for few
 objects like in the example dataset. But when it comes to many
 features (eg some 1000 and more - I'm trying to visualize clustered
 microarray data) it renders a tree, that one cannot inspect, because
 of overlapping text and lines. My question is, is there a way or a
 plotting parameter for plotting a tree which is wide enough to have
 all leaves separated and readable labels even for that many objects?
 This would produce a very big image, so I think scrolling is
 essential.

 I believe the answer is simple, but I'm unable to figure it out.

 Thanks in advance

 --
 BEST

 Paul Murima

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html

 and provide commented, minimal, self-contained, reproducible code.


--
Thomas Petzoldt
Technische Universitaet Dresden
Institut fuer Hydrobiologie[EMAIL PROTECTED]
01062 Dresden  http://tu-dresden.de/hydrobiologie/
GERMANY

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] lines, ecdf and colors

[Thomas Petzoldt]

eric lee wrote:

Hi.  I'm trying to plot two ecdf's on the same graph using two
different colors.  I can plot using the same color, but it doesn't
work when I change colors?  Any suggestions?  Thanks in advance for
your help.

x - c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46)
y - c(1.15, 0.88, 0.90, 0.74, 1.21)
plot(ecdf(x))
# it works without col='blue', but doesn't with it
lines(ecdf(y), col = blue)


Does the following what you expect:

x - c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46)
y - c(1.15, 0.88, 0.90, 0.74, 1.21)
plot(ecdf(x))
lines(ecdf(y), col.points  = blue, col.hor = blue)


See ?plot.ecdf and especially ?plot.stepfun for details.

HTH Thomas P.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] getting the p-value from lm as a list object

[Thomas Petzoldt]

eric lee wrote:

Hi,

I'm trying to get the p-value from the 'lm' regression function as a list
object.  For example, I can get r^2 from the following code by entering
summary(fm)$r.squared.  Is there a way to get the p-value?  If not, is there
a function where I can enter the f-value and degrees of freedom to get the
p-value?  Thanks.

x - c(1,2,3,4,5,6,7,8,9,10)
y - c(1,2,3,4,4,5,6,8,1,9)

fm - lm(y ~ x)
str(summary(fm))


What about the following (taken from  stats:::print.summary.lm):


x - c(1,2,3,4,5,6,7,8,9,10)
y - c(1,2,3,4,4,5,6,8,1,9)

fm - lm(y ~ x)

summary(fm) # for comparison only

sfm - summary(fm)
pf(sfm$fstatistic[1], sfm$fstatistic[2], sfm$fstatistic[3],
  lower.tail = FALSE)



Thomas P.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Help needed with Waterfall plot

[Thomas Petzoldt]

Philip Twumasi-Ankrah schrieb:

Hi friends,
I need suggestions/directions on how to producing a waterfall plot for present extend of change in tumour size for a set of respondents in a study.  Example of use of waterfall plot is in the following slides presented at ASCO 2007 by Axel Grothey. Link is 


http://media.asco.org/player/default.aspx?LectureID=AG265conferenceFolder=GI2007SessionFolder=Posterslideonly=yesTrackID=N929LectureTitle=Waterfall%20plots%20provide%20detailed%20information%20on%20magnitude%20of%20response%20to%20conventional%20chemotherapy%20in%20colorectal%20cancer%3a%20Lessons%20learned%20from%20N9741.Key=vm_45_3_26_265SpeakerName=%3b%20Presenter%3a%20Axel%20Grothey%2c%20MDmediaURL=%2fmediaServerName=media.asco.orgmax=12ext=jpguseASX=falseplaytype=playtype=playtype=,

The link is pretty long but it takes you right to the presentation.


Hi Phillip,

is this a waterfall plot:

## a few data
x - 0:99
y - sort(rnorm(100), decreasing=TRUE)

# the plot
plot(y, type=n)
polygon(c(min(x), x, max(x), 0), c(0, y, 0, 0), col=green)

Thomas P.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] pretty formatting of lists

[Thomas Petzoldt]

Hello,

is there already a function in any R package which does
source code formatting of deparsed lists?

Let's create the following list:

x - list(a = round(rnorm(3), 2),
   b = round(rnorm(3), 2))

xx -c(aa = round(rnorm(30)), f = function(a) a + b, list(x, x))


Now, I want deparse it in a way that yields something like:


list(
   aa = c(0.25, 0.18, 0.84, -1.25, 0.09, -0.99, 1.64,
 1.42, -1.29, -0.14, -1.07, -1.05, 0.98, 0.33,
 1.76, -1.66, 0.96, -0.21, -1.29, 0.78, -0.4,
 -1.63, -0.78, -1.05, 1.27, -1.44, 0.12, -0.4,
 0.02, 1.03),
   f = function (a) a + b,
   list(
 a = c(2.22, 0.36, 0.74),
 b = c(0.46, 0.41, 1.46)
   ),
   list(
 a = c(2.22, 0.36, 0.74),
 b = c(0.46, 0.41, 1.46)
   )
)


Thanks a lot!

Thomas P.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] pretty formatting of lists

[Thomas Petzoldt]

jim holtman wrote:
 Is this basically what you want?

 xx -list(aa = round(rnorm(30),2), f = function(a) a + b, a=x, b=x)
 xx
 $aa
  [1]  1.34 -0.21 -0.18 -0.10  0.71 -0.07 -0.04 -0.68 -0.32  0.06 -0.59
  0.53 -1.52  0.31 -1.54 -0.30
 [17] -0.53 -0.65 -0.06 -1.91  1.18 -1.66 -0.46 -1.12 -0.75  2.09  0.02
 -1.29 -1.64  0.45

 $f
 function(a) a + b

 $a
 $a$a
 [1] -0.06 -0.16 -1.47

 $a$b
 [1] -0.48  0.42  1.36


 $b
 $b$a
 [1] -0.06 -0.16 -1.47

 $b$b
 [1] -0.48  0.42  1.36



No, I want a neat source code representation, not simply the printed
content -- technically something like the following, but more reader
friendly:

 cat(deparse(xx), fill=60)

structure(list(aa = c(0.25, 0.18, 0.84, -1.25, 0.09, -0.99, 1.64,
1.42, -1.29, -0.14, -1.07, -1.05, 0.98, 0.33, 1.76, -1.66, 0.96,
-0.21, -1.29, 0.78, -0.4, -1.63, -0.78, -1.05, 1.27, -1.44, 0.12,
-0.4, 0.02, 1.03), f = function (a)
a + b, structure(list(a = c(2.22, 0.36, 0.74), b = c(0.46, 0.41,
1.46)), .Names = c(a, b)), structure(list(a = c(2.22, 0.36,
0.74), b = c(0.46, 0.41, 1.46)), .Names = c(a, b))), .Names =
c(aa,  f, , ))


I wonder whether a function for some kind of pretty deparsing already
exists before I start to write my own solution.

Thomas P.


 
 On Sun, Mar 16, 2008 at 8:47 AM, Thomas Petzoldt [EMAIL PROTECTED] wrote:
 Hello,

 is there already a function in any R package which does
 source code formatting of deparsed lists?

 Let's create the following list:

 x - list(a = round(rnorm(3), 2),
   b = round(rnorm(3), 2))

 xx -c(aa = round(rnorm(30)), f = function(a) a + b, list(x, x))


 Now, I want deparse it in a way that yields something like:


 list(
   aa = c(0.25, 0.18, 0.84, -1.25, 0.09, -0.99, 1.64,
 1.42, -1.29, -0.14, -1.07, -1.05, 0.98, 0.33,
 1.76, -1.66, 0.96, -0.21, -1.29, 0.78, -0.4,
 -1.63, -0.78, -1.05, 1.27, -1.44, 0.12, -0.4,
 0.02, 1.03),
   f = function (a) a + b,
   list(
 a = c(2.22, 0.36, 0.74),
 b = c(0.46, 0.41, 1.46)
   ),
   list(
 a = c(2.22, 0.36, 0.74),
 b = c(0.46, 0.41, 1.46)
   )
 )


 Thanks a lot!

 Thomas P.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] pdf() device uses fonts to represent points - data alteration?

[Thomas Petzoldt]

Hello,

I had the same problem. When opening PDFs with a recent developer 
version of inkscape all circles were replaced by the letter q, see a 
screenshoot of the imported figure:

http://www.simecol.de/figs/R_pdf_inkscape.png

I spent at least two hours trying different development versions of 
inkscape, different versions of R, reading docs, trying different 
machines, installing fonts etc., finally giving up. Now, the two 
workarounds of Paul Murrell indeed solved the problem for *me*. Thank 
you. Here are example results of workarounds (i) and (ii):

http://www.simecol.de/figs/R_ps_pdf_inkscape.png

or

http://www.simecol.de/figs/R_pdftrans_inkscape.png

One problem remains. I wanted to demonstrate post-editing PDFs with 
inkscape to motivate students to use R for graphics even if they dont 
want to become programmers. However, double conversion (via 
postscript) or the magics of transparency and opaqueness are not yet the 
way to increase the trust of point-and-click users to R. Maybe this is a 
topic for r-devel?

Thomas


-- 
Thomas Petzoldt
Technische Universitaet Dresden
Institut fuer Hydrobiologie
01062 Dresden
GERMANY

http://tu-dresden.de/Members/thomas.petzoldt?set_language=en

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] pdf() device uses fonts to represent points - data alteration?

[Thomas Petzoldt]

Hi,

this was my former suggestion but as the RSvgDevice doc says:

This driver currently does not have any font metric information, so the 
use of plotmath is not supported.

Maybe, the future (or already recent) alternative will be:

library(Cairo)

that can among others also create svg and pdf. While I had some problems 
with a former version, it worked now out of the box and was able to do 
both, algebra (plotmath) and geometry (circles) ;-) in high quality. 
Let's see if there is any other drawback.

Thomas P.




jiho wrote:

 On 2007-November-01  , at 20:18 , Thomas Petzoldt wrote:
 I had the same problem. When opening PDFs with a recent developer 
 version of inkscape all circles were replaced by the letter q, see a 
 screenshoot of the imported figure:

 http://www.simecol.de/figs/R_pdf_inkscape.png

 I spent at least two hours trying different development versions of 
 inkscape, different versions of R, reading docs, trying different 
 machines, installing fonts etc., finally giving up. Now, the two 
 workarounds of Paul Murrell indeed solved the problem for *me*. Thank 
 you. 

[...]

 
 By the way, depending on what OS you are, you may find an entirely SVG 
 workflow more suitable:
 R with RSVvgDevice package -1- SVG figures -2- Inkscape -3- 
 whatever you like (SVG, PNG, PDF...)
 This gives all transparency, fonts etc to Inkscape so it is fine on this 
 side. The only problem with this workflow for me is that many of my 
 plots stay between stage 1 and 2 and I like to be able to view them 
 quickly. I would need a quick SVG viewer but there are none on OS X. If 
 you are on Linux, many documents viewers (eog, evince, gThumbs) can 
 display SVGs so you would be all set.
 
 JiHO
 ---
 http://jo.irisson.free.fr/
 
 


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.