[R] Apply rmarkdown::render() outside the RStudio don't find pandoc

2014-08-28 Thread walmes .
Hello,

I want to render a Rmd file using rmarkdown::render() function but I'm used
to edit files with Emacs and I open a linux terminal session to compile
this file with knitr::knit2html()

## On the linux terminal.
$ echo "require(knitr); knit2html('teste01.Rmd')" | R --vanilla

I want to use rmarkdown::render instead. When I tried this, with a R
session open in the linux terminal, I get the message

## In a R session inside the linux terminal
> rmarkdown::render("teste01.Rmd")
Error: pandoc version 1.12.3 or higher is required and was not found.

So I open the the teste01.Rmd file in RStudio and clicked on the "knit
HTML" button and it works. The output is bellow.

## Output in the RStudio after clicking on the "Knit HTML" button.

processing file: teste01.Rmd
  |.|  20%
  ordinary text without R code

  |..   |  40%
label: unnamed-chunk-1
  |...  |  60%
  ordinary text without R code

  | |  80%
label: unnamed-chunk-2 (with options)
List of 1
 $ echo: logi FALSE

  |.| 100%
  ordinary text without R code


output file: teste01.knit.md

*/usr/lib/rstudio/bin/pandoc/pandoc* teste01.utf8.md --to html --from
markdown+autolink_bare_uris+ascii_identifiers+tex_math_single_backslash-implicit_figures
--output teste01.html --smart --email-obfuscation none --standalone
--section-divs --table-of-contents --toc-depth 3 --template
/home/walmes/R/x86_64-pc-linux-gnu-library/3.1/rmarkdown/rmd/h/default.html
--css /home/walmes/Dropbox/ridiculas/markdown/ridiculas.css --variable
theme:cerulean --include-in-header
/tmp/RtmpkgSjDp/rmarkdown-str4aee5c225085.html --mathjax --variable
mathjax-url:
https://c328740.ssl.cf1.rackcdn.com/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML
--no-highlight --variable highlightjs=teste01_files/highlight --variable
highlightjs-theme=textmate

Output created: teste01.html

