[R] Apply rmarkdown::render() outside the RStudio don't find pandoc
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
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
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
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
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
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
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...
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
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
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"
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?
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()?
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
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
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
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
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)
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
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
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
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
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
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
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
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?
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
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
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
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
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)
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
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
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
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
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
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
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
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?
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?
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?
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
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
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
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
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?
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)
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
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()?
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()?
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
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
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
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
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
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
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
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?
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
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
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
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
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)) ???
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?
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
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?
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?
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
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.