[R] twoby2 (Odds Ratio) for variables with 3 or more levels
Dear all, I am using Epi package to calculate Odds ratio in my bivariate analysis. How can I make *twoby2 *in variables that have 3 or more levels. For example: I have 4 level var (Age) m=matrix(c(290, 100,232, 201, 136, 99, 182, 240), nrow=4, ncol=2) twoby2(m) R gives me only Comparing : Row 1 vs. Row 2 While I would like to have reference value in Row 1, and compare Row 2, Row 3 and Row 4 with it. Thanks for your help! [[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] woby2 (Odds Ratio) for variables with 3 or more levels
Dear all, I am using Epi package to calculate Odds ratio in my bivariate analysis. How can I make *twoby2 *in variables that have 3 or more levels. For example: I have 4 level var (Age) m=matrix(c(290, 100,232, 201, 136, 99, 182, 240), nrow=4, ncol=2) library (Epi) twoby2(m) R gives me only Comparing : Row 1 vs. Row 2 While I would like to have reference value in Row 1, and compare Row 2, Row 3 and Row 4 with it. Thanks for your help! [[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] network plot problem
Is it possible to write a numbers of link next to lines of the network? For exmple, I have 3 trials that have studied two drugs in network. So, now it seems all links are of same weight. 2012/7/27 Vlatka Matkovic Puljic v.matkovic.pul...@gmail.com Thank you! It works fine now! :) 2012/7/27 Rui Barradas ruipbarra...@sapo.pt Hello, As for the directed = FALSE, it works with me: g - graph(given, directed = FALSE) As for the label, it's easy. Create an attribute 'name' and then set the label to that attribute. For instance, using the alphabet's letters, first uppercase. V(g)$name - c(LETTERS, letters)[V(g)] # or any other names V(g)$label - V(g)$name Keep 'name' and 'label' separate, like this you may change the label if needed but the name will still be there. Then the plot functions recognise the label attribute. Hope this helps, Rui Barradas Em 27-07-2012 11:37, Vlatka Matkovic Puljic escreveu: Thank you Rui. With matrix works better. I got plot I have expected to have. I want to be undirected. But directed = FALSE or as.undirected(graph) are not working? Where I'd gone wrong? Shouldn't it be possible to add names instead of numbers of nods with vertex.label= ? 2012/7/26 Rui Barradas ruipbarra...@sapo.pt Hello, I don't see the problem. given - scan(text= 15 2 10 4 10 4 10 4 13 4 13 4 15 4 18 4 11 5 2 6 7 6 7 6 7 6 12 6 15 6 15 6 19 6 22 6 24 6 6 7 5 12 5 12 7 12 11 12 13 12 13 12 13 12 13 12 16 12 17 12 23 12 23 12 23 12 23 12 6 13 12 13 6 14 6 15 9 15 12 15 13 15 17 16 16 17 1 18 12 18 23 18 2 19 6 19 24 19 21 22 3 25 5 26 6 27 7 27 15 29 20 30 25 31 28 31 8 32 6 33 14 33 22 34 ) mat - matrix(given, ncol=2, byrow=TRUE) g - graph(given) # Or, not run (note the transpose) #g - graph(t(mat)) # The edges are exactly what is given E(g) V(g)$label - V(g) g$layout - layout.fruchterman.reingold plot(g, edge.arrow.size=0.5, edge.loop.angle=1, edge.curved=FALSE) So my guess is you've messed up the graph creation. Also, don't post datasets like that, use dput(), lik this: dput(given) c(15, 2, 10, 4, 10, 4, 10, 4, 13, 4, 13, 4, 15, 4, 18, 4, 11, 5, 2, 6, 7, 6, 7, 6, 7, 6, 12, 6, 15, 6, 15, 6, 19, 6, 22, 6, 24, 6, 6, 7, 5, 12, 5, 12, 7, 12, 11, 12, 13, 12, 13, 12, 13, 12, 13, 12, 16, 12, 17, 12, 23, 12, 23, 12, 23, 12, 23, 12, 6, 13, 12, 13, 6, 14, 6, 15, 9, 15, 12, 15, 13, 15, 17, 16, 16, 17, 1, 18, 12, 18, 23, 18, 2, 19, 6, 19, 24, 19, 21, 22, 3, 25, 5, 26, 6, 27, 7, 27, 15, 29, 20, 30, 25, 31, 28, 31, 8, 32, 6, 33, 14, 33, 22, 34) Now all anyone has to do is to copy that output and paste it into an R session. (Try it with the matrix to see the result) Hope this helps, Rui Barradas Em 25-07-2012 16:27, Vlatka Matkovic Puljic escreveu: **Hi, I wanted to create a network of drugs that are being studied together. So, I created a file as a graph. But plot of my network is not corect! It looks like to me that something else is lying behind this network links created. Could someone help me with this? Thanx! a-read.graph(file=file.choose(), format=edgelist) plot.igraph(a) File (number present ID of drug): 15 2 10 4 10 4 10 4 13 4 13 4 15 4 18 4 11 5 2 6 7 6 7 6 7 6 12 6 15 6 15 6 19 6 22 6 24 6 6 7 5 12 5 12 7 12 11 12 13 12 13 12 13 12 13 12 16 12 17 12 23 12 23 12 23 12 23 12 6 13 12 13 6 14 6 15 9 15 12 15 13 15 17 16 16 17 1 18 12 18 23 18 2 19 6 19 24 19 21 22 3 25 5 26 6 27 7 27 15 29 20 30 25 31 28 31 8 32 6 33 14 33 22 34 -- ** Vlatka Matkovic Puljic gsm: +32.474.894953 -- ** Vlatka Matkovic Puljic gsm: +32.474.894953 [[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] network plot problem
Thank you Rui. With matrix works better. I got plot I have expected to have. I want to be undirected. But directed = FALSE or as.undirected(graph) are not working? Where I'd gone wrong? Shouldn't it be possible to add names instead of numbers of nods with vertex.label= ? 2012/7/26 Rui Barradas ruipbarra...@sapo.pt Hello, I don't see the problem. given - scan(text= 15 2 10 4 10 4 10 4 13 4 13 4 15 4 18 4 11 5 2 6 7 6 7 6 7 6 12 6 15 6 15 6 19 6 22 6 24 6 6 7 5 12 5 12 7 12 11 12 13 12 13 12 13 12 13 12 16 12 17 12 23 12 23 12 23 12 23 12 6 13 12 13 6 14 6 15 9 15 12 15 13 15 17 16 16 17 1 18 12 18 23 18 2 19 6 19 24 19 21 22 3 25 5 26 6 27 7 27 15 29 20 30 25 31 28 31 8 32 6 33 14 33 22 34 ) mat - matrix(given, ncol=2, byrow=TRUE) g - graph(given) # Or, not run (note the transpose) #g - graph(t(mat)) # The edges are exactly what is given E(g) V(g)$label - V(g) g$layout - layout.fruchterman.reingold plot(g, edge.arrow.size=0.5, edge.loop.angle=1, edge.curved=FALSE) So my guess is you've messed up the graph creation. Also, don't post datasets like that, use dput(), lik this: dput(given) c(15, 2, 10, 4, 10, 4, 10, 4, 13, 4, 13, 4, 15, 4, 18, 4, 11, 5, 2, 6, 7, 6, 7, 6, 7, 6, 12, 6, 15, 6, 15, 6, 19, 6, 22, 6, 24, 6, 6, 7, 5, 12, 5, 12, 7, 12, 11, 12, 13, 12, 13, 12, 13, 12, 13, 12, 16, 12, 17, 12, 23, 12, 23, 12, 23, 12, 23, 12, 6, 13, 12, 13, 6, 14, 6, 15, 9, 15, 12, 15, 13, 15, 17, 16, 16, 17, 1, 18, 12, 18, 23, 18, 2, 19, 6, 19, 24, 19, 21, 22, 3, 25, 5, 26, 6, 27, 7, 27, 15, 29, 20, 30, 25, 31, 28, 31, 8, 32, 6, 33, 14, 33, 22, 34) Now all anyone has to do is to copy that output and paste it into an R session. (Try it with the matrix to see the result) Hope this helps, Rui Barradas Em 25-07-2012 16:27, Vlatka Matkovic Puljic escreveu: **Hi, I wanted to create a network of drugs that are being studied together. So, I created a file as a graph. But plot of my network is not corect! It looks like to me that something else is lying behind this network links created. Could someone help me with this? Thanx! a-read.graph(file=file.**choose(), format=edgelist) plot.igraph(a) File (number present ID of drug): 15 2 10 4 10 4 10 4 13 4 13 4 15 4 18 4 11 5 2 6 7 6 7 6 7 6 12 6 15 6 15 6 19 6 22 6 24 6 6 7 5 12 5 12 7 12 11 12 13 12 13 12 13 12 13 12 16 12 17 12 23 12 23 12 23 12 23 12 6 13 12 13 6 14 6 15 9 15 12 15 13 15 17 16 16 17 1 18 12 18 23 18 2 19 6 19 24 19 21 22 3 25 5 26 6 27 7 27 15 29 20 30 25 31 28 31 8 32 6 33 14 33 22 34 -- ** Vlatka Matkovic Puljic gsm: +32.474.894953 [[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] network plot problem
Thank you! It works fine now! :) 2012/7/27 Rui Barradas ruipbarra...@sapo.pt Hello, As for the directed = FALSE, it works with me: g - graph(given, directed = FALSE) As for the label, it's easy. Create an attribute 'name' and then set the label to that attribute. For instance, using the alphabet's letters, first uppercase. V(g)$name - c(LETTERS, letters)[V(g)] # or any other names V(g)$label - V(g)$name Keep 'name' and 'label' separate, like this you may change the label if needed but the name will still be there. Then the plot functions recognise the label attribute. Hope this helps, Rui Barradas Em 27-07-2012 11:37, Vlatka Matkovic Puljic escreveu: Thank you Rui. With matrix works better. I got plot I have expected to have. I want to be undirected. But directed = FALSE or as.undirected(graph) are not working? Where I'd gone wrong? Shouldn't it be possible to add names instead of numbers of nods with vertex.label= ? 2012/7/26 Rui Barradas ruipbarra...@sapo.pt Hello, I don't see the problem. given - scan(text= 15 2 10 4 10 4 10 4 13 4 13 4 15 4 18 4 11 5 2 6 7 6 7 6 7 6 12 6 15 6 15 6 19 6 22 6 24 6 6 7 5 12 5 12 7 12 11 12 13 12 13 12 13 12 13 12 16 12 17 12 23 12 23 12 23 12 23 12 6 13 12 13 6 14 6 15 9 15 12 15 13 15 17 16 16 17 1 18 12 18 23 18 2 19 6 19 24 19 21 22 3 25 5 26 6 27 7 27 15 29 20 30 25 31 28 31 8 32 6 33 14 33 22 34 ) mat - matrix(given, ncol=2, byrow=TRUE) g - graph(given) # Or, not run (note the transpose) #g - graph(t(mat)) # The edges are exactly what is given E(g) V(g)$label - V(g) g$layout - layout.fruchterman.reingold plot(g, edge.arrow.size=0.5, edge.loop.angle=1, edge.curved=FALSE) So my guess is you've messed up the graph creation. Also, don't post datasets like that, use dput(), lik this: dput(given) c(15, 2, 10, 4, 10, 4, 10, 4, 13, 4, 13, 4, 15, 4, 18, 4, 11, 5, 2, 6, 7, 6, 7, 6, 7, 6, 12, 6, 15, 6, 15, 6, 19, 6, 22, 6, 24, 6, 6, 7, 5, 12, 5, 12, 7, 12, 11, 12, 13, 12, 13, 12, 13, 12, 13, 12, 16, 12, 17, 12, 23, 12, 23, 12, 23, 12, 23, 12, 6, 13, 12, 13, 6, 14, 6, 15, 9, 15, 12, 15, 13, 15, 17, 16, 16, 17, 1, 18, 12, 18, 23, 18, 2, 19, 6, 19, 24, 19, 21, 22, 3, 25, 5, 26, 6, 27, 7, 27, 15, 29, 20, 30, 25, 31, 28, 31, 8, 32, 6, 33, 14, 33, 22, 34) Now all anyone has to do is to copy that output and paste it into an R session. (Try it with the matrix to see the result) Hope this helps, Rui Barradas Em 25-07-2012 16:27, Vlatka Matkovic Puljic escreveu: **Hi, I wanted to create a network of drugs that are being studied together. So, I created a file as a graph. But plot of my network is not corect! It looks like to me that something else is lying behind this network links created. Could someone help me with this? Thanx! a-read.graph(file=file.choose(), format=edgelist) plot.igraph(a) File (number present ID of drug): 15 2 10 4 10 4 10 4 13 4 13 4 15 4 18 4 11 5 2 6 7 6 7 6 7 6 12 6 15 6 15 6 19 6 22 6 24 6 6 7 5 12 5 12 7 12 11 12 13 12 13 12 13 12 13 12 16 12 17 12 23 12 23 12 23 12 23 12 6 13 12 13 6 14 6 15 9 15 12 15 13 15 17 16 16 17 1 18 12 18 23 18 2 19 6 19 24 19 21 22 3 25 5 26 6 27 7 27 15 29 20 30 25 31 28 31 8 32 6 33 14 33 22 34 -- ** Vlatka Matkovic Puljic gsm: +32.474.894953 [[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] network plot problems
Hi, I wanted to create a network of drugs that are being studied together. So, I created a file as a graph. But plot of my network is not corect! It looks like to me that something else is lying behind this network links created. Could someone help me with this? Thanx! a-read.graph(file=file.choose(), format=edgelist) plot.igraph(a) File (number present ID of drug): 15 2 10 4 10 4 10 4 13 4 13 4 15 4 18 4 11 5 2 6 7 6 7 6 7 6 12 6 15 6 15 6 19 6 22 6 24 6 6 7 5 12 5 12 7 12 11 12 13 12 13 12 13 12 13 12 16 12 17 12 23 12 23 12 23 12 23 12 6 13 12 13 6 14 6 15 9 15 12 15 13 15 17 16 16 17 1 18 12 18 23 18 2 19 6 19 24 19 21 22 3 25 5 26 6 27 7 27 15 29 20 30 25 31 28 31 8 32 6 33 14 33 22 34 -- ** Vlatka Matkovic Puljic +32/ 485/ 453340 [[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] network plot problem
**Hi, I wanted to create a network of drugs that are being studied together. So, I created a file as a graph. But plot of my network is not corect! It looks like to me that something else is lying behind this network links created. Could someone help me with this? Thanx! a-read.graph(file=file.choose(), format=edgelist) plot.igraph(a) File (number present ID of drug): 15 2 10 4 10 4 10 4 13 4 13 4 15 4 18 4 11 5 2 6 7 6 7 6 7 6 12 6 15 6 15 6 19 6 22 6 24 6 6 7 5 12 5 12 7 12 11 12 13 12 13 12 13 12 13 12 16 12 17 12 23 12 23 12 23 12 23 12 6 13 12 13 6 14 6 15 9 15 12 15 13 15 17 16 16 17 1 18 12 18 23 18 2 19 6 19 24 19 21 22 3 25 5 26 6 27 7 27 15 29 20 30 25 31 28 31 8 32 6 33 14 33 22 34 -- ** Vlatka Matkovic Puljic +32/ 485/ 453340 [[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] betareg help
Yes, data is extremely skewed (as all AIDS behavioral data). Thank you very much for time and solutions! 2011/3/13 Achim Zeileis achim.zeil...@uibk.ac.at Ben, thanks for your analysis...I also just sent a message with some similar (and some different) ideas. On Sun, 13 Mar 2011, Ben Bolker wrote: The problem seems to be that the algorithm for coming up with a starting guess for the phi (dispersion) parameter is getting a negative number. Yes, as I had written earlier today for the regression on an intercept only. It's not all that easy to figure this out ... Indeed. I've already added a better warning and an ad hoc workaround to the devel version of betareg. thx, Z The data set is a little bit nasty (lots of points stacked on the equivalent of (0,0)), but not pathological. I'm cc'ing the maintainer of the package -- adding the lines if (sum(phi_y)0) { stop(bad estimated start value for phi: consider setting start values manually\n, (see ?betareg.control)\n, Estimated starting values for mean:\n,paste(beta,collapse=,)) } at line 162 of betareg.R in the current CRAN version provides a more informative error message in this case. See below for solutions. = results - read.csv2(betareg_tmp.csv) results$drugcat - cut(results$drug,c(0,0.005,0.06,0.17)) table(results$drug) table(results$alcoh) table(results$cond) ## shows that fairly large fractions of the data are ## in the lowest category: ## 165/209=0.001 'drug' ## 54/209=0.001 'alcoh' ## 38/209=0.001 'cond' ## so this will be a fairly challenging problem in any case library(ggplot2) ggplot(results,aes(x=alcoh,y=cond))+stat_sum(aes(size=..n..),alpha=0.7)+ facet_wrap(~drugcat)+theme_bw() library(betareg) ## set phi link to logarithmic ## basic problem (digging through betareg.fit etc.) is ## that initial estimate of phi, based on ## linear model of logit(cond) ~ alcoh + cond, is NEGATIVE ... ## doesn't seem to be any way to override this starting value ## brute force try(gyl-betareg(cond ~ alcoh + drug, data=results, link.phi=log)) ## pick through, debugging ... find starting values used svec - c(-1.6299469,0.8048446,1.7071124,0) gyl-betareg(cond ~ alcoh + drug, data=results, link.phi=log, control=betareg.control(start=svec)) ## would work fine with more generic starting values svec2 - c(qlogis(mean(results$cond)),0,0,0) gyl2-betareg(cond ~ alcoh + drug, data=results, link.phi=log, control=betareg.control(start=svec2)) ## before I got that to work, I also tried this (which ## will be slower and less efficient but is a useful ## alternative library(bbmle) ## define a variant parameterization of the beta distribution with ## m=a/(a+b), phi=(a+b) dbeta2 - function(x,m,phi,log=FALSE) { a - m*phi b - phi*(1-m) dbeta(x,shape1=a,shape2=b,log=log) } m1 - mle2(cond~dbeta2(m=plogis(mu),phi=exp(logphi)), parameters=list(mu~alcoh+drug), data=results, start=list(mu=qlogis(mean(results$cond)),logphi=0)) summary(m1) p1 - profile(m1) plot(p1,show.points=TRUE) confint(p1) confint(m1,method=quad) ## not much difference coef(m1) coef(gyl) On 11-03-13 12:59 PM, Vlatka Matkovic Puljic wrote: Sorry, here is my data (attached). 2011/3/12 Ben Bolker bbol...@gmail.com mailto:bbol...@gmail.com Vlatka Matkovic Puljic v.matkovic.puljic at gmail.com http://gmail.com writes: That was also my first thought. But I guess it has something to do with W and phihat (which I'm struggling to check Again, it would help to post a reproducible example ... hard to debug/diagnose by remote control. If you can't possibly post the data to the list, or put them on a web site somewhere, or randomize them a bit so you're not giving anything away, or find a simulated example that shows the same problem, you could as a last resort send them to me. Ben Bolker __ R-help@r-project.org mailto:R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- ** Vlatka Matkovic Puljic +32/ 485/ 453340 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- ** Vlatka Matkovic Puljic +32/ 485/ 453340 [[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
Re: [R] betareg help
Sorry, here is my data (attached). 2011/3/12 Ben Bolker bbol...@gmail.com Vlatka Matkovic Puljic v.matkovic.puljic at gmail.com writes: That was also my first thought. But I guess it has something to do with W and phihat (which I'm struggling to check Again, it would help to post a reproducible example ... hard to debug/diagnose by remote control. If you can't possibly post the data to the list, or put them on a web site somewhere, or randomize them a bit so you're not giving anything away, or find a simulated example that shows the same problem, you could as a last resort send them to me. Ben Bolker __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- ** Vlatka Matkovic Puljic +32/ 485/ 453340 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] betareg help
http://dl.dropbox.com/u/21595123/Book1.csv 2011/3/13 David Winsemius dwinsem...@comcast.net Nothing came through. You need to read the posting guide. On Mar 13, 2011, at 12:59 PM, Vlatka Matkovic Puljic wrote: Sorry, here is my data (attached). 2011/3/12 Ben Bolker bbol...@gmail.com Vlatka Matkovic Puljic v.matkovic.puljic at gmail.com writes: That was also my first thought. But I guess it has something to do with W and phihat (which I'm struggling to check Again, it would help to post a reproducible example ... hard to debug/diagnose by remote control. If you can't possibly post the data to the list, or put them on a web site somewhere, or randomize them a bit so you're not giving anything away, or find a simulated example that shows the same problem, you could as a last resort send them to me. Ben Bolker David Winsemius, MD West Hartford, CT -- ** Vlatka Matkovic Puljic +32/ 485/ 453340 [[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] betareg help
Dear R users, I'm trying to do betareg on my dataset. Dependent variable is not normally distributed and is proportion (of condom use (0,1)). But I'm having problems: gyl-betareg(cond ~ alcoh + drug, data=results) Error in optim(par = start, fn = loglikfun, gr = gradfun, method = method, : initial value in 'vmmin' is not finite Why is R returning me error in optim()? What should I do? In advance, thank you for advice! -- ** Vlatka Matkovic Puljic +32/ 485/ 453340 [[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] betareg help
Maybe I should include data: results$cond [1] 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 [13] 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 [25] 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 [37] 0.001 0.001 0.020 0.020 0.020 0.030 0.030 0.050 0.050 0.050 0.050 0.050 [49] 0.050 0.050 0.050 0.050 0.060 0.060 0.060 0.070 0.070 0.070 0.070 0.080 [61] 0.090 0.100 0.100 0.100 0.100 0.100 0.100 0.120 0.120 0.130 0.130 0.140 [73] 0.160 0.160 0.170 0.170 0.180 0.180 0.200 0.200 0.210 0.210 0.210 0.220 [85] 0.230 0.250 0.250 0.250 0.250 0.270 0.280 0.300 0.300 0.300 0.300 0.300 [97] 0.310 0.320 0.330 0.340 0.350 0.350 0.350 0.360 0.400 0.430 0.430 0.450 [109] 0.450 0.450 0.450 0.450 0.450 0.450 0.450 0.450 0.450 0.460 0.470 0.470 [121] 0.470 0.470 0.480 0.490 0.490 0.500 0.500 0.500 0.500 0.500 0.500 0.500 [133] 0.500 0.500 0.500 0.520 0.530 0.550 0.550 0.550 0.560 0.600 0.600 0.600 [145] 0.600 0.600 0.620 0.640 0.650 0.650 0.650 0.650 0.660 0.680 0.680 0.680 [157] 0.680 0.700 0.700 0.700 0.700 0.700 0.700 0.710 0.740 0.750 0.750 0.750 [169] 0.750 0.760 0.760 0.770 0.780 0.800 0.800 0.800 0.800 0.800 0.810 0.820 [181] 0.820 0.830 0.830 0.840 0.850 0.850 0.850 0.850 0.860 0.870 0.870 0.870 [193] 0.900 0.900 0.900 0.900 0.900 0.910 0.920 0.920 0.930 0.930 0.950 0.950 [205] 0.950 0.960 0.980 0.980 0.999 2011/3/12 Vlatka Matkovic Puljic v.matkovic.pul...@gmail.com Dear R users, I'm trying to do betareg on my dataset. Dependent variable is not normally distributed and is proportion (of condom use (0,1)). But I'm having problems: gyl-betareg(cond ~ alcoh + drug, data=results) Error in optim(par = start, fn = loglikfun, gr = gradfun, method = method, : initial value in 'vmmin' is not finite Why is R returning me error in optim()? What should I do? In advance, thank you for advice! -- ** Vlatka Matkovic Puljic +32/ 485/ 453340 -- ** Vlatka Matkovic Puljic +32/ 485/ 453340 [[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] betareg help
That was also my first thought. But I guess it has something to do with W and phihat (which I'm struggling to check?) 2011/3/12 Ben Bolker bbol...@gmail.com Vlatka Matkovic Puljic v.matkovic.puljic at gmail.com writes: My first guess would have been that you had zeros and ones in your data, but you seem to have taken care of that. Are you willing to post the whole data set (use dput() to get it in a useful form)? Ben Bolker __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- ** Vlatka Matkovic Puljic +32/ 485/ 453340 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] difference in pairs ?
Well, it should be difference by ID and TIME for q1: something like: for ID 1187 in TIME 1 q1=3 and TIME 2 (for same ID) q1=3 so diff would be 3-3=0 TIME ID q1 1 1187 3 1 1187 3 And I don't know how to make R to find pairs and calculate diff? 2011/2/21 Dennis Murphy djmu...@gmail.com Hi: Assuming dd is the name of your data frame, dd$diff - with(dd, q2 - q1) dd TIME ID q1 q2 diff 11 1187 3 2 -1 21 1706 3 30 31 1741 2 42 42 1187 3 2 -1 52 1706 3 30 62 1741 2 42 is one way to do it. HTH, Dennis On Mon, Feb 21, 2011 at 11:05 AM, Vlatka Matkovic Puljic v.matkovic.pul...@gmail.com wrote: Dear all, I want to perform paired Wilcoxon signed ranks test on my data. I have pairs defined by ID and TIME variables. How can I calculate difference in variables q1, q2 in each pair? TIME ID q1 q2 1 1187 3 2 1 1706 3 3 1 1741 2 4 2 1187 3 2 2 1706 3 3 2 1741 2 4 Please, any clue! :) -- ** Vlatka Matkovic Puljic +32/ 485/ 453340 [[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. -- ** Vlatka Matkovic Puljic +32/ 485/ 453340 [[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] difference in pairs ?
Dear all, I want to perform paired Wilcoxon signed ranks test on my data. I have pairs defined by ID and TIME variables. How can I calculate difference in variables q1, q2 in each pair? TIME ID q1 q2 1 1187 3 2 1 1706 3 3 1 1741 2 4 2 1187 3 2 2 1706 3 3 2 1741 2 4 Please, any clue! :) -- ** Vlatka Matkovic Puljic +32/ 485/ 453340 [[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] compare variables and take largest value?
Thanx a lot! That's it! 2010/4/22 Duncan Murdoch murdoch.dun...@gmail.com On 22/04/2010 11:07 AM, Vlatka Matkovic Puljic wrote: Hi! I have v1, v2 and v3 (data in v. ranges 1 to 4) v1 [1] 1 1 1 1 NaN 1 1 4 1 1 3 v2 [1] 1 1 1 1 1 1 1 1 1 1 1 v3 [1] 1 1 1 1 1 1 1 1 1 1 1 I want R to compare entries in v1,v2,v3 and in v4 to put largest value (NA=0). For example. first line v4 would be 1, but last line v4 would be 3. What function should I use? pmax. Duncan Murdoch __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- ** Vlatka Matkovic Puljic 095/8618 171 [[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] compare variables and take largest value?
Only thing that when I put na.rm=FALSE, R is puting in v4 NA (as is largest) and I want to see NA as zero (0, smallest value) 2010/4/22 Vlatka Matkovic Puljic v.matkovic.pul...@gmail.com Thanx a lot! That's it! 2010/4/22 Duncan Murdoch murdoch.dun...@gmail.com On 22/04/2010 11:07 AM, Vlatka Matkovic Puljic wrote: Hi! I have v1, v2 and v3 (data in v. ranges 1 to 4) v1 [1] 1 1 1 1 NaN 1 1 4 1 1 3 v2 [1] 1 1 1 1 1 1 1 1 1 1 1 v3 [1] 1 1 1 1 1 1 1 1 1 1 1 I want R to compare entries in v1,v2,v3 and in v4 to put largest value (NA=0). For example. first line v4 would be 1, but last line v4 would be 3. What function should I use? pmax. Duncan Murdoch __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- ** Vlatka Matkovic Puljic 095/8618 171 -- ** Vlatka Matkovic Puljic 095/8618 171 [[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] median Q?
Hi, I have dataset n1 and v1 (years). when i ask median(year) [1] NA but if i put summary of dataset n1: summary(n1) R produces median (together with min/max/mean) why it is so? -- ** Vlatka Matkovic Puljic 095/8618 171 [[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] calculating age from year of birth and date?
Hi, I have v1 (date of test) and v2 (year of birth). v1 v2 15.5.2008 88 18.6.200954 I want R to use only year in v1; and v2 see as 1988 and calculate age in v3. any ideas how to do that? -- ** Vlatka MatkoviÄ PuljiÄ 095/8618 171 [[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] median Q?
Thank you! 2010/4/21 Joshua Wiley jwiley.ps...@gmail.com Hello, summary() removes NAs by default. You can get the same results using median(year, na.rm=TRUE) see ?median HTH, Josh On Wed, Apr 21, 2010 at 6:40 AM, Vlatka Matkovic Puljic v.matkovic.pul...@gmail.com wrote: Hi, I have dataset n1 and v1 (years). when i ask median(year) [1] NA but if i put summary of dataset n1: summary(n1) R produces median (together with min/max/mean) why it is so? -- ** Vlatka Matkovic Puljic 095/8618 171 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Joshua Wiley Senior in Psychology University of California, Riverside http://www.joshuawiley.com/ -- ** Vlatka Matkovic Puljic 095/8618 171 [[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] calculating age from year of birth and date?
Hm, well.. I can do it in xls but I wanted to see if there is some (simple) way to do this in R. 2010/4/21 Dimitri Liakhovitski ld7...@gmail.com And what is the format for year of birth in Excel? Can't you use concatenate to concatenate a new column (that contains only 19) and the column that contains a 2-digit year of birth? And then change all those who are born (in the new column) before 1910 to XX+100 - which will give you 2009 (or 2001), etc.? On Wed, Apr 21, 2010 at 9:57 AM, Vlatka Matkovic Puljic v.matkovic.pul...@gmail.com wrote: Hi, well...when imported to R (from xls file) it transform to number print(D1) [1] 39804 39527 39917 39860 39489 ??? but in xls is dd/mm/ and year of birth is as yy those born after 1999 are removed Dana 21. travnja 2010. 15:52 Dimitri Liakhovitski ld7...@gmail.com je napisao/la: In what format are you dates? Are they always like this: two digits for day, then one or two digits for month, and four digits for the year? Also - are any of your people born after 1999? Dimitri 2010/4/21 Vlatka MatkoviÄ PuljiÄ vlatk...@gmail.com: Hi, I have v1 (date of test) and v2 (year of birth). v1 v2 15.5.2008 88 18.6.200954 I want R to use only year in v1; and v2 see as 1988 and calculate age in v3. any ideas how to do that? -- ** Vlatka MatkoviÄ PuljiÄ 095/8618 171 [[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. -- Dimitri Liakhovitski Ninah.com dimitri.liakhovit...@ninah.com -- ** Vlatka Matkovic Puljic 095/8618 171 -- Dimitri Liakhovitski Ninah.com dimitri.liakhovit...@ninah.com -- ** Vlatka Matkovic Puljic 095/8618 171 [[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] median Q?
I have additional Q: v1 is gender (M=1 and F=2) v2 is age I want R to calculate median only for M (1), but my comand is not good :) while(v1=1){median(v2,na.rm=TRUE)} Error: unexpected '=' in while(Q2= 2010/4/21 Vlatka Matkovic Puljic v.matkovic.pul...@gmail.com Thank you! 2010/4/21 Joshua Wiley jwiley.ps...@gmail.com Hello, summary() removes NAs by default. You can get the same results using median(year, na.rm=TRUE) see ?median HTH, Josh On Wed, Apr 21, 2010 at 6:40 AM, Vlatka Matkovic Puljic v.matkovic.pul...@gmail.com wrote: Hi, I have dataset n1 and v1 (years). when i ask median(year) [1] NA but if i put summary of dataset n1: summary(n1) R produces median (together with min/max/mean) why it is so? -- ** Vlatka Matkovic Puljic 095/8618 171 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Joshua Wiley Senior in Psychology University of California, Riverside http://www.joshuawiley.com/ -- ** Vlatka Matkovic Puljic 095/8618 171 -- ** Vlatka Matkovic Puljic 095/8618 171 [[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] median Q?
Perfect! Thanx a lot! :) 2010/4/21 Joshua Wiley jwiley.ps...@gmail.com Look at ?by for example by(data=v2, INDICES=v1, FUN=median, na.rm=TRUE) This will calculate the median of v2 (age) for each level of the indices v1 (in your case M and F). If you are only interested the median for a single level, Mohamed's solution is simpler. Josh On Wed, Apr 21, 2010 at 8:10 AM, Vlatka Matkovic Puljic v.matkovic.pul...@gmail.com wrote: I have additional Q: v1 is gender (M=1 and F=2) v2 is age I want R to calculate median only for M (1), but my comand is not good :) while(v1=1){median(v2,na.rm=TRUE)} Error: unexpected '=' in while(Q2= 2010/4/21 Vlatka Matkovic Puljic v.matkovic.pul...@gmail.com Thank you! 2010/4/21 Joshua Wiley jwiley.ps...@gmail.com Hello, summary() removes NAs by default. You can get the same results using median(year, na.rm=TRUE) see ?median HTH, Josh On Wed, Apr 21, 2010 at 6:40 AM, Vlatka Matkovic Puljic v.matkovic.pul...@gmail.com wrote: Hi, I have dataset n1 and v1 (years). when i ask median(year) [1] NA but if i put summary of dataset n1: summary(n1) R produces median (together with min/max/mean) why it is so? -- ** Vlatka Matkovic Puljic 095/8618 171 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Joshua Wiley Senior in Psychology University of California, Riverside http://www.joshuawiley.com/ -- ** Vlatka Matkovic Puljic 095/8618 171 -- ** Vlatka Matkovic Puljic 095/8618 171 -- Joshua Wiley Senior in Psychology University of California, Riverside http://www.joshuawiley.com/ -- ** Vlatka Matkovic Puljic 095/8618 171 [[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] recoding variables-recode not working
Dear, my variable is numerical indicating how many times smb done test summary(Q12) Min. 1st Qu. MedianMean 3rd Qu.Max.NA's 0. 0. 0. 0.7989 1. 30. 66. I want to change this to categories-- 0=none testing; 1:30=done testing (NA excluded) John, *cut(Q12, breaks=c(0,1,30), include.lowest=TRUE, labels=c('0', '1 and more'))* is working. Only, when I want frequencies (for ex.) for this new 2 categories, R is returning me freq for all 0 to 30. 2010/4/7 John Fox j...@mcmaster.ca Dear Vlatka, It's impossible to know what the problem is without knowing something about your data, which you didn't tell us either in this message or your subsequent one. The recode command should work: (x - c(rep(0, 5), sample(1:30, 5, replace=TRUE))) [1] 0 0 0 0 0 17 27 19 19 2 recode(x, 0='A'; 1:30='B') [1] A A A A A B B B B B The cut command requires include.lowest=TRUE and it helps to spell the labels argument correctly: cut(x, breaks=c(0,1,30), include.lowest=TRUE, labels=c('0', '1 and more')) [1] 0 0 0 0 0 1 and more [7] 1 and more 1 and more 1 and more 1 and more Levels: 0 1 and more I hope this helps, John John Fox Senator William McMaster Professor of Social Statistics Department of Sociology McMaster University Hamilton, Ontario, Canada web: socserv.mcmaster.ca/jfox -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Vlatka Matkovic Puljic Sent: April-07-10 1:31 PM To: r-help@r-project.org Subject: [R] recoding variables-recode not working Hi, I have numerical variable that I want to recode into categories '0' and '1 and more' and do analysis with that data. I have tried various of possibilities to do so, but I am sucked and nothing is working. recode(Q12, 0='A';1:30='B') cut(Q12, breaks=c(0,1,30), lables=c('0', '1 and more')) cat(Q12, 0=0;1-33=1) What should I do to make it right? -- ** Vlatka [[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. -- ** Vlatka MatkoviÄ PuljiÄ 095/8618 171 [[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] recoding variables-recode not working
Now I have :) Thanx a lot! 2010/4/8 David Winsemius dwinsem...@comcast.net On Apr 8, 2010, at 8:13 AM, Vlatka Matkovic Puljic wrote: Dear, my variable is numerical indicating how many times smb done test summary(Q12) Min. 1st Qu. MedianMean 3rd Qu.Max.NA's 0. 0. 0. 0.7989 1. 30. 66. I want to change this to categories-- 0=none testing; 1:30=done testing (NA excluded) John, *cut(Q12, breaks=c(0,1,30), include.lowest=TRUE, labels=c('0', '1 and more'))* is working. Only, when I want frequencies (for ex.) for this new 2 categories, R is returning me freq for all 0 to 30. It what fashion are you telling R that you want frequencies? Did you assign the output of cut(Q12, ...) to a new variable? -- David. 2010/4/7 John Fox j...@mcmaster.ca Dear Vlatka, It's impossible to know what the problem is without knowing something about your data, which you didn't tell us either in this message or your subsequent one. The recode command should work: (x - c(rep(0, 5), sample(1:30, 5, replace=TRUE))) [1] 0 0 0 0 0 17 27 19 19 2 recode(x, 0='A'; 1:30='B') [1] A A A A A B B B B B The cut command requires include.lowest=TRUE and it helps to spell the labels argument correctly: cut(x, breaks=c(0,1,30), include.lowest=TRUE, labels=c('0', '1 and more')) [1] 0 0 0 0 0 1 and more [7] 1 and more 1 and more 1 and more 1 and more Levels: 0 1 and more I hope this helps, John John Fox Senator William McMaster Professor of Social Statistics Department of Sociology McMaster University Hamilton, Ontario, Canada web: socserv.mcmaster.ca/jfox -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org ] On Behalf Of Vlatka Matkovic Puljic Sent: April-07-10 1:31 PM To: r-help@r-project.org Subject: [R] recoding variables-recode not working Hi, I have numerical variable that I want to recode into categories '0' and '1 and more' and do analysis with that data. I have tried various of possibilities to do so, but I am sucked and nothing is working. recode(Q12, 0='A';1:30='B') cut(Q12, breaks=c(0,1,30), lables=c('0', '1 and more')) cat(Q12, 0=0;1-33=1) What should I do to make it right? -- ** Vlatka [[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. -- ** Vlatka Matkoviæ Puljiæ 095/8618 171 [[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. David Winsemius, MD West Hartford, CT -- ** Vlatka Matkoviæ Puljiæ 095/8618 171 [[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] recoding variables-recode not working
Hi, I have numerical variable that I want to recode into categories '0' and '1 and more' and do analysis with that data. I have tried various of possibilities to do so, but I am sucked and nothing is working. recode(Q12, 0='A';1:30='B') cut(Q12, breaks=c(0,1,30), lables=c('0', '1 and more')) cat(Q12, 0=0;1-33=1) What should I do to make it right? -- ** Vlatka [[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] recoding variables-recode not working
atomic [1:1578] 0 0 0 0 0 0 0 0 0 0 ... - attr(*, levels)= chr 0=0,1:33=1 2010/4/7 David Winsemius dwinsem...@comcast.net On Apr 7, 2010, at 1:31 PM, Vlatka Matkovic Puljic wrote: Hi, I have numerical variable that I want to recode into categories '0' and '1 and more' and do analysis with that data. I have tried various of possibilities to do so, but I am sucked and nothing is working. recode(Q12, 0='A';1:30='B') cut(Q12, breaks=c(0,1,30), lables=c('0', '1 and more')) cat(Q12, 0=0;1-33=1) What should I do to make it right? What does str(Q12) return? -- ** Vlatka [[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. David Winsemius, MD West Hartford, CT -- ** Vlatka Matkoviæ Puljiæ 095/8618 171 [[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.