Its possible to see that RStudio call the pandoc that is at
*/usr/lib/rstudio/bin/pandoc/* and not the old version of pandoc in my
system.

I try overcome this by creating an alias
*pandoc='/usr/lib/rstudio/bin/pandoc/pandoc'*, bus this doesn't work.
Someone has any idea to solve this? I want to use rmarkdown::render()
outside the RStudio.

Thanks.
Walmes.

======
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
skype: walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] Apply rmarkdown::render() outside the RStudio don't find pandoc

2014-08-28 Thread walmes .
Yihui,

Thank you so much. I read the PANDOC.md as you suggested and this solved my
issue. I've created the soft-links and then I run in linux terminal

$ echo "rmarkdown::render('teste01.Rmd')" | R --vanilla

R version 3.1.1 (2014-07-10) -- "Sock it to Me"
Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
... skip ...

> rmarkdown::render('teste01.Rmd')

processing file: teste01.Rmd
  |.|  20%

label: unnamed-chunk-2 (with options)
List of 1
 $ echo: logi FALSE

  |.| 100%
  ordinary text without R code


output file: teste01.knit.md

/usr/local/bin/pandoc teste01.utf8.md --to html --from
markdown+autolink_bare_uris+ascii_identifiers+tex_math_single_backslash-implicit_figures
--output teste01.html --smart --email-obfuscation none --standalone
--section-divs --table-of-contents --toc-depth 3 --template
/home/walmes/R/x86_64-pc-linux-gnu-library/3.1/rmarkdown/rmd/h/default.html
--css /home/walmes/Dropbox/ridiculas/markdown/ridiculas.css --variable
theme:bootstrap --include-in-header
/tmp/RtmpW1cQMQ/rmarkdown-str14a07b76054a.html --mathjax --variable
mathjax-url:
https://c328740.ssl.cf1.rackcdn.com/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML
--no-highlight --variable highlightjs=teste01_files/highlight

Output created: teste01.html

Everything works! Thanks you.

W.


​

[[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] Unexpected behaviour of plyr::ddply

2014-09-16 Thread walmes .
Hello R users,

I'm writing a brief tutorial of getting statistical measures by splitting
according strata and over columns. When I used plyr::ddply I got and
unexpected result, with NA/NaN for non existing cells. Below is a minimal
reproducible code with the result that I got. For comparison, the result of
aggregate is showed. Why this behaviour? What I can do to avoid it?

> require(plyr)
>
> hab <-
+ read.table("http://www.leg.ufpr.br/~walmes/data/ipea_habitacao.csv";,
+header=TRUE, sep=",", stringsAsFactors=FALSE, quote="",
+encoding="utf-8")
>
> hab <- hab[,-ncol(hab)]
> names(hab) <- c("sig", "cod", "mun", "agua", "ener", "tel", "carro",
+ "comp", "tot")
> hab <- transform(hab, sig=factor(sig))
> hab$siz <- cut(hab$tot, breaks=c(-Inf, 5000, Inf),
+labels=c("P","G"))
> str(hab, ve.len=1)
'data.frame':5596 obs. of  10 variables:
 $ sig  : Factor w/ 27 levels "AC","AL","AM",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ cod  : int  1200013 1200054 1200104 1200138 1200179 1200203 1200252
1200302 1200328 1200336 ...
 $ mun  : chr  "Acrelândia" "Assis Brasil" "Brasiléia" "Bujari" ...
 $ agua : num  21.5 27.4 26.9 17.3 13.1 ...
 $ ener : num  56.2 65.3 55.9 43.9 35.9 ...
 $ tel  : num  8.85 26.71 22.73 12.28 9.19 ...
 $ carro: num  9.3 6.03 7.47 6.49 5.73 ...
 $ comp : num  0.947 1.637 1.857 0.127 0.088 ...
 $ tot  : int  1878 807 4114 1365 1267 14368 2807 5268 740 2308 ...
 $ siz  : Factor w/ 2 levels "P","G": 1 1 1 1 1 2 1 2 1 1 ...
>
> xtabs(~siz+sig, hab)
   sig
siz  AC  AL  AM  AP  BA  CE  DF  ES  GO  MA  MG  MS  MT  PA  PB  PE  PI
PR  RJ
  P  17  72  47  13 266 103   0  43 187 152 671  52 100  80 195  97 199
310  27
  G   5  29  15   3 149  81   1  34  55  65 182  25  26  63  28  88  22
89  64
   sig
siz  RN  RO  RR  RS  SC  SE  SP  TO
  P 147  34  14 355 236  56 396 129
  G  19  18   1 112  57  19 249  10
>
> x <- ddply(hab, ~sig+siz,
+colwise(.fun=mean,
+.cols=~agua+ener+tel+carro+comp, na.rm=TRUE))
> head(x)
  sig  siz agua ener   telcarro  comp
1  ACP 19.30229 51.30535 12.857118 5.395824 0.7028235
2  ACG 28.39300 65.95740 26.322800 8.942000 2.3806000
3  ALP 42.14337 81.20935  4.801500 7.069958 0.5332639
4  ALG 54.03966 87.47834  9.771724 9.428586 1.3583793
5  *AL   NaN  NaN   NaN  NaN   NaN*
6  AMP 25.66202 61.12709  5.749596 1.980362 0.7629362
>
> x <- ddply(hab, ~sig+siz,
+colwise(.fun=sum,
+.cols=~agua+ener+tel+carro+comp, na.rm=TRUE))
> head(x)
  sig  siz agua ener tel   carro   comp
1  ACP  328.139  872.191 218.571  91.729 11.948
2  ACG  141.965  329.787 131.614  44.710 11.903
3  ALP 3034.323 5847.073 345.708 509.037 38.395
4  ALG 1567.150 2536.872 283.380 273.429 39.393
5  *AL 0.0000.000   0.000   0.000  0.000*
6  AMP 1206.115 2872.973 270.231  93.077 35.858
>
> y <- aggregate(as.matrix(hab[,4:8])~sig+siz, data=hab, FUN=mean,
+na.rm=TRUE)
> head(y)
  sig siz agua ener   tel carro  comp
1  AC   P 19.30229 51.30535 12.857118  5.395824 0.7028235
2  AL   P 42.14337 81.20935  4.801500  7.069958 0.5332639
3  AM   P 25.66202 61.12709  5.749596  1.980362 0.7629362
4  AP   P 37.20362 82.04515 14.929154  5.775923 0.7922308
5  BA   P 41.23200 66.45516  4.524045 10.387203 0.8404624
6  CE   P 36.78599 78.20176  6.339990  7.768981 0.6446990
>
> sessionInfo()
R version 3.1.1 (2014-07-10)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8   LC_NUMERIC=C
 [3] LC_TIME=pt_BR.UTF-8LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=pt_BR.UTF-8LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=pt_BR.UTF-8   LC_NAME=C
 [9] LC_ADDRESS=C   LC_TELEPHONE=C
[11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] plyr_1.8.1

loaded via a namespace (and not attached):
[1] compiler_3.1.1 Rcpp_0.11.2tools_3.1.1
>

Thanks in advance.
Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
skype: walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] Confidence intervals nls

2010-02-15 Thread Walmes Zeviani

The problem was because attr(predict(model, list(DOY=days))) didn't provide
any gradient, it returned NULL. To correct this we need to specify the
gradient by using deriv3() function. Then, when you use predict(model,
list(DOY=days)) will be provided the gradient at this new experimental
points.

# toy data ---
da <- data.frame(DOY=seq(1,120,l=30))
da$CET <- 9.5+(-6.5*sin(((2*pi)/365)*(da$DOY+65)))+rnorm(da$DOY,0,0.2)
plot(CET~DOY, data=da)

# to get the gradient 
ff <- deriv3(~a+(b*sin(((2*pi)/365)*(DOY+t))),
 c("a","b","t"), function(a, b, t, DOY) NULL)

# model fit --
model <- nls(CET~ff(a,b,t,DOY), data=da, start=list(a=9.5, b=-6.5, t=65))

# prediction -
days <- seq(0,365,1)
str(predict(model, list(DOY=days)))
attr(predict(model, list(DOY=days)),"gradient")

se.fit <- sqrt(apply(attr(predict(model, list(DOY=days)),"gradient"), 1,
 function(x) sum(vcov(model)*outer(x,x

# plot ---
matplot(days, predict(model,list(DOY = days))+
outer(se.fit, qnorm(c(.5, .025,.975))),
type="l", col=c(1,2,2), lty=c(1,2,2))

#-----


Sincerely.
Walmes Zeviani.

-
..oooO
......
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/Confidence-intervals-nls-tp1556487p1556702.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] HELP on Non-Linera Mixed-Effect model

2010-02-16 Thread Walmes Zeviani

I never had seen spline() inside some nls()/nlme() function. Are you sure
that this fit is possible? Splines makes a lot or successive local
polynomial fits, so they need a lot of parameters. I don't think that is
possible a good spline fit with only four parameters. In this case you could
use loess() or locfit() in locfit library.

Walmes - Brasil.

-
..oooO
..
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/HELP-on-Non-Linera-Mixed-Effect-model-tp1556953p1557173.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Triangular filled contour plot

2010-02-16 Thread Walmes Zeviani

Johannes,

Some months ago it was posted on R-help a list of packages that handle with
ternaryplots, altough none of them can handle surfaceplots, just
scatterplots, on a triangular area. The packages were

plot.acomp in compositions
tri in cwhmisc.cwhtool
triax in plotrix
ternary in StatDA
ternaryplot in vcd
ternaryplot in Zelig 

See http://n4.nabble.com/triangular-plot-td893434.html#a893435

Due my needs I implemented something a little bit useful (at least to me) to
make a prediction surface triangle plot using Lattice package. My
reproducible code is the following:

paint <- data.frame(mono=c(17.5, 10, 15, 25, 5, 5, 11.25, 5, 18.13, 8.13,
25, 15, 10, 5),
cross=c(32.5, 40, 25, 25, 25, 32.5, 32.5, 40, 28.75,
28.75, 25, 25, 40, 25),
resin=c(50, 50, 60, 50, 70, 62.5, 56.25, 55, 53.13,
63.13, 50, 60, 50, 70),
hardness=c(29, 26, 17, 28, 35, 31, 21, 20, 29, 25, 19,
14, 30, 23))

pseudo <- with(paint, data.frame(mono=(mono-5)/(25-5),
resin=(resin-50)/(70-50)))
pseudo$cross <- with(pseudo, 1-mono-resin)
pseudo$hardness <- paint$hardness

m1 <- lm(hardness~(mono+cross+resin)^2-mono, data=pseudo)

trian <- expand.grid(base=seq(0,1,l=100*2), high=seq(0,sin(pi/3),l=87*2))
trian <- subset(trian, (base*sin(pi/3)*2)>high)
trian <- subset(trian, ((1-base)*sin(pi/3)*2)>high)

new2 <- data.frame(cross=trian$high*2/sqrt(3))
new2$resin <- trian$base-trian$high/sqrt(3)
new2$mono <- 1-new2$resin-new2$cross

trian$yhat <- predict(m1, newdata=new2)

grade.trellis <- function(from=0.2, to=0.8, step=0.2, col=1, lty=2,
lwd=0.5){
  x1 <- seq(from, to, step)
  x2 <- x1/2
  y2 <- x1*sqrt(3)/2
  x3 <- (1-x1)*0.5+x1
  y3 <- sqrt(3)/2-x1*sqrt(3)/2
  panel.segments(x1, 0, x2, y2, col=col, lty=lty, lwd=lwd)
  panel.text(x1, 0, label=x1, pos=1)
  panel.segments(x1, 0, x3, y3, col=col, lty=lty, lwd=lwd)
  panel.text(x2, y2, label=rev(x1), pos=2)
  panel.segments(x2, y2, 1-x2, y2, col=col, lty=lty, lwd=lwd)
  panel.text(x3, y3, label=rev(x1), pos=4)
}

levelplot(yhat~base*high, trian, aspect="iso", xlim=c(-0.1,1.1),
ylim=c(-0.1,0.96),
  xlab=NULL, ylab=NULL, contour=TRUE,
  par.settings=list(axis.line=list(col=NA), axis.text=list(col=NA)),
  panel=function(..., at, contour=TRUE, labels=NULL){
panel.levelplot(..., at=at, contour=contour, #labels=labels,
lty=3, lwd=0.5, col=1)
  })
trellis.focus("panel", 1, 1, highlight=FALSE)
panel.segments(c(0,0,0.5), c(0,0,sqrt(3)/2), c(1,1/2,1), c(0,sqrt(3)/2,0))
grade.trellis()
panel.text(0, 0, label="mono", pos=2)
panel.text(1/2, sqrt(3)/2, label="cross", pos=3)
panel.text(1, 0, label="resin", pos=4)
trellis.unfocus()

http://n4.nabble.com/file/n1557735/one.png 

I think (and hope :-) ) Deepayan Sarkar will implement this on Lattice soon
because R doesn't have this kind of graphic yet. My code isn't perfect so I
expect suggestions.

Walmes Zeviani, Brasil.

-
..ooo0
...
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/Triangular-filled-contour-plot-tp1557386p1557735.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Normality in split-plot design

2010-02-22 Thread Walmes Zeviani

Silvano,

To test normality of the erros you can use

m0 <- aov(...)
shapiro.test(residuals(m0))

or another test like Kolmogorov-Smirnov (ks.test), Jarque-Bera
(jarque.test{moments}, jarque.bera.test{tseries}) among others.
Fell free to join our R brazilian group on yahoo groups named R_STAT. You
will welcome.

Bests.
Walmes Zeviani. Lavras, MG, Brasil.

-
..ooo0
...
..()... 0ooo...      Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/Normality-in-split-plot-design-tp1564486p1564631.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] help with lattice boxplots...

2010-03-01 Thread Walmes Zeviani

In complement to Dallazuanna's solution, use box.umbrella=list() inside
par.settings() to change the default color and line type specification:

bwplot(y~x, data=ex, pch="|",
   par.settings=list(box.rectangle=list(col="black"),
 box.umbrella=list(lty=1, col="black"))
   )


Sincerely.
Walmes Zeviani. Lavras - MG, Brasil.

-
..ooo0
...
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/help-with-lattice-boxplots-tp1573781p1573896.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Correct nested design for GLM

2010-03-03 Thread Walmes Zeviani

I guess your model specification is incoherent. It's seems to me that you
incorporated some kind of lmer() style with "(0 + Micro/Trt/Year)" term.
Consider the toy data and code below. Nesting model is m2. Models m1 and m3
are the same so showing that terms in brackets don't are used.

da <- expand.grid(A=factor(1:3), B=factor(1:3))[rep(1:9, 3),]
da$y <- rnbinom(da$A, size=10, prob=0.5)

m1 <- glm.nb(y~A*B, data=da)
anova(m1)

m2 <- glm.nb(y~A/B, data=da)
anova(m2)

m3 <- glm.nb(y~A*B+(A/B), data=da)
anova(m3)

Bests.
Walmes Zeviani.

-
..ooo0
...
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/Correct-nested-design-for-GLM-tp1576038p1576483.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] confidence intervals for non-linear regression

2010-03-14 Thread Walmes Zeviani

In the paste I had provided one possible solution to this, see

http://n4.nabble.com/Confidence-intervals-nls-td1556487.html#a1556702

Walmes.

-
..ooo0
...
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/confidence-intervals-for-non-linear-regression-tp1592521p1592552.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] testing parallelism of does-response curves using package "drc"

2010-03-17 Thread Walmes Zeviani

Using just nls() you can set the complete model and the alternative model in
this way

da <- expand.grid(x=20:50, trat=c(0.9,1))
da$y <- 10*da$trat*da$x/(12+da$x)+rnorm(da$x,0,0.1)
plot(y~x, da)
da$trat <- as.factor(da$trat)

n0 <- nls(y~A[trat]*x/(B[trat]+x), data=da,
  start=list(A=c(9,10), B=c(12,12)))
n1 <- nls(y~A*x/(B[trat]+x), data=da,
  start=list(A=c(9.5), B=c(12,12)))
anova(n1,n0)

Look at gnls() function in the nlme package for a easier way to specify the
model.

Walmes Zeviani.

-
..ooo0
...
..()... 0ooo...      Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/testing-parallelism-of-does-response-curves-using-nls-tp1591356p1596256.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Expected pairwise.student.t and TukeyHSD behavior?

2010-03-24 Thread Walmes Zeviani

This is because pairwise.student.t and TukeyHSD use differents estimates for
error term. TukeyHSD use the parameters covariance matrix and t use sample
variance. In this case a treatment with only one value doesn't have sample
variance but has a estimated standard error in a covariance matrix. The
lines below will illustrate the point:

# toy data
da <- data.frame(A=factor(rep(c("a","b","c"), c(1,5,7
da$y <- rnorm(da$A)

# sample variance
tapply(da$y, da$A, var)
pairwise.t.test(da$y, da$A, p.adj="none") # problem!

# model and covariance matriz
m0 <- aov(y~A, da)
vcov(m0)
TukeyHSD(m0) # no problem!

Walmes.



-
..ooo0
...
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/Expected-pairwise-student-t-and-TukeyHSD-behavior-tp1690100p1690136.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] lattice: how to add points to the plot generated by levelplot()?

2010-03-31 Thread Walmes Zeviani

You can use

levelplot(runif(100)~rep(1:10, each=10)+rep(1:10, times=10))
trellis.focus("panel", 1, 1, highlight=FALSE)
 lpoints(runif(100,0,10), runif(100,0,10), pch=2, col=2, cex=2)
trellis.unfocus()

# or

levelplot(runif(100)~rep(1:10, each=10)+rep(1:10, times=10),
  panel=function(...){
panel.levelplot(...)
grid.points(runif(100,0,10), runif(100,0,10), pch=2)
  }
  ) 

Walmes

-
..ooo0
...
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/lattice-how-to-add-points-to-the-plot-generated-by-levelplot-tp1747260p1747499.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Line Graph Labels

2010-04-01 Thread Walmes Zeviani

You can adapt the following example

da <- data.frame(y=rnorm(50), x=1:50)
plot(y~x, data=da)
abline(h=c(-2,2), lty=3)
with(da, text(x[abs(y)>2], y[abs(y)>2], label=x[abs(y)>2], pos=2))

Walmes.

-
..ooo0
...
..()... 0ooo...      Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/Line-Graph-Labels-tp1748218p1748431.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Line Graph Labels

2010-04-01 Thread Walmes Zeviani

You can set the number of extreme points to be labeled instead of define a
cutoff. Look:

da <- data.frame(y=rnorm(50), x=1:50)
plot(y~x, data=da)
abline(h=c(-2,2), lty=3)
with(da, text(x[abs(y)>2], y[abs(y)>2], label=x[abs(y)>2], pos=2))

da <- da[order(da$y),]
plot(y~x, data=da)
# five small and big numbers
num <- 5
with(da, points(x[c(1:num, nrow(da):(nrow(da)-num))],
y[c(1:num, nrow(da):(nrow(da)-num))], pch=3, col=2))
with(da, text(x[c(1:num, nrow(da):(nrow(da)-num))],
  y[c(1:num, nrow(da):(nrow(da)-num))],
  label=x[c(1:5,nrow(da):(nrow(da)-5))], pos=rep(c(1,3),
each=num)))

Hope that helps.
Walmes.

-
..ooo0
...
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/Line-Graph-Labels-tp1748218p1748558.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] mouse-clicking on xy-plot

2010-04-02 Thread Walmes Zeviani

You can use identify() to obtain coordinates from plotted points but if you
want any coordinates you could use locator():

> plot(1:10)
> loc <- locator(n=3)
> str(loc)
List of 2
 $ x: num [1:3] 2.3 5.4 8.29
 $ y: num [1:3] 6.15 8.33 2.6
> points(loc$x, loc$y, col=2)
> 

Walmes.

-
..ooo0
...
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/mouse-clicking-on-xy-plot-tp1749562p1749586.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] predict.lm with NAs

2010-04-15 Thread Walmes Zeviani

You can use predict() by specifying a complete data.frame() for prediction to
the argument newdata=. Look:

da <- expand.grid(x1=LETTERS[1:4], x2=1:9)
da$y <- rnorm(da$x1)
da$y[sample(length(da$y), 5)] <- NA

m0 <- lm(y~x1+x2, data=da)
predict(m0) # NA not predicted
predict(m0, newdata=da) # NA predicted

Sincerely.
Walmes.

-
..ooo0
...
..()... 0ooo...      Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/predict-lm-with-NAs-tp1840661p1886457.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] using nls for gamma distribution (a,b,d)

2010-04-15 Thread Walmes Zeviani

I never had seen someone using a density function for a mean function in
nonlinear regression. The convergence problem observed can be due the model
derivatives respect to the parameters. How could computationally be
d.gamma/d.alpha? digamma function? Does R understand this?

If your phenomena exhibit a gamma density form you could use a mean function
like:

fun <- function(x, alpha, beta, theta){
  x*(alpha+beta*x)^(-1/theta)
}

da <- data.frame(x=1:20)
da$y <- fun(da$x, 1, 2, 0.9)+rnorm(da$x,0,0.001)

plot(y~x, da)
curve(fun(x, 1,2,0.9), 1, 100, add=TRUE)

n0 <- nls(y~fun(x, alpha, beta, theta), data=da,
  start=list(alpha=1, beta=2, theta=0.9))

Or you can use a gamma density without the integration equal one. So, you
can replace the gamma(alpha) for B term:

gammafun <- function(x, alpha, B, b){
  (1/B)*alpha^b*x^(b-1)*exp(-x/alpha)
}

curve(gammafun(x, 5, 1, 2), 0, 20)

da <- data.frame(x=1:20)
da$z <- gammafun(da$x, 5, 1, 2)+rnorm(da$x,0,1)

plot(z~x, da)
curve(gammafun(x, 5, 1, 2), add=TRUE)

n1 <- nls(z~gammafun(x, alpha, B, b), data=da,
  start=list(alpha=5, B=1, b=2))

Sincerely.
Walmes.

-
..ooo0
...
..()... 0ooo...      Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/using-nls-for-gamma-distribution-a-b-d-tp1899967p1909193.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] basic table statistics

2010-04-24 Thread Walmes Zeviani

You can use fBasics package

sample(c("A","B","C"),size=20,replace = T) -> type
rnorm(20) -> value
data.frame(ty=type,val=value) -> test

require(fBasics)

nam <- rownames(basicStats(test$val))
result <- do.call("cbind", with(test, tapply(val, ty, basicStats)))
rownames(result) <- nam
result

Bests.

-
..ooo0
.......
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://r.789695.n4.nabble.com/basic-table-statistics-tp2062829p2063832.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] NLS "Singular Gradient" Error

2010-04-28 Thread Walmes Zeviani

I don't know any about you research but I agree with Gabor: the model is
identifiable.

No*(1-exp(a*(b*Ne-T))) can be reparametrized to No*(1-exp(C*NeD)), where
C=a*b and D=a*T. This reduces the model to 3 parameters and can see shown
that is a reperametrization of the SSasympOff() defined in R (with default
start values). "No" is a parameter too, so you need provide start values in
start list. Look:

# indentifiability
No <- 100; a <- 1; b <- -1; T <- 2
Ne <- seq(1, 10, l=8)

curve(No*(1-exp(a*(b*x-T))), 0, 10)
abline(h=No*(1-exp(a*(b*0-T # intercept

C <- a*b; D <- a*T
curve(No*(1-exp(C*x-D)), add=TRUE, lty=2, col=2, lwd=2)

help(SSasympOff, help_type="html")

# model speficication
nls(Ne~No*(1-exp(a*(b*Ne-T))), start=list(a=1.2, b=0.015, T=24))
^ ^
^ ^
y y?

> y <- No*(1-exp(a*(b*Ne-T)))+rnorm(Ne,0,0.1)
> plot(y~Ne)
> nls(y~No*(1-exp(a*(b*Ne-T))), start=list(No=No, a=a, b=b, T=T))
Erro em numericDeriv(form[[3L]], names(ind), env) : 
  Obtido valor faltante ou infinito quando avaliando o modelo
> nls(y~No*(1-exp(C*Ne-D)), start=list(No=No, C=C, D=D))
Nonlinear regression model
  model:  y ~ No * (1 - exp(C * Ne - D)) 
   data:  parent.frame() 
 No   C   D 
99.9763 -0.9972  2.0197 
 residual sum-of-squares: 0.02746

Number of iterations to convergence: 2 
Achieved convergence tolerance: 6.673e-07 
> 

Bests.

-
..ooo0
...
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://r.789695.n4.nabble.com/NLS-Singular-Gradient-Error-tp2069029p2073452.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] superscript in ylab

2009-10-29 Thread Walmes Zeviani

Try

ylab=expression(Temperature*degree*C))

type demo(mathplot) at R prompt for more customizations.

Walmes Zeviani
Lavras - MG, Brasil.



Lathouri, Maria wrote:
> 
> Dear all
> 
> I am doing some plots in R.
> 
> I want to have as label in y-axis Temperature (oC). I have used
> ylab=expression(paste({Temperature} ^o*C)) but what I get is
> TemperatureoC.
> How can I have a space between Temperature and the units and also the
> units to be in brackets?
> 
> Many thanks
> Maria
> 
> 
> 
>   [[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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/superscript-in-ylab-tp26099304p26099677.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Titles on lattice colorkey

2009-11-02 Thread Walmes Zeviani



r-help.20.trevva wrote:
> 
> Dear R-ers,
> 
> I'm not sure if this is a missing feature, a support request, or stupidity
> on my part, but nevertheless, its a question. Is it possible to add titles
> to colorkey legends? As far as I can tell, there is a command to do it for
> normal "key" legends, but not for "colorkeys". 
> 
> eg it works for a normal key, created through auto.key
> 
> xyplot(decrease ~ treatment, OrchardSprays, groups = rowpos,
>type = "a",
>auto.key = list(space = "right", points = FALSE, lines =
> TRUE,title="Key title"))
> 
> but there is no comparable command for a colorkey
> 
> x <- seq(pi/4, 5 * pi, length = 100)
> y <- seq(pi/4, 5 * pi, length = 100)
> r <- as.vector(sqrt(outer(x^2, y^2, "+")))
> grid <- expand.grid(x=x, y=y)
> grid$z <- cos(r^2) * exp(-r/(pi^3))
> levelplot(z~x*y, grid, cuts = 50, scales=list(log="e"), xlab="",
>   ylab="", main="Weird Function", sub="with log scales",
>   region = TRUE,
>   colorkey = list(space="right",title="Doesn't work"))
> 
> 
> Cheers,
> 
> 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.
> 
> 

Maybe using library(grid) we can solve. Try

levelplot(z~x*y, grid, cuts = 50, scales=list(log="e"), xlab="",
  ylab="", main="Weird Function", sub="with log scales",
  region = TRUE,
  colorkey=list(space="right", title="Doesn't work"),
  legend=list(top=list(fun=grid::textGrob("Size\n(m)", x=1.09

Walmes Zeviani.
-- 
View this message in context: 
http://old.nabble.com/Titles-on-lattice-colorkey-tp22868692p26156677.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Trellis Plot

2009-11-24 Thread Walmes Zeviani

You could try adapt this code:

yy <- c(rnorm(20,2),rnorm(35,3),
rnorm(30,2),rnorm(20,3),rnorm(4,2),rnorm(10,3))
xx <- c(1:20,1:35,1:30,1:20,1:4,1:10)
gg <- rep(c('A','B','A','B','A','B'), c(20,35,30,20,4,10))
pp <- rep(c('Cond 1','Cond 2','Cond 3'), c(55, 50, 14))

#pdf("teste.pdf")
xyplot(yy ~ xx | pp, groups=gg)
 trellis.focus('strip', 1, 1, highlight=FALSE)
  ltext(0, .5, '20', col='red', pos=4)
  ltext(1, .5, '35', col='black', pos=2)
 trellis.unfocus()
 trellis.focus('strip', 2, 1, highlight=FALSE)
  ltext(0,.5,'30',col='red', pos=4)
  ltext(1,.5,'20',col='black', pos=2)
 trellis.unfocus()
 trellis.focus('strip', 1, 2, highlight=FALSE)
  ltext(0,.5,'4',col='red', pos=4)
  ltext(1,.5,'10',col='black', pos=2)
 trellis.unfocus()
#dev.off()

Walmes Zeviani, Brasil.





ychu066 wrote:
> 
> anyone know how to add text in the Trellis plot panel ?? i want to add
> things eg: dot dot dot. in the headrer of the panel.
> 
> eg:  http://old.nabble.com/file/p26486579/hist1.png hist1.png 
> 

-- 
View this message in context: 
http://old.nabble.com/Trellis-Plot-tp26486579p26494827.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Piecewise regression in lmer

2010-01-04 Thread Walmes Zeviani



AD Hayward wrote:
> 
> Dear all,
> 
> I'm attempting to use a piecewise regression to model the trajectory  
> of reproductive traits with age in a longitudinal data set using a  
> mixed model framework. The aim is to find three slopes and two points-  
> the slope from low performance in early age to a point of high  
> performance in middle age, the slope (may be 0) of the plateau from  
> the start of high performance to the end of high performance , and the  
> slope of the decline from the end of high performance to the end of  
> life.
> 
> I've found the segmented package useful, but it cannot be implemented  
> in a mixed model framework. I've also attempted piecewise regression  
> using this formula in lmer:
> 
> m<-lmer(repro ~ OTHER FIXED EFFECTS + age*(age < 2) + age*(age >= 2 &  
> age < 8) + age*(age >= 8) + (1|id) + (1|yr), data = reproduction,  
> family = binomial, link = "logit", GHQ = TRUE)
> 
> However, this gives the warning:
> 
> Warning message:
> In mer_finalize(ans) : gr cannot be computed at initial par (65)
> 
> which is not apparent if I use just two break points or I implement  
> the model in glm.
> 
> My question is essentially whether anyone can recommend a method for  
> performing piecewise regression in lmer or another mixed model  
> framework. Any advice would be greatly appreciated.
> 
> Regards,
> 
> Adam
> 
> 
> 
> Adam Hayward
> PhD Student
> Wild Evolution Group
> Institute of Evolutionary Biology
> Room 133 Ashworth Laboratories
> King's Buildings
> University of Edinburgh
> West Mains Road
> Edinburgh
> EH9 3JT
> 
> adam.hayw...@ed.ac.uk
> http://wildevolution.biology.ed.ac.uk/lkruuk/AdamHayward.html
> 
> 
> -- 
> The University of Edinburgh is a charitable body, registered in
> Scotland, with registration number SC005336.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 


Adam,

A segmented linear model, for estimation purposes, is a nonlinear model. It
requires a iteractive procedure for estimation of fixed effects. You could
use nlmer() for this.

Walmes Zeviani, Lavras - MG, Brasil.


-
..oooO
..
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/Piecewise-regression-in-lmer-tp998154p998229.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] multivariate group means

2010-01-04 Thread Walmes Zeviani



Ondřej Mikula wrote:
> 
> Hello,
> I look for a simple command computing multivariate group means and
> returning an object of class "matrix" rather than "list". Does any
> such function exist in standard packages?
> I'm beginning with R, so I'm sorry if the solution is trivial.
> Ondra Mikula
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

Ondra,

Look at kmeans() procedure (help(kmeans)). I think this is what you're
looking for.



-
..oooO
......
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/multivariate-group-means-tp990945p998292.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Adding a distance scale to a plot?

2010-01-04 Thread Walmes Zeviani



JustADude wrote:
> 
> Do you know what steps I need to take to add a scale to a plot? 
> 
> I'm pulling my example out of "An Introduction to R: Software for
> StatisticalModelling & Computing" (see the R code around Figure 76).
> 
> 
> I would like to add a scale to the image produced by the following code. 
> I would like to the scale to list the distance and units: 
> 
> data(volcano)
> x <- 10*(1:nrow(volcano))
> y <- 10*(1:ncol(volcano))
> 
> image(x, y, volcano,col = terrain.colors(100),axes = FALSE)
> contour(x, y, volcano, levels = seq(90, 200, by=5),add = TRUE, col =
> "peru")
> axis(1, at = seq(100, 800, by = 100))
> axis(2, at = seq(100, 600, by = 100))
> box()
> title(main = "Maunga Whau Volcano",font.main = 4)
> 
> 
> Thanks for any feedback and insights.
> 
> 
> 
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 


JustADude,

Take a look at help documentation of spplot() function on the sp package.
This function offer the scale mark on the plot.



-
..oooO
..
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/Adding-a-distance-scale-to-a-plot-tp998399p998444.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Barchart bar lengths not proportionate

2010-01-14 Thread Walmes Zeviani

Rex,

I think this problem can be solved using xlim()/ylim() argument. Look at the
follwing code:

require(lattice)
da <- expand.grid(A=c("a","b"), x=1:4)
da$y <- c(1,5,6,3,2,0,6,0)

barchart(y~x|A, data=da, horizontal=FALSE)

barchart(y~x|A, data=da, horizontal=FALSE,
 ylim=c(0, 1.05*max(da$y)))

At your disposal.
Walmes.



Rex C. Eastbourne wrote:
> 
> When I use barchart (with default formatting options), I get bars whose
> lengths/heights are not proportional to their value. For example:
> 
> http://drop.io/wbagm6s/asset/capture-png
> 
> Many of the values in this chart are 1; however, because the blue bars
> extend to the left of the "0" tick mark, those bars appear to represent
> higher numeric values. Is there a way to make the length of the bar
> proportional to the data value, so that people looking at my chart are not
> misled?
> 
> Thanks,
> 
> Rex
> 
>   [[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.
> 
> 


-
..oooO
......
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/Barchart-bar-lengths-not-proportionate-tp1013702p1013918.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Hierarchical Linear Model using lme4's lmer

2010-01-16 Thread Walmes Zeviani

Doug,

It appears you are mixing nlme and lme4 formulation type.
On nlme library you type

lme(y~x, random=~1|subjetc)

On lme4 library you type

lmer(y~x+(1|subject))

You mixed them.

At your disposal.
Walmes.


Doug Adams wrote:
> 
> Hi,
> 
> I was wondering:  I've got a dataset where I've got student 'project's
> nested within 'school's, and 'division' (elementary, junior, or
> senior) at the student project level.  (Division is at the student
> level and not nested within schools because some students are
> registered as juniors & others as seniors within the same school.)
> 
> So schools are random, division is fixed, and the student Score is the
> outcome variable.  This is what I've tried:
> 
> lmer(data=Age6m, Score ~ division + (1|school), random=~1 | school)
> 
> Am I on the right track?  Thanks everyone,   :)
> 
> Doug Adams
> MStat Student
> University of Utah
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 


-
..oooO
......
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/Hierarchical-Linear-Model-using-lme4-s-lmer-tp1015485p1015621.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Remove term from formula for predict.lm

2010-01-19 Thread Walmes Zeviani

Werner,

You could set 0 to that regressors you don't want to consider for
prediction.

da <- expand.grid(A=1:20, B=rnorm(20, 4, 0.2), C=10^seq(1,2,l=20))
da$y <- rnorm(da$A, 0, 0.3)

m0 <- lm(y~A+B+C, data=da)

new <- da
new$C <- 0

predict(m0)[1:5]
predict(m0, newdata=new)[1:5]

At your disposal.
Walmes.

-
..oooO
..
..()... 0ooo...      Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/Remove-term-from-formula-for-predict-lm-tp1017686p1017759.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Model

2010-01-19 Thread Walmes Zeviani

Jim,

Did you read the posting guide?
Did you do a google search, for example, with terms like "[R] generalized
linear models", "[R] count models", "[R] poisson regression"?
I think you should do.

Walmes.

-
..oooO
..
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/Re-Model-tp1017712p1017772.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] confidence intervals for mean (GLM)

2010-01-22 Thread Walmes Zeviani

Access the information in

help(predict)
help(fitted)

With your reproducible code, run the following

predict(model1)
predict(model1, type="response")
predict(model1, type="response", se.fit=TRUE)

require(lattice)

xyplot(value+fitted(model1)~treatments, data=mydata1,
   distribute.type=TRUE, type=c("p","a"))

Explore the documentation and examples showed there.

At your disposal.
Walmes.

-
..oooO
..
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/confidence-intervals-for-mean-GLM-tp1100091p1100132.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] About LU decomposition in R

2010-01-23 Thread Walmes Zeviani

Read the Matrix package documentation

require(Matrix)
help(lu, html=TRUE)

Walmes.

-
..oooO
..
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/About-LU-decomposition-in-R-tp1229137p1288264.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] problem with "nls" function

2010-01-25 Thread Walmes Zeviani

I supose you are following the SAS formulation style. R has a different
formulation style, such as:

da <- expand.grid(A=factor(1:3), x=1:10)
da$y <- as.numeric(da$A)*da$x/(1.2+da$x)+rnorm(da$x, 0, 0.1)

m0 <- nls(y~Asym[A]*x/(Time[A]+x), data=da, start=list(Asym=c(1,2,3),
Time=c(1,1,1)))
summary(m0)

At your disposal.
Walmes.

-
..oooO
..
..()... 0ooo...      Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://n4.nabble.com/problem-with-nls-function-tp1290020p1290034.html
Sent from the R help mailing list archive at Nabble.com.

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


[R] Diference in results from doBy::popMeans, multcomp::glht and contrast::contrast for a lme model

2012-11-05 Thread Walmes Zeviani
Hello R users,

I'm analyzing an experiment in a balanced incomplet block design (BIB). The
effect of blocks are assumed to be random, so I'm using nlme::lme for this.
I'm analysing another more complex experiments and I notice some diferences
from doBy::popMeans() compared multcomp::glht() and contrast::contrast().
In my example, glht() and contrast() were equal I suspect popMeans() are
doing something diferent. This a proprosital difference? Does popMeans
account the variance of random effects or something similar? The code below
is reproducible, for easy evaluation the principal results are appended.

Thank you.

# reading BIB data
da <- read.table("http://www.leg.ufpr.br/~walmes/data/pimentel_bib3.txt";,
+  header=TRUE, sep="\t",
colClasses=c("factor","factor",NA))
str(da)
## 'data.frame':30 obs. of  3 variables:
##  $ bloco: Factor w/ 10 levels "1","10","2","3",..: 1 1 1 3 3 3 4 4 4
5 ...
##  $ variedade: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 1 2 4 1 2 5 1
...
##  $ y: int  35 28 27 30 20 22 28 16 18 36 ...

# fixed effect model
m0 <- lm(y~bloco+variedade, da)
summary(m0)
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept)  35.7556 1.0004  35.742  < 2e-16 ***
## bloco10  -0.3778 1.2729  -0.297 0.770447
## ...
## variedade2   -9.2000 0.9262  -9.933 3.01e-08 ***
## variedade3   -8.0667 0.9262  -8.710 1.81e-07 ***
## ...

require(contrast)

contr <- contrast(m0, type="average", list(bloco=levels(da$bloco),
variedade="1"))
contr
##   Contrast  S.E.LowerUpper t df Pr(>|t|)
## 1 32.76667 0.6438886 31.40168 34.13165 50.89 160

require(doBy)

pop <- popMeans(m0, effect="variedade")
pop
##   beta0 Estimate Std.Error  t.value DF Pr(>|t|)LowerUpper
variedade
## 1 0 32.76667 0.6438886 50.88872 160 31.40168
34.13165 1

summary(glht(m0, linfct=matrix(contr$X, nrow=1)))
## Linear Hypotheses:
##Estimate Std. Error t value Pr(>|t|)
## 1 == 0  32.7667 0.6439   50.89   <2e-16 ***

# ok, results are the same!

# random effects model
require(nlme)

mm0 <- lme(y~variedade, random=~1|bloco, da)
summary(mm0)
## Fixed effects: y ~ variedade
##Value Std.Error DF   t-value p-value
## (Intercept) 32.67807 1.5887794 16 20.568037   0
## variedade2  -9.08144 0.9236918 16 -9.831679   0

contr <- contrast(mm0, list(variedade="1"))
contr
##  Contrast S.E.LowerUpper t df Pr(>|t|)
##  32.67807 1.588779 29.56412 35.79202 20.57 230

pop <- popMeans(mm0, effect="variedade")
pop
##   beta0 Estimate Std.Error  t.value DF Pr(>|t|)LowerUpper
variedade
## 1 0 30.5  1.918043 15.90162 250 26.54972
34.45028 1

summary(glht(mm0, linfct=matrix(contr$X, nrow=1)))
## Linear Hypotheses:
##Estimate Std. Error z value Pr(>|z|)
## 1 == 0   32.678  1.589   20.57   <2e-16 ***

# sdt error calculated by hand
sqrt(contr$X%*%vcov(mm0)%*%t(contr$X))
##  1
## 1 1.588779

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
skype: walmeszeviani
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] barplot groups of different size i.e. height is NOT a matrix

2011-05-25 Thread Walmes Zeviani
Victor,

I agree with Marc's point of view. So, if you can use another representation
of you data, like points, considering looking at
http://lmdvr.r-forge.r-project.org/figures/figures.html figures 10.20 and
10.21 for a start point.

Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] Adjusted Rate Ratios in R

2011-05-25 Thread Walmes Zeviani
Matthew,

You can change the matrix (restriction) involved. Start from

help(contr.sum)

to know how specify this.

Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] barplot groups of different size i.e. height is NOT a matrix

2011-05-25 Thread Walmes Zeviani
You can produce a graph similar to the ggplot with lattice::barchart,

require(lattice)
dataset <- data.frame(Main=c("A","A","A","B","B"),
  Detail=c("a","b","c","1","2"),
  value=runif(5, min= 0.5, max=1))
barchart(value~Detail|Main, data=dataset,
 scales=list(x=list(relation="free")))

Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==


On Wed, May 25, 2011 at 12:04 PM, Marc Schwartz wrote:

> On May 25, 2011, at 7:56 AM, Victor Gabillon wrote:
>
> > Hello,
> >
> > I want to use the function barplot do display several group of bars.
> > A standard example is given at this link
> >
> http://onertipaday.blogspot.com/2007/05/make-many-barplot-into-one-plot.html
> >
> > But in their example the 4 groups of bars are all composed of 8 bars.
> > I want to be able do display the same kind of graph but where the number
> of bars in each group are not the same. For example the first group of bars
> would have 2 bars and the second group of bars would have 10 bars.
> >
> > barplot function has a first parameter named height which is a matrix
> where each line  are the values  for the bars of one particular group.
> > One solution could be to have a height matrix with NA values but then the
> space occupied by each group is equal to the size of the largest group!! So
> you end up with gaps (empty) where there are NAs.
> >
> > Do you know how to solve this problem?
> > Do i have to consider multiple barplots in the same plot with the same
> axis? (btw, i don't know how to do that)
> >
> > In fact the bar would represent the performance of an algorithm.
> > A group of bars would be the performance of an algorithms with different
> parameters.
> > But when comparing different algorithms it is possible that we don't want
> to display the same number of parameters for each algorithm.
> >
> > Thanks for your help.
> > Victor
>
>
> barplot() is fundamentally built upon the use of rect() to construct the
> bars, so you could always create your own variant to allow for the
> flexibility that you desire.
>
> That being said, if your performance measures (the bar heights) are other
> than discrete counts or proportions, I would advise you to consider using
> other visual presentation forms, as these are really the only two types of
> data for which barplots are generally considered satisfactory. A key to
> barplots of course is that they are based at 0 for proper visual comparison.
> Thus, if you need to have the minima of the relevant axis at a value other
> than 0, this is another reason to not use them.
>
> Even then, many folks have moved away from barplots to use point or dot
> plots and similar formats, especially where you also need to include some
> type of confidence interval for each measure.
>
> HTH,
>
> Marc Schwartz
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[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] Using deriv3() in a separated nonlinear regression model

2011-05-25 Thread Walmes Zeviani
Hi all,

I'm adjusting a nonlinear regression model for data that has a categorigal
variable present. So, I can use nls() to do this considering the categorical
variable, like this

#
da <- expand.grid(tr=gl(2,1,la=c("tr")), x=1:12)
da$y <- 10*da$x/(3+da$x)+rnorm(da$x,0,0.1)
plot(y~x, da)

n0 <- nls(y~A[tr]*x/(B[tr]+x), data=da,
  start=list(A=c(10,10), B=c(3,3)), trace=TRUE)
summary(n0)
attr(n0$m$fitted(), "gradient")
attr(n0$m$fitted(), "hessian")
#

but I want specify a model using deriv3() function because I use the
gradient and hessian of the model fit

#
model <- deriv3(~A*x/(B+x), c("A","B"), function(x, A, B) NULL)
n1 <- nls(y~model(x, A, B), data=da,
  start=list(A=c(10), B=c(3)), trace=TRUE)
summary(n1)

attr(n1$m$fitted(), "gradient")
attr(n1$m$fitted(), "hessian")
#

How can I declare that categorical variable in my adjust using deriv3()
method? My try was

#
> n2 <- nls(y~model(x, A[tr], B[tr]), data=da,
+   start=list(A=c(10,10), B=c(3,3)), trace=TRUE)
0.2446993 :  10 10  3  3
0.2366604 :  9.887467 9.916289 3.00 3.00
Erro em nls(y ~ model(x, A[tr], B[tr]), data = da, start = list(A = c(10,  :

  fator de passos 0.000488281 reduzido abaixo de 'minFactor' de 0.000976562
#

That indicates the categorical variable was accounted. But values of B
didn't change during iteration. Why?

Thanks,
Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] 3 Y axis possible?

2011-05-26 Thread Walmes Zeviani
You can use mtext()

par(mar=c(5.1,4.1,4.1,5.1))
plot(x$Time, x$y1, type='l', bty = 'c', col = 'red')
par(new = TRUE)
plot(x$Time, x$y2, type = 'l', axes = FALSE, xlab = '', ylab = '', col=
'green')
axis(4, col='green')
mtext(side=4, text="label green", line=2)
par(new = TRUE)
plot(x$Time, x$y3, type = 'l', axes = FALSE, xlab = '', ylab = '', col =
'blue')
axis(4, col='blue', line = -4)
mtext(side=4, text="label green", line=-2)

Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] 3 Y axis possible?

2011-05-26 Thread Walmes Zeviani
You can pass a vector of ticks to the axis() function that expand the actual
range,

par(mar=c(5.1,4.1,4.1,5.1))
plot(x$Time, x$y1, type='l', bty = 'c', col = 'red')
par(new = TRUE)
plot(x$Time, x$y2, type = 'l', axes = FALSE, xlab = '', ylab = '', col=
'green')
##axis(4, col='green')
axis(4, at=seq(0,12,2), col='green')
mtext(side=4, text="label green", line=2)
par(new = TRUE)
plot(x$Time, x$y3, type = 'l', axes = FALSE, xlab = '', ylab = '', col =
'blue')
##axis(4, col='blue', line = -4)
axis(4, at=seq(0,0.12,0.02), col='blue', line = -4)
mtext(side=4, text="label green", line=-2)

ti <- axTicks(4)
ti.delta <- diff(ti)[1]
ti.expan <- seq(ti[1]-ti.delta, ti[length(ti)]+ti.delta, ti.delta)
ti.expan

axis(4, at=ti.expan, col='blue', line=-8)


Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==


On Thu, May 26, 2011 at 5:22 PM, Jun Shen  wrote:

> Thanks a bunch, Walmes.
>
> One more concern, the new Y axes added do not extend all the way down to
> cross with x axis. Is there anyway to make them look like the very first Y
> axis on the left?
>
> Jun
>
> On Thu, May 26, 2011 at 1:24 PM, Walmes Zeviani 
> wrote:
>
>> You can use mtext()
>>
>> par(mar=c(5.1,4.1,4.1,5.1))
>> plot(x$Time, x$y1, type='l', bty = 'c', col = 'red')
>> par(new = TRUE)
>> plot(x$Time, x$y2, type = 'l', axes = FALSE, xlab = '', ylab = '', col=
>> 'green')
>> axis(4, col='green')
>> mtext(side=4, text="label green", line=2)
>> par(new = TRUE)
>> plot(x$Time, x$y3, type = 'l', axes = FALSE, xlab = '', ylab = '', col =
>> 'blue')
>> axis(4, col='blue', line = -4)
>> mtext(side=4, text="label green", line=-2)
>>
>> Walmes.
>>
>> ==
>> Walmes Marques Zeviani
>> LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
>> Departamento de Estatística - Universidade Federal do Paraná
>> fone: (+55) 41 3361 3573
>> VoIP: (3361 3600) 1053 1173
>> e-mail: wal...@ufpr.br
>> twitter: @walmeszeviani
>> homepage: http://www.leg.ufpr.br/~walmes
>> linux user number: 531218
>> ==
>>
>>[[alternative HTML version deleted]]
>>
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
>

[[alternative HTML version deleted]]

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


Re: [R] 3 Y axis possible?

2011-05-26 Thread Walmes Zeviani
Jun,

The par() command is to give extra margin on right side to accomodate the y
axis label.
I recognize, like Rolf, that three y axis can be cumbersome, confuse. So, I
would adopt another approach using lattice::xyplot()

require(reshape)
x <- melt(x, id="Time")
str(x)

require(lattice)
xyplot(value~Time|variable, data=x, layout=c(3,1),
   type="l", scales=list(y=list(relation="free")))

Walmes.

======
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[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] xtable with conditional formatting using \textcolor

2011-06-01 Thread Walmes Zeviani
Hello list,

I'm doing a table with scores and I want include colors to represent status
of an individual. I'm using sweave <>= and xtable but I can't
get a result I want. My attemps are

#-
# code R

da <- data.frame(id=letters[1:5], score=1:5*2)

col <- function(x){
  ifelse(x>7,
 paste("\textcolor{blue}{", formatC(x, dig=2, format="f"), "}"),
 paste("\textcolor{red}{", formatC(x, dig=2, format="f"), "}"))
}

da$score.string <- col(da$score)

require(xtable)
xtable(da[,c("id","score.string")])

#-

actual result
#-
\begin{tabular}{rll}
  \hline
 & id & score.string \\
  \hline
1 & a & extcolor\{red\}\{ 2.00 \} \\
  2 & b & extcolor\{red\}\{ 4.00 \} \\
  3 & c & extcolor\{red\}\{ 6.00 \} \\
  4 & d & extcolor\{blue\}\{ 8.00 \} \\
  5 & e & extcolor\{blue\}\{ 10.00 \} \\
   \hline
\end{tabular}
#-

desired result (lines omited to save space)
#-
1 & a & \textcolor{red}{ 2.00 } \\
2 & b & \textcolor{red}{ 4.00} \\
#-----

Any contribution will be useful. Thanks.
Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] xtable with conditional formatting using \textcolor

2011-06-01 Thread Walmes Zeviani
Marc,

Thank you very much. You gave exactly what I wanted.

Bests.
Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[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] Use line break at scrip but avoid line break on graphics

2011-06-02 Thread Walmes Zeviani
Hello list,

I have plots with long strings in main=, ylab= or xlab=. So, in I my script
I use break long lines to avoid lines hiden on my monitor and in sweave
document pages. I use graphics like this

plot(1, main="
  ")

but I would like a plot result like this

plot(1, main=" ")

I remember once I saw a meta character like "\n" that avoid this breack line

plot(1, main="\(?)
  ")

Does someone know that?

Thanks.
Walmes

======
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] logistic growth model

2011-06-06 Thread Walmes Zeviani
If you use RStudio (www.rstudio.org) you can find good initial start values
by interactive plot using manipulate() function. Look the simple code below.

age <- 1:20
k <- 3; b0 <- -5; b1 <- .5
y <-  k*exp(b0+b1*age)/(1+exp(b0+b1*age))+rnorm(age,0,0.1)
plot(y~age)
start <- list()

require(manipulate)
manipulate(
   {
 plot(y~age)
 k <- kk; b0 <- b00; b1 <- b10
 curve(k*exp(b0+b1*x)/(1+exp(b0+b1*x)), add=TRUE)
 start <<- list(k=k, b0=b0, b1=b1)
   },
   kk=slider(0,5,step=0.1,initial=2.7),
   b00=slider(-5,0,step=0.1,initial=-2.5),
   b10=slider(0,1,step=0.1,initial=0.5))

start

n0 <- nls(y~k*exp(b0+b1*age)/(1+exp(b0+b1*age)), start=start)
summary(n0)

If you don't use RStudio you can follow the procedures in this blog (that
uses gWidgetsRGtk2).

http://ridiculas.wordpress.com/2011/04/09/metodo-grafico-interativo-para-valores-iniciais-em-regressao-nao-linear/

Bests.
Walmes.

======
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==


On Mon, Jun 6, 2011 at 4:20 AM, Rubén Roa  wrote:

>
> Write the growth formula in an R script.
> Define initial par values.
> Input the size and age data.
> Plot the size and age data as points.
> Plot the growth model with the initial par values as a line.
> Play with the initial par values until you see a good agreement between the
> model (the line) and the data (the points).
> Optimise.
> Re-plot.
> Plot a residual histogram.
> Plot a residual scatterplot.
> Plot a Q-Q residual plot.
>
> HTH
>
> Rubén
>
> 
>
> Dr. Rubén Roa-Ureta
> AZTI - Tecnalia / Marine Research Unit
> Txatxarramendi Ugartea z/g
> 48395 Sukarrieta (Bizkaia)
> SPAIN
>
>
>
> > -Mensaje original-
> > De: r-help-boun...@r-project.org
> > [mailto:r-help-boun...@r-project.org] En nombre de Renalda
> > Enviado el: sábado, 04 de junio de 2011 6:17
> > Para: r-help@r-project.org
> > Asunto: [R] logistic growth model
> >
> > I want to Fit a logistic growth model for y = k
> > *eb0+b1(age)/1 + eb0+b1(age), can some one help on how to get
> > the initial coefficients b0 and b1? I need to estimate in
> > order to do the regression analysis. When I run using b0=0.5
> > and b1=3.4818, I get the following error
> >
> > 397443.8 :  0.5 3.4818
> > Error in nls(Height ~ k * exp(b1 + b2 * Age)/(1 + exp(b1 + b2
> > * Age)),  :
> >singular gradient
> > "please tell me what is wrong with my initials values, and
> > how to get the initial values"
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> > http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

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


Re: [R] nls() and nls2() behavior?

2010-05-11 Thread Walmes Zeviani

I think the problem with nls() is the model specification. Look:

> y <- c(0.4334,0.3200,0.5848,0.6214,0.3890,0.5233,0.4753,
+0.2104,0.3240,0.2827,0.3847,0.5571,0.5432,0.1326,0.3481)
> x <- c(0.3521,0.4334,0.3200,0.5848,0.6214,0.3890,0.5233,
+0.1379,0.2104,0.3240,0.3404,0.3847,0.5571,0.5432,0.1326)
> dummy <- rep(c("m1","m2","m3"), c(7,3,5))
> 
> d <- data.frame(y=y, x=x, dummy=dummy)
> 
> n0 <- nls(y~b0[dummy]+b1*x, data=d,
+   start=list(b0=c(1,1,1), b1=1))
> summary(n0)

Formula: y ~ b0[dummy] + b1 * x

Parameters:
Estimate Std. Error t value Pr(>|t|)   
b01  0.486160.13993   3.474  0.00520 **
b02  0.276260.09970   2.771  0.01820 * 
b03  0.399940.12597   3.175  0.00884 **
b1  -0.017350.28352  -0.061  0.95230   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.1331 on 11 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 2.713e-07 

I don't know abou nls2().

Walmes.

-
..ooo0
...
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://r.789695.n4.nabble.com/nls-and-nls2-behavior-tp2173931p2173978.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] splitting a factor in an analysis of deviance table (negative binomial model)

2010-05-28 Thread Walmes Zeviani

Sorry by the delay. You could do:

> my.data <- expand.grid(A=factor(1:4), B=factor(1:4), rep=1:4)
> my.data$y <- rbinom(my.data$A, 10, 0.5)
> 
> model <- glm(cbind(y, 10-y)~A*B, family=binomial, data=my.data)
> anova(model, test="Chisq" )
Analysis of Deviance Table

Model: binomial, link: logit

Response: cbind(y, 10 - y)

Terms added sequentially (first to last)


 Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL63 53.528  
A 3   1.679260 51.8490.6416
B 3   2.938057 48.9110.4013
A:B   9  13.838048 35.0730.1282
> 
> X <- model.matrix(~A/B, my.data)
> 
> A <- X[,2:4]
> B.A1 <- X[,grep("A1:", colnames(X))]
> B.A2 <- X[,grep("A2:", colnames(X))]
> B.A3 <- X[,grep("A3:", colnames(X))]
> B.A4 <- X[,grep("A4:", colnames(X))]
> 
> model2 <- glm(cbind(y, 10-y)~A+B.A1+B.A2+B.A3+B.A4, family=binomial,
> data=my.data)
> anova(model2, test="Chisq")
Analysis of Deviance Table

Model: binomial, link: logit

Response: cbind(y, 10 - y)

Terms added sequentially (first to last)


 Df Deviance Resid. Df Resid. Dev P(>|Chi|)  
NULL63 53.528
A 3   1.679260 51.849   0.64156  
B.A1  3   1.292457 50.556   0.73093  
B.A2  3   4.312654 46.244   0.22962  
B.A3  3   1.886551 44.357   0.59629  
B.A4  3   9.284448 35.073   0.02574 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 

Walmes.

-
..ooo0
...
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://r.789695.n4.nabble.com/splitting-a-factor-in-an-analysis-of-deviance-table-negative-binomial-model-tp1017744p2235345.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] using the design matrix to correctly configure contrasts

2010-06-02 Thread Walmes Zeviani

Take a look at contrast package. It has functions to handle with user
contrasts. Below you have a reproducible and minimal example using
contrast's functions.

da <- expand.grid(A=factor(paste("A", 1:3, sep="")),
  B=factor(paste("B", 1:4, sep="")), rep=1:4)
eta <- model.matrix(~A*B, da)%*%matrix(c(0, 1,2, -1,0,1, 0,1,2, 0,0,0))
da$y <- rnorm(da$A, mean=eta, sd=1)

g0 <- lm(y~A*B, data=da)
anova(g0)

require(contrast)

c0 <- contrast(g0,
   list(B="B1", A="A1"),
   list(B="B2", A="A1"))
c0
c0$X

Walmes.

-
..ooo0
...
..()... 0ooo...  Walmes Zeviani
...\..(.(.)... Master in Statistics and Agricultural
Experimentation
\_). )../   walmeszevi...@hotmail.com, Lavras - MG, Brasil

(_/
-- 
View this message in context: 
http://r.789695.n4.nabble.com/using-the-design-matrix-to-correctly-configure-contrasts-tp2238719p2240372.html
Sent from the R help mailing list archive at Nabble.com.

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


[R] How set lm() to don't return NA in summary()?

2011-08-05 Thread Walmes Zeviani
Hi,

I've data from an incomplete fatorial design. One level of a factor doesn't
has the levels of the other. When I use lm(), the summary() return NA for
that non estimable parameters. Ok, I understant it. But I use
contrast::contrast(), gmodels::estimable(), multcomp::glht() and all these
fail when model has NA estimates. This is becouse vcov() and coef() has
different dimensions. Is possible set lm() to omit NA? Below same toy data
and code.

> # toy data
> adi <- expand.grid(cult=gl(1,3,la=LETTERS[1]), fert=101)
> fat <- expand.grid(cult=gl(2,3,la=LETTERS[2:3]), fert=seq(50,150,50))
> da <- rbind(adi, fat)
> da$y <- da$fert+rnorm(nrow(da),0,10)
>
> # plot
> require(lattice)
> xyplot(y~fert|cult, da)
>
> # adjust
> m0 <- lm(y~cult*fert, da)
> summary(m0)

...
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept)  7.55401   10.18956   0.7410.469
cultB   -8.04486   13.54672  -0.5940.561
cultC   -3.836446.74848  -0.5680.578
fert 0.874860.08265  10.586 1.24e-08 ***
cultB:fert   0.135090.11688   1.1560.265
cultC:fertNA NA  NA   NA
...
> require(gmodels)
> estimable(m0, cm=c(1,0,0,100,0,0))
Erro em estimable.default(m0, cm = c(1, 0, 0, 100, 0, 0)) :
 Dimension of structure(c(1, 0, 0, 100, 0, 0), .Dim = c(1L, 6L)): 1x6, not
compatible with no of parameters in m0: 5

Thanks.
Walmes.

======
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[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] multcomp::glht() doesn't work for an incomplete factorial using aov()?

2011-08-06 Thread Walmes Zeviani
Hi R users,

I sent a message yesterday about NA in model estimates (
http://r.789695.n4.nabble.com/How-set-lm-to-don-t-return-NA-in-summary-td3722587.html).
If I use aov() instead of lm() I get no NA in model estimates and I use
gmodels::estimable() without problems. Ok!
Now I'm performing a lot of contrasts and I need correcting for
multiplicity. So, I can use multcomp::glht() for this. However, glht()
return an error message that is not compatible with my expectations. Someone
know I or has a suggestion for? Below some reproducible code.

# toy data
adi <- expand.grid(cult=gl(1,3,la=LETTERS[1]), fert=101)
fat <- expand.grid(cult=gl(2,3,la=LETTERS[2:3]), fert=seq(50,150,50))
da <- rbind(adi, fat)
da$y <- da$fert+rnorm(nrow(da),0,10)

# plot
require(lattice)
xyplot(y~fert|cult, da)

# adjust complete factorial
m0 <- aov(y~cult*fert, subset(da, cult!="A"))
summary(m0)
coef(m0)

# adjust incomplete factorial
m1 <- aov(y~cult*fert, da)
summary(m1)
coef(m1)

require(multcomp)
glht(m0, linfct=matrix(c(1,1,10,0), nrow=1))   # work
glht(m1, linfct=matrix(c(1,1,0,10,0), nrow=1)) # don't work

Thank you.
Walmes.

======
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] Legend in lattice

2011-06-15 Thread Walmes Zeviani
Júlio,

Your code is not reproducible, you doesn't provide any data. So I did a
minimal code that illustrates a possible procedure is the following

n <- 30
da <- data.frame(x=runif(n), y=runif(n), z=runif(n))
da$z <- cut(da$z, seq(0,1,0.25))

require(lattice)
xyplot(y~x, da, cex=as.numeric(da$z), col=1,
   key=list(points=list(cex=1:4, pch=1), text=list(levels(da$z)),
columns=4))

Bests,
Walmes.

======
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] Error "singular gradient matrix at initial parameter estimates" in nls

2011-06-30 Thread Walmes Zeviani
I think that these aren't good initial values. If you do a plot of data and
add a curve, the curve don't approximate the data. Frequently I use
interactive procedures to get good initial values. Using playwith() you can
handle sliders to adjust values and use in nls(), look the following

r <-
c(1.16,1.143,1.109,1.093,1.079,1.066,1.053,1.040,1.027,1.015,1.004,0.994,0.985,0.977)
D <-
c(0.1806551,0.2703113,0.3757225,0.5271811,0.8665835,1.0812568,1.0612762,1.0726612,
   1.167927,1.191092,1.1336938,1.1215107,0.9619603,0.8315467)
names <- c(La,Ce,Pr,Nd,Sm,Eu,Gd,Tb,Dy,Ho,Er,Tm,Yb,Lu)

plot(log(D)~r)
Do <- 0.8; En <- 390; ro <- 1.03
curve(log10(Do)*exp(((-4*pi*En*Na)*((ro/2)*(x-ro)^2+(1/3)*(x-ro)^3))/(R*T)),
add=TRUE)

Na <- 6.0221415*10^23 # Avrogadro's number
T <- 1010 # Temp in K
R <- 8.3144   # Gas constant [J mol^-1 K^-1]

## Brice-Model
bricemod <-
nls(log10(D)~log10(Do)*exp(((-4*pi*En*Na)*((ro/2)*(r-ro)^2+(1/3)*(r-ro)^3))/(R*T)),
start=list(Do=0.8, En=390, ro=1.03), trace=TRUE)

require(playwith)

start <- list()
playwith(
 {
   plot(log(D)~r)

curve(log10(Do)*exp(((-4*pi*En*Na)*((ro/2)*(x-ro)^2+(1/3)*(x-ro)^3))/(R*T)),
add=TRUE)
   start <<- list(Do=Do, En=En, ro=ro)
 },
 parameters=list(
   Do=seq(0.7,0.8,by=0.001),
   En=seq(350,450,by=1),
   ro=seq(0.5,1.5,by=0.1)))
start

bricemod <-
nls(log10(D)~log10(Do)*exp(((-4*pi*En*Na)*((ro/2)*(r-ro)^2+(1/3)*(r-ro)^3))/(R*T)),
start=start, trace=TRUE)

Bests.
Walmes.

======
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[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] Tables and merge

2011-07-05 Thread Walmes Zeviani
Silvano,

I have some examples using merge() from my class notes in
http://www.leg.ufpr.br/doku.php/disciplinas:ce223-2011-01. See "aula11.R".
For the moment, this minimal reproducible code can be useful

id <- 1:30
n <- 20

a1 <- data.frame(id=sample(id, n), v1=rnorm(n))
a2 <- data.frame(id=sample(id, n), v2=rpois(n,10))
a3 <- data.frame(id=sample(id, n), v3=runif(n))

merge(a1, a2, by="id") # just two data.frame at once

a0 <- list(a1, a2, a3)

Reduce(function(x, y) merge(x, y, by="id"), a0, accumulate=FALSE)
# font: http://rwiki.sciviews.org/doku.php?id=tips:data-frames:merge

You can also join to the R-br mailing list (brazilian R-help list).
Instructions in http://www.leg.ufpr.br/doku.php/software:rbr

Bests.
Walmes.

======
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] Odp: Multiple plots in one subplot

2011-12-16 Thread Walmes Zeviani
You can use a different way of split the plotting area that is by means of
layout() function.

x <- rnorm(100)
M <- matrix(c(rep(1:5, e=2), 6, 7), byrow=TRUE, nrow=2)
layout(M)
plot(x)
hist(x)
qqnorm(x)
boxplot(x)
plot(density(x))
plot(abs(x))
hist(abs(x))

Bests.
Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[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] Lattice densityplot with semitransparent filled regions

2012-04-11 Thread Walmes Zeviani
Hello,

I'm doing some graphics for a paper and a need customize such with filled
region above the density curve. My attempts I get something very near what
I need, but I don't solve the problem of use semitransparent filled. Below
a minimal reproducible code. Someone has any idea?

require(lattice)

# toy data...
dt <- expand.grid(A=1:2, B=1:3, y=1:50)
dt$y <- rnorm(nrow(dt), dt$B, dt$A)

# regular plot...
densityplot(~y|B, groups=A, data=dt, plot.points="rug")

# the actual panel...
panel.densityplot

# so, I edit this...
my.panel.densityplot <-
function (x, darg = list(n = 30), plot.points = "jitter", ref = FALSE,
groups = NULL, weights = NULL, jitter.amount = 0.01 *
diff(current.panel.limits()$ylim),
type = "p", ..., identifier = "density")
{
if (ref) {
reference.line <- trellis.par.get("reference.line")
panel.abline(h = 0, col = reference.line$col, lty =
reference.line$lty,
lwd = reference.line$lwd, identifier = paste(identifier,
"abline"))
}
if (!is.null(groups)) {
panel.superpose(x, darg = darg, plot.points = plot.points,
ref = FALSE, groups = groups, weights = weights,
panel.groups = panel.densityplot, jitter.amount =
jitter.amount, # alterei para my.panel
type = type, ...)
}
else {
switch(as.character(plot.points), `TRUE` = panel.xyplot(x = x,
y = rep(0, length(x)), type = type, ..., identifier =
identifier),
rug = panel.rug(x = x, start = 0, end = 0, x.units = c("npc",
"native"), type = type, ..., identifier = paste(identifier,
"rug")), jitter = panel.xyplot(x = x, y = jitter(rep(0,
length(x)), amount = jitter.amount), type = type,
..., identifier = identifier))
density.fun <- function(x, weights, subscripts = TRUE,
darg, ...) {
do.call("density", c(list(x = x, weights =
weights[subscripts]),
darg))
}
if (sum(!is.na(x)) > 1) {
h <- density.fun(x = x, weights = weights, ..., darg = darg)
lim <- current.panel.limits()$xlim
id <- h$x > min(lim) & h$x < max(lim)
panel.lines(x = h$x[id], y = h$y[id], ..., identifier =
identifier)
## line above was added
panel.polygon(x=h$x[id], y = h$y[id], ..., identifier =
identifier, alpha=0.2)
}
}
}

# my customized plot, I want semitransparent colors
# and use the colors of trellis.par.set("superpose.polygon") to fill
densityplot(~y|B, groups=A, data=dt,
plot.points="rug", col=2:3,
panel=panel.superpose,
panel.groups=my.panel.densityplot)

Thanks!
Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] Lattice densityplot with semitransparent filled regions

2012-04-11 Thread Walmes Zeviani
Thank you Ilai. Problem solved. There is a small detail, alpha affects the
rug and the curve line opacity too. Is possible to specify it just to
polygon?

Bests.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] Encoding of Sweave file error message

2012-04-11 Thread Walmes Zeviani
I had the same problem! So, as I'm a linux user, I prefer use linux
terminal. On terminal I type this to compile

R CMD Sweave --encoding=utf-8 myfile.Rnw

and the compilation is successful. Try to set the encoding option in
Sweave().

Bests.
Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] How to test if a slope is different than 1?

2012-04-25 Thread Walmes Zeviani
You can also run two nls() models, one under h0 restriction, other under no
restriction or h1, and compare them (if they are nested) by likelihood
ratio test using anova() method, look

x1 <- seq(0,10,l=15)
x2 <- runif(x1)
set.seed(1)
y <- x1+0.5*x2+rnorm(x1,0,0.01)

nls.h0 <- nls(y~b0+x1+b2*x2, start=c(b0=1, b2=1))
nls.h1 <- nls(y~b0+b1*x1+b2*x2, start=c(b0=1, b1=1, b2=1))
summary(nls.h1)
anova(nls.h0, nls.h1)

another option is to adjust a model parametrized according to test h0, like

nls.h2 <- nls(y~b0+(b1+1)*x1+b2*x2, start=c(b0=1, b1=-1, b2=1))
summary(nls.h2)

Bests.
Walmes

======
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] VarCorr procedure from lme4

2012-05-01 Thread Walmes Zeviani
It could be a bad coexistence between packages in the same R session. Are
you using nlme and/or doBy packages too?

Bests.
Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] error bars for a barchart

2012-05-01 Thread Walmes Zeviani
I have a repoducibe example here

http://ridiculas.wordpress.com/2011/11/23/media-e-desvio-padrao-de-muitas-variaveis-separado-por-grupos/

Sorry for it be in Portuguese.

Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] confidence intervals for nls or nls2 model

2012-05-16 Thread Walmes Zeviani
If you want a confidence based in new x values you can do. I have this post
with steps to do this. It's written in Portuguese but the R code is useful.

http://ridiculas.wordpress.com/2011/05/19/bandas-de-confianca-para-modelo-de-regressao-nao-linear/

Bests.
Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


Re: [R] question how to add Standard Deviation as "Whiskers" in a simple plot

2012-05-28 Thread Walmes Zeviani
This post are useful.

http://myowelt.blogspot.com.br/2008/03/beautiful-error-bars-in-r.html
http://mapas.mma.gov.br/i3geo/pacotes/rlib/win/gplots/html/plotCI.html

Walmes.

==
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==

[[alternative HTML version deleted]]

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


[R] How to provide a result from D(f(x), "x") to a curve(f'(x)) ???

2010-02-17 Thread Walmes Marques Zeviani

Hi all,

I want to provide the result from D() to curve(), because I want to plot the 
k-th derivative of some functions. Actually, I copy from console the result 
given by D() and paste inside curve(). With a lot of functions and high degree 
differentiation this process is tedious. Can I provide directly??

# what I actually have done (very simple function)
D(expression(x^3), "x") # copy this result
curve(3 * x^2) # paste inside

# my failed attempts
curve(as.expression(D(expression(x^3), "x")))
curve(as.character(as.expression(D(expression(x^3), "x"
curve(noquote(as.character(as.expression(D(expression(x^3), "x")

Thanks in advance.
Walmes Zeviani, Brasil.   
_
Quer deixar seus vídeos mais divertidos? Com o Movie Maker isso fica fácil.

ndows Live:Dicas - Movie Maker:Hotmail:Tagline:1x1:Titulo Legendas Creditos
[[alternative HTML version deleted]]

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


[R] How get the single bar x coordinate in barchart when groups is used?

2010-02-20 Thread Walmes Marques Zeviani

Hi all,

I am using barchart() to plot values above the bars. When using groups argument 
we get bars grouped arround a given x level. By placing values above this bars 
we need to know the respective x coordinates. How can I get it?

require(lattice)
da <- expand.grid(x=1:5, z=1:3, w=1:2)
da$y <- rpois(da$x, lambda=23)

barchart(y~x|w, groups=z, data=da, horizontal=FALSE,
 panel=function(x, y, subscripts, groups, ...){
   panel.barchart(x, y, subscripts=subscripts, groups=groups, ...)
   d <- 0.22 # <-- how obtain "d" or coordinates?
   panel.text(x+c(-d,0,d), y, label=y, pos=3)
 })

Thanks in advance.
Walmes Zeviani, Lavras - MG, Brasil.
  
_

os.

dium=Tagline&utm_campaign=InfuseSocial
[[alternative HTML version deleted]]

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


[R] R algorithm for maximum curvatures measures of nonlinear models

2010-03-05 Thread Walmes Marques Zeviani

Hi all,

I'm looking for the algorithm to calculate maximum intrinsic and parameter 
effects curvature of Bates & Watts (1980). I have Bates & Watts (1980) original 
article, Bates et al (1983) article, Seber & Wild (1989) and Ratkowsky (1983) 
books. Those sources show steps and algorithm to get this measures but I don't 
know translate C code to R language and I've no success until now.

I know and I use rms.curv() in MASS package but I would like the maximum 
curvatures measures.

Does someone have this implemented in R? Or know some paper illustrating this 
calculation with real data (but not so trivial) that can able me to create my 
functions?

Thanks a lot.
Walmes Zeviani.

Bates; Hamilton; Watts. Calculation onf intrinsic and parameter-effects 
curvatures for nonlinear models. Communications in statistics - simulations and 
computation, 469-477, 1983.
Bates; Watts. Relative curvature measures of nonlinearity. Journal Royal 
Statistical Society - B. 40, 1-25, 1980.
Seber; Wild. Nonlinear regression. Wiley, 1989.
Ratkowsky. Nonlinear regression modeling. Marcel Dekker, 1983.
  
_



[[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] Is possible a mini-plot into a big plot with Lattice?

2010-05-26 Thread Walmes Marques Zeviani

Hello,

I want to do with Lattice functions (qqmath, histogram) a figure like this 
below.
n <- 1000
x <- rnorm(n)
qqnorm(x); qqline(x)
op <- par(fig=c(.02,.5,.5,.98), new=TRUE)
hist(x, xlab="", ylab="", main="", axes=FALSE)
box()
par(op)

Is possible?

Thanks.
Walmes Zeviani.

  
_
CANSADO DE ENTRAR EM TODAS AS SUAS DIFERENTES CONTAS DE EMAIL? JUNTE TODAS 
AGORA.

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


Re: [R] Is possible a mini-plot into a big plot with Lattice?

2010-05-26 Thread Walmes Marques Zeviani

Mr. Murrell,

It is exactly what I want. Thanks very much.

Walmes.

> Date: Thu, 27 May 2010 10:27:02 +1200
> From: p.murr...@auckland.ac.nz
> To: walmeszevi...@hotmail.com
> CC: r-help@r-project.org
> Subject: Re: [R] Is possible a mini-plot into a big plot with Lattice?
> 
> Hi
> 
> On 27/05/2010 8:09 a.m., Walmes Marques Zeviani wrote:
> >
> > Hello,
> >
> > I want to do with Lattice functions (qqmath, histogram) a figure like this 
> > below.
> > n<- 1000
> > x<- rnorm(n)
> > qqnorm(x); qqline(x)
> > op<- par(fig=c(.02,.5,.5,.98), new=TRUE)
> > hist(x, xlab="", ylab="", main="", axes=FALSE)
> > box()
> > par(op)
> >
> > Is possible?
> 
> Something like ... ?
> 
> library(lattice)
> library(grid)
> qqmath(x,
> panel=function(...) {
> panel.qqmath(...)
> panel.abline(a=0, b=1)
> pushViewport(viewport(.02, .5, .48, .48,
>   just=c("left", "bottom")))
> print(histogram(x, xlab=NULL, ylab=NULL,
> scales=list(draw=FALSE)),
>   newpage=FALSE)
> popViewport()
> })
> 
> Paul
> 
> > Thanks.
> > Walmes Zeviani.
> >
> > 
> > _
> > CANSADO DE ENTRAR EM TODAS AS SUAS DIFERENTES CONTAS DE EMAIL? JUNTE TODAS 
> > AGORA.
> >
> > [[elided Hotmail 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.
> 
> -- 
> Dr Paul Murrell
> Department of Statistics
> The University of Auckland
> Private Bag 92019
> Auckland
> New Zealand
> 64 9 3737599 x85392
> p...@stat.auckland.ac.nz
> http://www.stat.auckland.ac.nz/~paul/
  
_
DEIXE SUAS CONVERSAS MAIS DIVERTIDAS. TRANSFORME AQUI SUAS FOTOS EM EMOTICONS, 
É GRÁTIS.

[[alternative HTML version deleted]]

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


Re: [R] ternary contour plot

2011-02-15 Thread Walmes Marques Zeviani

Colin,

If your propose is to create a ternary plot with points and vectors, I think 
easier do this with graphics based plots instead of trellis based plots. 
Although, with little work you can do with trellis too. I gave you a 
reproducible code to put an arrow in a ternary plot. I use the function 
locator() to extract coordinates. The code is the following

#--
# package and related documentaion

require(plotrix)
help.search("triax")
help(triax.plot, html=TRUE)

#--
# toy data

data(soils)
str(soils)

#--
# creating a ternary plot

triax.plot(soils[1:10,], main="DEFAULT")

#--
# extracting coodinates by mouse click on the plot (cartesian coordinates, x  
and y axis)

id <- locator(n=2) # click on the plot to extract 2 coodinates

#--
# draw an arrow with the coordinates extracted

arrows(id$x[1], id$y[1], id$x[2], id$y[2])

#------

At your disposal.
Walmes.

====


Walmes Marques Zeviani


LEG (Laboratório de Estatística e Geoinformação)


Departamento de Estatística - Universidade Federal do Paraná


fone: (+55) 41 3361 3573


VoIP: (3361 3600) 1053 1173


e-mail: wal...@ufpr.br / @walmeszeviani


homepage: http://www.leg.ufpr.br/~walmes







> Date: Mon, 14 Feb 2011 20:05:54 +
> From: colin.bl...@bristol.ac.uk
> To: walmeszevi...@hotmail.com
> Subject: ternary contour plot
> 
> Dear Walmes,
> 
> Firstly, thank you for distributing your ternary contour plot code.
> 
> could I ask you a question about it, hopefully you will be able to 
> answer it quite simply and so save me some time trying to work out what 
> it is that the function is actually doing.
> 
> Essentially I wish to overlay points and vectors onto the ternary 
> contour plot. Now can I do this by:
> 
> a): extracting the adjusted x y coordinates and using it with triax.fill
> 
> or
> 
> b) can I plot my vectors (arrows) and points onto the grid that is 
> created by levelplot.ternary
> 
> if either or both are true do you know how I would go about a) 
> extracting or b) inputing the relevent data?
> 
> I would be grateful for your advice,
> 
> best regards,
> 
> colin
> 
> -- 
> **
> **
> 
> Dr Colin Bleay
> Station d'Ecologie Experimentale du CNRS,
> 09200 Moulis,
> France.
> 
> Tel: +33 5 61 04 03 61
> Fax: +33 5 61 96 08 51
> email: colin.bl...@ecoex-moulis.cnrs.fr
> Webpage: http://www.ecoex-moulis.cnrs.fr/Staffpages/ColinBleay.htm
> 
  
[[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.