[R] problems perfroming the vuong test
Dear All, I am using the function vuong of the package pscl to compare 2 non nested glm models with a numeric response. I did the following m1-glm(y ~x ,data=xxx) m2-glm(y ~z , data=xxx) When calling the vuong function I get the following message: vuong(m1,m2) Error in predprob.glm(m1) : your object of class glm is unsupported by predprob.glmyour object of class lm is unsupported by predprob.glm My guess is that this function does not support numeric response! How can I solve the problem? Any help will be really appreciated. mirko __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] gamma distribution don't allow negative value in GLMs?
Dear friends, when i use glm() to fit my data, i use glm(formula = snail ~ vegtype + mhveg + humidity + elevation + soiltem, *family = Gamma(link = inverse),* data =a,)) It shows: error in eval(expr, envir, enclos) : *gamma distribution don't allow negative value*. But i use result-glm(formula = snail ~ vegtype + mhveg + humidity + elevation + soiltem, family = poisson, data =a) #this works In fact , there isn't any negative value in my dataset, who can tell me the reason? Thanks very much! I copy my data here so you can check it: vegtype mhveg humidity soiltem elevation snail 1 diluo 35.0 0.27985121.1 low 162 2 diluo 25.0 0.31609223.1 low 113 3 yuhao 35.0 0.29723821.7 low 105 4 huanghuacai 1.5 0.31068723.1 low 5 5 huanghuacai 2.0 0.26786828.3 low 1 6 yuhao 25.0 0.29013521.9 low10 7 huanghuacai 1.0 0.28520727.7 low 6 8 huanghuacai 2.0 0.25297328.3 low 1 9 huanghuacai 1.5 0.2728.1 low 1 10 huanghuacai 2.5 0.3029.1 low 1 11 huanghuacai 2.0 0.29615429.1 low 0 12 huanghuacai 2.0 0.30287427.5 low 3 13 huanghuacai 1.5 0.30149928.9 low 0 14 huanghuacai 3.0 0.29151330.3 low 1 15 huanghuacai 1.0 0.27343831.1 low 3 16 huanghuacai 1.5 0.29011627.9 low19 17 huanghuacai 2.5 0.19893231.9 low 0 18 huanghuacai 2.0 0.3930.5 high 4 19 huanghuacai 2.5 0.28259530.7 high 0 20 huanghuacai 1.0 0.26609724.7 high14 21yuhao 30.0 0.24051626.9 high51 22yuhao 35.0 0.22754126.7 high84 23yuhao 20.0 0.25283328.3 low30 24diluo 40.0 0.30303027.9 low91 25hucao 80.0 0.30386724.5 low 114 26diluo 25.0 0.33494826.7 low 115 27hucao 60.0 0.30689726.5 low23 28hucao 75.0 0.31446525.7 low43 29yuhao 30.0 0.25178326.1 low77 30diluo 10.0 0.2826.1 low62 31yuhao 25.0 0.29171626.1 low78 32hucao 90.0 0.28880024.5 low35 33diluo 25.0 0.33783026.3 high75 34yuhao 13.0 0.29659927.7 high23 35hucao 70.0 0.27949826.3 high 116 36diluo 3.0 0.28148128.1 high25 37hucao 70.0 0.29600023.7 high83 38diluo 10.0 0.27266227.7 low56 39hucao 70.0 0.28979625.3 high 112 40diluo 5.0 0.33971627.9 high84 41yuhao 35.0 0.23142724.9 high88 42hucao 80.0 0.27381024.1 high 134 43yuhao 40.0 0.27278925.1 high53 44yuhao 45.0 0.22603625.1 high88 45yuhao 55.0 0.28549523.9 high76 46hucao 80.0 0.25218523.9 high 106 47diluo 15.0 0.28993324.5 high 194 48hucao 95.0 0.26175623.1 high35 49hucao 55.0 0.23981924.7 high21 50hucao 75.0 0.25430723.9 high41 51 huanghuacai 1.0 0.28643223.7 low18 52 huanghuacai 2.0 0.30134223.1 low 2 53 huanghuacai 2.0 0.36956523.3 low 5 54 huanghuacai 1.5 0.24583324.3 low 4 55 huanghuacai 1.0 0.31567924.1 low 4 56 huanghuacai 2.5 0.29612423.7 low 4 57 huanghuacai 2.0 0.31266725.7 low 3 58 huanghuacai 3.0 0.30087025.7 low 0 59 huanghuacai 2.0 0.30374326.5 low 2 60 huanghuacai 1.0 0.26979925.3 low 7 61hucao 75.0 0.28125022.5 low14 62yuhao 35.0 0.35035023.3 low63 63hucao 65.0 0.30454522.7 low17 64diluo 7.0 0.31005624.9 low45 65hucao 80.0 0.28800022.9 low27 66hucao 80.0 0.28421122.7 low46 67diluo 25.0 0.28137923.5 low 161 68hucao 80.0 0.29053323.3 low 117 69yuhao 27.0 0.31656824.1 low 106 70yuhao 28.0 0.28515625.1 low82 71yuhao 30.0 0.2724.5 low55 72hucao 85.0 0.29034523.9 low54 73yuhao 35.0 0.31578924.1 low81 74diluo 15.0 0.28659828.3 low 102 75yuhao 45.0 0.31421124.1 low85 76yuhao 25.0 0.26879425.1 low63 77hucao 80.0 0.27569123.9 low59 78hucao 100.0 0.31661424.1 low46 79yuhao 40.0 0.33668325.5 low70 80diluo 20.0 0.27087426.1 high 167 81
Re: [R] combinatorics
Robin Hankin [EMAIL PROTECTED] writes: Hi How do I generate all ways of ordering sets of indistinguishable items? suppose I have two A's, two B's and a C. Then I want AABBC AABCB AACBC ABABC . . .snip... BBAAC . . .snip... CBBAA [there are 5!/(2!*2!) = 30 arrangements. Note AABBC != BBAAC] How do I do this? See Knuth, The Art of Computer Programming Vol 4, Fascicle 3 and 4. Jens __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] gamma distribution don't allow negative value in GLMs?
I think the 0 values for snail are hurting you. Kees On Sunday 15 October 2006 13:10, zhijie zhang wrote: Dear friends, when i use glm() to fit my data, i use glm(formula = snail ~ vegtype + mhveg + humidity + elevation + soiltem, *family = Gamma(link = inverse),* data =a,)) It shows: error in eval(expr, envir, enclos) : *gamma distribution don't allow negative value*. But i use result-glm(formula = snail ~ vegtype + mhveg + humidity + elevation + soiltem, family = poisson, data =a) #this works In fact , there isn't any negative value in my dataset, who can tell me the reason? Thanks very much! I copy my data here so you can check it: vegtype mhveg humidity soiltem elevation snail 1 diluo 35.0 0.27985121.1 low 162 2 diluo 25.0 0.31609223.1 low 113 3 yuhao 35.0 0.29723821.7 low 105 4 huanghuacai 1.5 0.31068723.1 low 5 5 huanghuacai 2.0 0.26786828.3 low 1 6 yuhao 25.0 0.29013521.9 low10 7 huanghuacai 1.0 0.28520727.7 low 6 8 huanghuacai 2.0 0.25297328.3 low 1 9 huanghuacai 1.5 0.2728.1 low 1 10 huanghuacai 2.5 0.3029.1 low 1 11 huanghuacai 2.0 0.29615429.1 low 0 12 huanghuacai 2.0 0.30287427.5 low 3 13 huanghuacai 1.5 0.30149928.9 low 0 14 huanghuacai 3.0 0.29151330.3 low 1 15 huanghuacai 1.0 0.27343831.1 low 3 16 huanghuacai 1.5 0.29011627.9 low19 17 huanghuacai 2.5 0.19893231.9 low 0 18 huanghuacai 2.0 0.3930.5 high 4 19 huanghuacai 2.5 0.28259530.7 high 0 20 huanghuacai 1.0 0.26609724.7 high14 21yuhao 30.0 0.24051626.9 high51 22yuhao 35.0 0.22754126.7 high84 23yuhao 20.0 0.25283328.3 low30 24diluo 40.0 0.30303027.9 low91 25hucao 80.0 0.30386724.5 low 114 26diluo 25.0 0.33494826.7 low 115 27hucao 60.0 0.30689726.5 low23 28hucao 75.0 0.31446525.7 low43 29yuhao 30.0 0.25178326.1 low77 30diluo 10.0 0.2826.1 low62 31yuhao 25.0 0.29171626.1 low78 32hucao 90.0 0.28880024.5 low35 33diluo 25.0 0.33783026.3 high75 34yuhao 13.0 0.29659927.7 high23 35hucao 70.0 0.27949826.3 high 116 36diluo 3.0 0.28148128.1 high25 37hucao 70.0 0.29600023.7 high83 38diluo 10.0 0.27266227.7 low56 39hucao 70.0 0.28979625.3 high 112 40diluo 5.0 0.33971627.9 high84 41yuhao 35.0 0.23142724.9 high88 42hucao 80.0 0.27381024.1 high 134 43yuhao 40.0 0.27278925.1 high53 44yuhao 45.0 0.22603625.1 high88 45yuhao 55.0 0.28549523.9 high76 46hucao 80.0 0.25218523.9 high 106 47diluo 15.0 0.28993324.5 high 194 48hucao 95.0 0.26175623.1 high35 49hucao 55.0 0.23981924.7 high21 50hucao 75.0 0.25430723.9 high41 51 huanghuacai 1.0 0.28643223.7 low18 52 huanghuacai 2.0 0.30134223.1 low 2 53 huanghuacai 2.0 0.36956523.3 low 5 54 huanghuacai 1.5 0.24583324.3 low 4 55 huanghuacai 1.0 0.31567924.1 low 4 56 huanghuacai 2.5 0.29612423.7 low 4 57 huanghuacai 2.0 0.31266725.7 low 3 58 huanghuacai 3.0 0.30087025.7 low 0 59 huanghuacai 2.0 0.30374326.5 low 2 60 huanghuacai 1.0 0.26979925.3 low 7 61hucao 75.0 0.28125022.5 low14 62yuhao 35.0 0.35035023.3 low63 63hucao 65.0 0.30454522.7 low17 64diluo 7.0 0.31005624.9 low45 65hucao 80.0 0.28800022.9 low27 66hucao 80.0 0.28421122.7 low46 67diluo 25.0 0.28137923.5 low 161 68hucao 80.0 0.29053323.3 low 117 69yuhao 27.0 0.31656824.1 low 106 70yuhao 28.0 0.28515625.1 low82 71yuhao 30.0 0.2724.5 low55 72hucao 85.0 0.29034523.9 low54 73yuhao 35.0 0.31578924.1 low81 74diluo 15.0 0.28659828.3 low 102 75yuhao 45.0 0.31421124.1 low85 76yuhao 25.0 0.26879425.1 low63 77hucao 80.0
[R] executing strings
I'd like to excute character strings such as z-plot( objects()[1]; eval(z) and viola I'd have a plot of my first dataframe in the first frame. Unfortunately this approach no longer works. Help? Lloyd L [EMAIL PROTECTED] [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Random efftects distribution for nlme
On 10/14/06, Cam [EMAIL PROTECTED] wrote: For nonlinear mixed effects models, can you specify the use of mass points to approximate the random effects distribution? How? If you mean in nonlinear mixed effects models fit by the nlme function the answer is no. The nlme function in the nlme package only allows for a Gaussian distribution of the random effects. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] executing strings
Try this: eval(parse(text = pi + 3)) On 10/14/06, Lloyd Lubet [EMAIL PROTECTED] wrote: I'd like to excute character strings such as z-plot( objects()[1]; eval(z) and viola I'd have a plot of my first dataframe in the first frame. Unfortunately this approach no longer works. Help? Lloyd L [EMAIL PROTECTED] [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] glm cannot find valid starting values
Em Sábado 14 Outubro 2006 11:15, Gabor Grothendieck escreveu: Try using OLS starting values: glm(Y~X,family=gaussian(link=log), start = coef(lm(Y~X))) Ok, using a starting value the model work. But, sometimes ago it work without any starting values. Why now I need a starting values? Thanks Ronaldo -- Nunca faça previsões, especialmente sobre o futuro. -- Samuel Goldwyn -- Prof. Ronaldo Reis Júnior | .''`. UNIMONTES/Depto. Biologia Geral/Lab. Ecologia Evolutiva | : :' : Campus Universitário Prof. Darcy Ribeiro, Vila Mauricéia | `. `'` CP: 126, CEP: 39401-089, Montes Claros - MG - Brasil | `- Fone: (38) 3229-8190 | [EMAIL PROTECTED] | ICQ#: 5692561 | LinuxUser#: 205366 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Hide line ends behind unfilled circles?
Dear r-helpers, xx - c(0.000, 0.210, 0.714, 0.514, 1.000, 0.190, 0.590, 0.152) yy - c(0.000, 0.265, 0.256, 0.521, 0.538, 0.761, 0.821, 1.000) aa - c(19, 19, 19, 21, 19, 21, 21, 21) x0 - xx[c(1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 6, 6, 7)] y0 - yy[c(1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 6, 6, 7)] x1 - xx[c(2, 3, 3, 4, 6, 4, 5, 5, 6, 7, 7, 7, 8, 8)] y1 - yy[c(2, 3, 3, 4, 6, 4, 5, 5, 6, 7, 7, 7, 8, 8)] plot(yy ~ xx, pch = aa, cex = 3) segments(x0, y0, x1, y1) Can anyone suggest a way of insuring that the lines are hidden behind the unfilled circles? _ Professor Michael Kubovy University of Virginia Department of Psychology USPS: P.O.Box 400400Charlottesville, VA 22904-4400 Parcels:Room 102Gilmer Hall McCormick RoadCharlottesville, VA 22903 Office:B011+1-434-982-4729 Lab:B019+1-434-982-4751 Fax:+1-434-982-4766 WWW:http://www.people.virginia.edu/~mk9y/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] split-plot analysis with lme()
The problem in your example is that 'lme' doesn't know how to handle the Variety*nitro interaction when all 12 combinations are not present. The error message singularity in backsolve means that with data for only 11 combinations, which is what you have in your example, you can only estimate 11 linearly independent fixed-effect coefficients, not the 12 required by this model: 1 for intercept + (3-1) for Variety + (4-1) for nitro + (3-1)*(4-1) for Variety*nitro = 12. Since 'nitro' is a fixed effect only, you can get what you want by keeping it as a numeric factor and manually specifying the (at most 5, not 6) interaction contrasts you want, something like the following: fit2. - lme(yield ~ Variety+nitro+I(nitro^2)+I(nitro^3) +Variety:(nitro+I(nitro^2)), data=Oats, random=~1|Block/Variety, subset=!(Variety == Golden Rain nitro == 0)) NOTE: This gives us 4 degrees of freedom for the interaction. With all the data, we can estimate 6. Therefore, there should be some way to get 5, but so far I haven't figured out an easy way to do that. Perhaps someone else will enlighten us both. Even without a method for estimating an interaction term with 5 degrees of freedom, I hope I've at least answered your basic question. Best Wishes, Spencer Graves i.m.s.white wrote: Dear R-help, Why can't lme cope with an incomplete whole plot when analysing a split-plot experiment? For example: R : Copyright 2006, The R Foundation for Statistical Computing Version 2.3.1 (2006-06-01) library(nlme) attach(Oats) nitro - ordered(nitro) fit - lme(yield ~ Variety*nitro, random=~1|Block/Variety) anova(fit) numDF denDF F-value p-value (Intercept) 145 245.14333 .0001 Variety 210 1.48534 0.2724 nitro 345 37.68560 .0001 Variety:nitro 645 0.30282 0.9322 # Excellent! However --- fit2 - lme(yield ~ Variety*nitro, random=~1|Block/Variety, subset= + !(Variety == Golden Rain nitro == 0)) Error in MEEM(object, conLin, control$niterEM) : Singularity in backsolve at level 0, block 1 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Hide line ends behind unfilled circles?
Michael Kubovy wrote: Dear r-helpers, xx - c(0.000, 0.210, 0.714, 0.514, 1.000, 0.190, 0.590, 0.152) yy - c(0.000, 0.265, 0.256, 0.521, 0.538, 0.761, 0.821, 1.000) aa - c(19, 19, 19, 21, 19, 21, 21, 21) x0 - xx[c(1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 6, 6, 7)] y0 - yy[c(1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 6, 6, 7)] x1 - xx[c(2, 3, 3, 4, 6, 4, 5, 5, 6, 7, 7, 7, 8, 8)] y1 - yy[c(2, 3, 3, 4, 6, 4, 5, 5, 6, 7, 7, 7, 8, 8)] plot(yy ~ xx, pch = aa, cex = 3) segments(x0, y0, x1, y1) Can anyone suggest a way of insuring that the lines are hidden behind the unfilled circles? For example, draw the cirles filled (with white background): plot(yy ~ xx, pch = aa, cex = 3, type = n) segments(x0, y0, x1, y1) points(yy ~ xx, pch = aa, cex=3, bg = white) Uwe Ligges _ Professor Michael Kubovy University of Virginia Department of Psychology USPS: P.O.Box 400400Charlottesville, VA 22904-4400 Parcels:Room 102Gilmer Hall McCormick RoadCharlottesville, VA 22903 Office:B011+1-434-982-4729 Lab:B019+1-434-982-4751 Fax:+1-434-982-4766 WWW:http://www.people.virginia.edu/~mk9y/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] commands overwritten
Yuval Sapir wrote: Hi, I am often re-use commands in R by using the arrows to retype them and modify for the new needs. This way, I don't need to write commands again, and I can modify the command by inserting/deleting/re-writing the command. Today, suddenly, re-writing change the commands as the insert function does in Word. for example - if I am standing in the middle of model-glm(seeds~RIL+height); summary(model) trying to re-write it to model-glm(seeds~(RIL+height)^2); summary(model) I get instead model-glm(seeds~(IL+height)^2)ummary(model) How can I get rid of this awkwardiness? (YES, I tried already to hit insert key) Which OS is this? If Windows: Ctrl + O Uwe Ligges Thanks, Yuval __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Hide line ends behind unfilled circles?
Michael Kubovy [EMAIL PROTECTED] writes: Dear r-helpers, xx - c(0.000, 0.210, 0.714, 0.514, 1.000, 0.190, 0.590, 0.152) yy - c(0.000, 0.265, 0.256, 0.521, 0.538, 0.761, 0.821, 1.000) aa - c(19, 19, 19, 21, 19, 21, 21, 21) x0 - xx[c(1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 6, 6, 7)] y0 - yy[c(1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 6, 6, 7)] x1 - xx[c(2, 3, 3, 4, 6, 4, 5, 5, 6, 7, 7, 7, 8, 8)] y1 - yy[c(2, 3, 3, 4, 6, 4, 5, 5, 6, 7, 7, 7, 8, 8)] plot(yy ~ xx, pch = aa, cex = 3) segments(x0, y0, x1, y1) Can anyone suggest a way of insuring that the lines are hidden behind the unfilled circles? plot(yy ~ xx, type = n) segments(x0, y0, x1, y1) points(yy ~ xx, pch = aa, cex = 3, bg=white) -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Hide line ends behind unfilled circles?
On Sun, 2006-10-15 at 11:21 -0400, Michael Kubovy wrote: Dear r-helpers, xx - c(0.000, 0.210, 0.714, 0.514, 1.000, 0.190, 0.590, 0.152) yy - c(0.000, 0.265, 0.256, 0.521, 0.538, 0.761, 0.821, 1.000) aa - c(19, 19, 19, 21, 19, 21, 21, 21) x0 - xx[c(1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 6, 6, 7)] y0 - yy[c(1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 6, 6, 7)] x1 - xx[c(2, 3, 3, 4, 6, 4, 5, 5, 6, 7, 7, 7, 8, 8)] y1 - yy[c(2, 3, 3, 4, 6, 4, 5, 5, 6, 7, 7, 7, 8, 8)] plot(yy ~ xx, pch = aa, cex = 3) segments(x0, y0, x1, y1) Can anyone suggest a way of insuring that the lines are hidden behind the unfilled circles? Try this: # Set up the plot region plot(yy ~ xx, type = n) # Draw the segments first segments(x0, y0, x1, y1) # Set the circle (pch) background colors col - c(rep(black, 3), white, black, rep(white, 3)) # Now draw the circles over the line intersections points(xx, yy, pch = 21, bg = col, cex = 3) The unfilled circles in this case are actually solid white or black. So instead of altering the point character (pch), we alter the colors. HTH, Marc Schwartz __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] extracting nested variances from lme4 model
On 10/14/06, Spencer Graves [EMAIL PROTECTED] wrote: You want to estimate a different 'cs' variance for each level of 'trth', as indicated by the following summary from your 'fake data set': tapply(dtf$x, dtf$trth, sd) FALSE TRUE 1.532251 8.378206 Since var(x[trth]) var(x[!trth]), I thought that the following should produce this: mod1.-lmer( x ~ (1|rtr)+ (trth|cs) , data=dtf) Unfortunately, this generates an error for me: Error in lmer(x ~ (1 | rtr) + (trth | cs), data = dtf) : the leading minor of order 1 is not positive definite I do not understand this error, and I don't have other ideas about how to get around it using 'lmer'. Hmm. It's an interesting problem. I can tell you where the error message originates but I can't yet tell you why. The error message originates when Cholesky decompositions of the variance-covariance matrices for the random effects are generated at the initial estimates. It looks as if one of the model matrices is not being formed correctly for some reason. I will need to do more debugging. It seems to me that 'lme' in library(nlme) should also be able to handle this problem, but 'lme' is designed for nested factors, and your 'rtr' and 'cs' are crossed. If it were my problem, I think I'd study ch. 4 and possibly ch. 5 of Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer) on 'pdMat' and 'corrStruct' classes, because I believe those may support models with crossed random effects like in your example AND might also support estimating different variance components for different levels of a fixed effect like 'trth' in your example. If we are lucky, someone else might educate us both. I'm sorry I couldn't answer your question, but I hope these comments might help in some way. Spencer Graves Frank Samuelson wrote: I have a model: mod1-lmer( x ~ (1|rtr)+ trth/(1|cs) , data=dtf) # Here, cs and rtr are crossed random effects. cs 1-5 are of type TRUE, cs 6-10 are of type FALSE, so cs is nested in trth, which is fixed. So for cs I should get a fit for 1-5 and 6-10. This appears to be the case from the random effects: mean( ranef(mod1)$cs[[1]][1:5] ) [1] -2.498002e-16 var( ranef(mod1)$cs[[1]][1:5] ) [1] 23.53083 mean( ranef(mod1)$cs[[1]][6:10] ) [1] 2.706169e-16 var( ranef(mod1)$cs[[1]][6:10] ) [1] 1.001065 However VarCorr(mod1) gives me only a single variance estimate for the cases. How do I get the other variance estimate? Thanks for any help. -Frank My fake data set: --- data: $frame x rtr trth cs 1 -4.4964808 1 TRUE 1 2 -0.6297254 1 TRUE 2 31.8062857 1 TRUE 3 42.7273275 1 TRUE 4 5 -0.6297254 1 TRUE 5 6 -1.3331767 1 FALSE 6 7 -0.3055796 1 FALSE 7 81.3331767 1 FALSE 8 90.1574314 1 FALSE 9 10 -0.1574314 1 FALSE 10 11 -3.0958803 2 TRUE 1 12 12.4601615 2 TRUE 2 13 -5.2363532 2 TRUE 3 14 1.5419810 2 TRUE 4 15 7.7071005 2 TRUE 5 16 -0.3854953 2 FALSE 6 17 0.7673098 2 FALSE 7 18 2.9485984 2 FALSE 8 19 0.3854953 2 FALSE 9 20 -0.3716558 2 FALSE 10 21 -6.4139940 3 TRUE 1 22 -3.8162344 3 TRUE 2 23 -11.0518554 3 TRUE 3 24 2.7462631 3 TRUE 4 25 16.2543990 3 TRUE 5 26 -1.7353668 3 FALSE 6 27 1.3521369 3 FALSE 7 28 1.3521369 3 FALSE 8 29 0.2054067 3 FALSE 9 30 -1.7446691 3 FALSE 10 31 -6.7468356 4 TRUE 1 32 -11.3228243 4 TRUE 2 33 -18.0088554 4 TRUE 3 34 4.6264891 4 TRUE 4 35 9.2973991 4 TRUE 5 36 -4.1122576 4 FALSE 6 37 -0.5091994 4 FALSE 7 38 1.2787085 4 FALSE 8 39 -1.6867089 4 FALSE 9 40 -0.5091994 4 FALSE 10 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help 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
On Sat, 2006-10-14 at 23:52 -0400, Mohsen Jafarikia wrote:
Hello everyone,
I have the following program to draw a barplot.
MP-read.table(file='AR.out')
names(MP)-c('BN','BL','LR','Q')
Graph- barplot(MP$LR, main='LR Value', col='orange', border='black', space=
0.05, width=(MP$BL), xlab='Length', ylab='LR each')
axis(1, at=Graph, sprintf('%0.2f',MP$BL))
mtext(1, at=Graph, text=sprintf('%0.2f',MP$Q), line=2)
mtext(1, at=par('usr')[1], text='BL', line=1)
mtext(1, at=par('usr')[1], text='Var', line=2)
abline(h=3.841,col='blue')
abline(h=6.635,col='red')
I have two more questions about the graph that I have:
1) I want to write the 'Q' only when it is not equal to zero.
2) I would like to change the bars when 'Q' is not zero. For example, from
orange to green.
I would appreciate your input to this question.
Thanks,
Mohsen
It would be helpful to have the data that you are working with so that
we can provide each other a working example.
However, here are some hints:
1.
mtext(1, at = Graph,
text = ifelse(MP$Q != 0, sprintf('%0.2f',MP$Q), ),
line=2)
2.
cols - ifelse(MP$Q != 0, orange, green)
barplot(..., col = cols, ...)
See ?ifelse
HTH,
Marc Schwartz
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] MASS need R = 2.4.0!!!
Hi, I have R 2.3.1 installed by debian páckage. I install only the base and recommended R packages from Debian source, all others packages I install from the R source at CRAN. But, I try to use library MASS, but I received this message: library(MASS) Error: This is R 2.3.1, package 'MASS' needs = 2.4.0 What is the problem? Inte Ronaldo -- It ain't over until it's over. -- Casey Stengel -- Prof. Ronaldo Reis Júnior | .''`. UNIMONTES/Depto. Biologia Geral/Lab. Ecologia Evolutiva | : :' : Campus Universitário Prof. Darcy Ribeiro, Vila Mauricéia | `. `'` CP: 126, CEP: 39401-089, Montes Claros - MG - Brasil | `- Fone: (38) 3229-8190 | [EMAIL PROTECTED] | ICQ#: 5692561 | LinuxUser#: 205366 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] MASS need R = 2.4.0!!!
See this thread on the R-Debian list: https://stat.ethz.ch/pipermail/r-sig-debian/2006-October/000126.html -Deepayan On 10/15/06, Ronaldo ReisJunior [EMAIL PROTECTED] wrote: Hi, I have R 2.3.1 installed by debian páckage. I install only the base and recommended R packages from Debian source, all others packages I install from the R source at CRAN. But, I try to use library MASS, but I received this message: library(MASS) Error: This is R 2.3.1, package 'MASS' needs = 2.4.0 What is the problem? Inte Ronaldo -- It ain't over until it's over. -- Casey Stengel -- Prof. Ronaldo Reis Júnior | .''`. UNIMONTES/Depto. Biologia Geral/Lab. Ecologia Evolutiva | : :' : Campus Universitário Prof. Darcy Ribeiro, Vila Mauricéia | `. `'` CP: 126, CEP: 39401-089, Montes Claros - MG - Brasil | `- Fone: (38) 3229-8190 | [EMAIL PROTECTED] | ICQ#: 5692561 | LinuxUser#: 205366 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Hide line ends behind unfilled circles?
Here is a completely different solution using gplot in sna. We create an edge matrix, edges, and plot it. library(sna) edges - replace(matrix(0, 8, 8), cbind(match(x0, xx), match(x1, xx)), 1) gplot(edges, coord = cbind(xx, yy), usearrows = FALSE, vertex.col = c(black, white)[factor(aa)]) On 10/15/06, Michael Kubovy [EMAIL PROTECTED] wrote: Dear r-helpers, xx - c(0.000, 0.210, 0.714, 0.514, 1.000, 0.190, 0.590, 0.152) yy - c(0.000, 0.265, 0.256, 0.521, 0.538, 0.761, 0.821, 1.000) aa - c(19, 19, 19, 21, 19, 21, 21, 21) x0 - xx[c(1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 6, 6, 7)] y0 - yy[c(1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 6, 6, 7)] x1 - xx[c(2, 3, 3, 4, 6, 4, 5, 5, 6, 7, 7, 7, 8, 8)] y1 - yy[c(2, 3, 3, 4, 6, 4, 5, 5, 6, 7, 7, 7, 8, 8)] plot(yy ~ xx, pch = aa, cex = 3) segments(x0, y0, x1, y1) Can anyone suggest a way of insuring that the lines are hidden behind the unfilled circles? _ Professor Michael Kubovy University of Virginia Department of Psychology USPS: P.O.Box 400400Charlottesville, VA 22904-4400 Parcels:Room 102Gilmer Hall McCormick RoadCharlottesville, VA 22903 Office:B011+1-434-982-4729 Lab:B019+1-434-982-4751 Fax:+1-434-982-4766 WWW:http://www.people.virginia.edu/~mk9y/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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 stop execution
I frequently run R code using source() or just paste large chunks of code into the console window. When there is an error in the code the execution does not stop. Rather, there is a blue error message and the remainder of the code continues to execute. Questions: 1) Is there a way to change the color of the error messages from blue to another color? There are other messages (not error messages) that are also blue and to have another color would make it easier to spot errors as the screen scrolls through the code. 2) Is there a way to make execution of the sourced or pasted code to stop immediatelyt if there is an error? I did some research and found a way to make the whole session end, but that is too radical. I also tried to change the option error from NULL to stop and that also does not do the trick: the remaining pasted or sourced code is still executed. Any help will be appreciated. I am running R version 2.3.0 on Windows XP Pro. FS __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] MASS need R = 2.4.0!!!
On 15 October 2006 at 14:45, Ronaldo ReisJunior wrote: | I have R 2.3.1 installed by debian páckage. | | I install only the base and recommended R packages from Debian source, all | others packages I install from the R source at CRAN. We some 50 or more r-cran-* packages that you could apt-get install too... | But, I try to use library MASS, but I received this message: | | library(MASS) | Error: This is R 2.3.1, package 'MASS' needs = 2.4.0 | | What is the problem? A mistake I made in specifying too weak a constraint when I said 'Depends: r-base-core ( 2.3.1)' as 2.3.1-1 satisfies this. So several CRAN et al packages compiled using the different R 2.4.0 release candidates in Debian unstable slipped into testing before R 2.4.0 itself did. As discussed on r-sig-debian from where I quote from a post from two days ago: - You can either 1) go to Debian unstable and install r-base-core, r-base-dev, r-mathlib, ... from the 2.4.0 release manually (or semi-automatically using the pinning feature of apt-get), or 2) go to snapshot.debian.net and get the previous version of the failing package (here r-cran-vr), or 3) use install.packages() or update.packages() within R to overwrite to package 4) wait for R 2.4.0 to hit testing. I recommend the first choice, esp with apt-get pinning which makes it quite easy. Choice 2 is probably the least scary if 'pinning' is a bit much. Choice 3 may be the easiest. Choice 4 obviously stinks. Again, sorry about this. I'll try hard not to do that again. - The R 2.4.0 packages may actually enter testing today, so 4) may not be so onerous. They may get blocked by apparent conflict with rpy which may require a manual push by the release team. I have contacted them. Again, sorry -- this shouldn't have happened. Dirk -- Hell, there are no rules here - we're trying to accomplish something. -- Thomas A. Edison __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help 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 stop execution
On 10/15/2006 1:19 PM, Fernando Saldanha wrote:
I frequently run R code using source() or just paste large chunks of
code into the console window. When there is an error in the code the
execution does not stop. Rather, there is a blue error message and the
remainder of the code continues to execute.
Questions:
1) Is there a way to change the color of the error messages from blue
to another color? There are other messages (not error messages) that
are also blue and to have another color would make it easier to spot
errors as the screen scrolls through the code.
No, there isn't.
2) Is there a way to make execution of the sourced or pasted code to
stop immediatelyt if there is an error? I did some research and found
a way to make the whole session end, but that is too radical. I also
tried to change the option error from NULL to stop and that also does
not do the trick: the remaining pasted or sourced code is still
executed.
source() should stop when you get an error. Rather than pasting a lot
of text into the clipboard, if you want it to stop, you could use
source('clipboard')
or
source('clipboard', echo=T)
Any help will be appreciated. I am running R version 2.3.0 on Windows
XP Pro.
My advice would be to upgrade to 2.4.0 (or at least 2.3.1).
Duncan Murdoch
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] mitools, multiple imputation
On Sat, 14 Oct 2006, John Sorkin wrote: R 2.2.0 windows XP I am beginning to explore the mitools package contributed by Thomas Lumley (thank you Thomas) and I have a few questions: (1) In the examples given in the mitools documentation, the only family argument used is family=binomial. Does the package support family=gaussian and other link functions? I ran the with function with family=gaussian and I obtained results, but I am not sure if gaussian is supported by the package, so I don't know if I should trust the values I obtained. Yes, it works with family=gaussian. It isn't even restricted to glm()s -- it works with anything that has coef() and vcov() methods (2) In the documentation, the smi dataset is used i.e. data(smi). Could someone tell me how I can print the data associated with smi? I'm not sure you really want to do that: it is 5 x 12 x 1170 More useful might be str(smi) or with(smi, fun=head) However, if you do want to, either unclass(smi) or with(smi,fun=print) will work. -thomas Thomas Lumley Assoc. Professor, Biostatistics [EMAIL PROTECTED] University of Washington, Seattle __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Compact presentation of multiple figures
Dear r-helpers, I have six panels: op(par(mfrow = c(2, 3), xaxt = 'n', yaxt = 'n', pty = 's'). Each of them has a variable main = 'i', and xlab = '', ylab = '' I would like to achieve two things: (1) a common x-axis label under panel 5, and one label to the left of panels 1 and 4 (2) minimal space between the panels. I have looked for examples, and even the thorough Les paramètres graphiques didn't help me: I haven't been able to reduce the space between the rows of panels nearly enough. My best attempt: par(pty = 's', xaxt = 'n', yaxt = 'n', mfrow = c(2, 3), mgp = c(1, 0, 0), mar = c(1.1, 2.1, 1.1, 0.1)) _ Professor Michael Kubovy University of Virginia Department of Psychology USPS: P.O.Box 400400Charlottesville, VA 22904-4400 Parcels:Room 102Gilmer Hall McCormick RoadCharlottesville, VA 22903 Office:B011+1-434-982-4729 Lab:B019+1-434-982-4751 Fax:+1-434-982-4766 WWW:http://www.people.virginia.edu/~mk9y/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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 can i compute the average of three blocks for each column ?
Dear all,
I want to compute the average of the three blocks for each x-variable which
is equal slide in the code below. How can I do that ?
block x1x2x3x4x5
12322 232423
12125 262139
123242223 23
220 21 232428
2 32 2334 24 26
2 19 34341334
3 12 32 ´ 233419
3 23 24 252627
3 12 78 232424
# read table of data for this slide=(x1)
a-read.table(file = slide[i],header=T,sep='\t',na.strings=NA)
#length(a$ID)
#Eleminate Neg. and Pos. controls from the dataset. The logical negation of
the %in% function,
#tells subset to only select those row where the ID column does not contain
either empty or none
new - subset(a,!ID %in% c(empty,none, ))
#length(new$ID)
#new[1:20,c(1,4,5,9)]
#five first columns give position identifiers, include a column with block
layout=new[,1:5]
layout[1:30,]
#9th columns which give the median foreground =values of x-variables
fg1=as.matrix(new[,9])
length(fg1)
mean(fg1) # calculate the mean of x1
I try to do something like :##
block1=fg1[layout$Block==1,]
block2=fg1[layout$Block==1,]
block2=fg1[layout$Block==1,]
average=(block1+block2+block3)/3
but it did not work.
## How can i calculate the means of remaining x_variables?
# Read data for the remaining slides =x2,x3,x4,x5 ###
for (i in 2:num.slides){
na1 - strsplit(na[[i]][k],.txt)
na2 - strsplit(na1[[1]][1],-)
bat=na2[[1]][1]
sli=na2[[1]][2]
nslide - cbind(nslide,as.numeric(sli))
# nslide is a vector giving the number of the slide in the batch
# read table of data for this slide
a-read.table(file=slide[i],header=T,sep='\t',na.strings=NA)
new- subset(a,!ID %in% c(empty,none, ))
# append FG data to the matrices containing the slides already read
fg1=cbind(fg1,as.matrix(new[,9]))
}
colnames(fg1)=nslide
fg-data.frame(peptide=c(new$Name),fg1)
fg - edit(fg)
# Another question : I have three graphs which are displayed one after
one with a large space between them. Can I move these graph closer each other
by making them bigger and how ? Below is the code that i have written for
plotting the graphs.
par(mfrow=c(3,1))
for (j in 1:3)
{
boxplot(split(pos$y[pos$Block==j],pos$Slide[pos$Block==j]), col=lightgray,
cex=.65, outline=TRUE, main=paste(Positive Controls Block,j))
}
Thank you for your help,
Regards,
Yen
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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
Hello again,
Thanks for answering my questions.
If my program is now:
pdf('Test.pdf')
BL-c(1.97,8.04,2.54,10.53,4.85,1.73)
LR-c(0.85,0.86,8.33,04.18,6.26,2.40)
Q-c(0.00,0.00,1.92,01.92,4.48,0.00)
cols - ifelse(Q!=0, orange, green)
Graph- barplot(LR, main='LR Value',col=cols, border='black', space=0.05,
width=(BL), xlab='Length', ylab='LR block')
axis(1, at=Graph, sprintf('%0.2f',BL))
mtext(1, at=Graph, text=ifelse(Q!=0, sprintf('%0.2f',Q), ), line=2)
mtext(1, at=par('usr')[1], text='BL', line=1)
mtext(1, at=par('usr')[1], text='Var', line=2)
abline(h=3.8, col='blue')I want to add alpha=5% at the end
of this line
And now I want to write alpha=5% at the end of the line that I have on my
graph.
Thanks,
Mohsen
On 10/15/06, Marc Schwartz [EMAIL PROTECTED] wrote:
On Sat, 2006-10-14 at 23:52 -0400, Mohsen Jafarikia wrote:
Hello everyone,
I have the following program to draw a barplot.
MP-read.table(file='AR.out')
names(MP)-c('BN','BL','LR','Q')
Graph- barplot(MP$LR, main='LR Value', col='orange', border='black',
space=
0.05, width=(MP$BL), xlab='Length', ylab='LR each')
axis(1, at=Graph, sprintf('%0.2f',MP$BL))
mtext(1, at=Graph, text=sprintf('%0.2f',MP$Q), line=2)
mtext(1, at=par('usr')[1], text='BL', line=1)
mtext(1, at=par('usr')[1], text='Var', line=2)
abline(h= 3.841,col='blue')
abline(h=6.635,col='red')
I have two more questions about the graph that I have:
1) I want to write the 'Q' only when it is not equal to zero.
2) I would like to change the bars when 'Q' is not zero. For example,
from
orange to green.
I would appreciate your input to this question.
Thanks,
Mohsen
It would be helpful to have the data that you are working with so that
we can provide each other a working example.
However, here are some hints:
1.
mtext(1, at = Graph,
text = ifelse(MP$Q != 0, sprintf('%0.2f',MP$Q), ),
line=2)
2.
cols - ifelse(MP$Q != 0, orange, green)
barplot(..., col = cols, ...)
See ?ifelse
HTH,
Marc Schwartz
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Compact presentation of multiple figures
This will get you started x - rnorm(20) y - rnorm(20) par(pty = 's', xaxt = 'n', yaxt = 'n', mfrow = c(2, 3), mgp= c(1, 0, 0), mar = c(1.1, 2.1, 1.1, 0.1)) plot(y ~ x, main=1, xlab=, ylab=y1) plot(y ~ x, main=2, xlab=, ylab=) plot(y ~ x, main=3, xlab=, ylab=) plot(y ~ x, main=4, xlab=, ylab=y4) plot(y ~ x, main=5, xlab=, ylab=, xaxt=s) plot(y ~ x, main=6, xlab=, ylab=) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help 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
Try adding
text(31,3.8,expression(paste(alpha==5,%)),pos=3)
On 15/10/06, Mohsen Jafarikia [EMAIL PROTECTED] wrote:
Hello again,
Thanks for answering my questions.
If my program is now:
pdf('Test.pdf')
BL-c(1.97,8.04,2.54,10.53,4.85,1.73)
LR-c(0.85,0.86,8.33,04.18,6.26,2.40)
Q-c(0.00,0.00,1.92,01.92,4.48,0.00)
cols - ifelse(Q!=0, orange, green)
Graph- barplot(LR, main='LR Value',col=cols, border='black', space=0.05,
width=(BL), xlab='Length', ylab='LR block')
axis(1, at=Graph, sprintf('%0.2f',BL))
mtext(1, at=Graph, text=ifelse(Q!=0, sprintf('%0.2f',Q), ), line=2)
mtext(1, at=par('usr')[1], text='BL', line=1)
mtext(1, at=par('usr')[1], text='Var', line=2)
abline(h=3.8, col='blue')I want to add alpha=5% at the end
of this line
And now I want to write alpha=5% at the end of the line that I have on my
graph.
Thanks,
Mohsen
On 10/15/06, Marc Schwartz [EMAIL PROTECTED] wrote:
On Sat, 2006-10-14 at 23:52 -0400, Mohsen Jafarikia wrote:
Hello everyone,
I have the following program to draw a barplot.
MP-read.table(file='AR.out')
names(MP)-c('BN','BL','LR','Q')
Graph- barplot(MP$LR, main='LR Value', col='orange', border='black',
space=
0.05, width=(MP$BL), xlab='Length', ylab='LR each')
axis(1, at=Graph, sprintf('%0.2f',MP$BL))
mtext(1, at=Graph, text=sprintf('%0.2f',MP$Q), line=2)
mtext(1, at=par('usr')[1], text='BL', line=1)
mtext(1, at=par('usr')[1], text='Var', line=2)
abline(h= 3.841,col='blue')
abline(h=6.635,col='red')
I have two more questions about the graph that I have:
1) I want to write the 'Q' only when it is not equal to zero.
2) I would like to change the bars when 'Q' is not zero. For example,
from
orange to green.
I would appreciate your input to this question.
Thanks,
Mohsen
It would be helpful to have the data that you are working with so that
we can provide each other a working example.
However, here are some hints:
1.
mtext(1, at = Graph,
text = ifelse(MP$Q != 0, sprintf('%0.2f',MP$Q), ),
line=2)
2.
cols - ifelse(MP$Q != 0, orange, green)
barplot(..., col = cols, ...)
See ?ifelse
HTH,
Marc Schwartz
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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 Barron
Said Business School
University of Oxford
Park End Street
Oxford OX1 1HP
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] Execution halting of lmer on UNIX when no problem on windows
Dear R-users, I have a frustrating problem that I am hoping has a simple fix. I am running a series of lmer models from the lme4 package of the general form: model-lmer(y~x1 + x2 . + xn + (1|site),data=dataframe,family=poisson,method=Laplace,control=list(usePQL=FALSE,msVerbose=TRUE)) where the same model is executed multiple times on a bootstrapped dataframe. For each bootstrapped model run the resulting model object is used to return the AIC (and models are then compared using a bootstrapped weight - frequency of runs a given model had the lowest AIC). This works just fine when running interactively on my windows machine (so there is nothing the matter with the code, hence I have not bored you with the details here) however when I submit the job as a batch to a UNIX system it usually (but not always) fails and the execution is halted after an error message is produced: Error in objective (.par,...): Leading minor of order 1 in downdated X'X is not positive definite In addition: There were 12 warnings (use warnings() to see them) Error in logLik(model) : no applicable method for logLik Execution halted On windows when the execution of a single model run fails to estimate the logLik (an unstable model...) R just continues past the error and still runs to the end of the script (i.e. runs through all the bootstraps), and I can then inspect the output and any errors at the end. My question, then is when using UNIX on batch mode, how can I get the job to NOT halt it's execution on the production of an unstable model (not positive definite) and continue running? If the models are not bootstrapped then the script is executed without any problem in UNIX (so the script is successfuly submitted as a batch job), so it seems that some of the bootstrap runs are producing unstable models in rare instances, but just one is sufficient to halt the execution of the script. I am running R.2.3.1 on both Windows and UNIX, Many thanks in advance, Toby Gardner School of Environmental Sciences University of East Anglia Norwich, NR4 7TJ United Kingdom Email: [EMAIL PROTECTED] Website: www.uea.ac.uk/~e387495 [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Need help with barplots
On 2006-10-13 10:55, David Barron wrote: First, produce two barplots for comparison: par(mfrow=c(2,1) ) barplot(VADeaths,beside=TRUE) barplot(VADeaths) The same information is in both plots; in the top, it is displayed as 5 separate bars for each group, and in the stacked plot it is shown as 5 separate regions in each of the four bars. The hight of each of these regions is the same as the hight of the corresponding bar in the side-by-side plot. The stacked plot enables you to see overall differences more easily (easier to see that the death rate is highest for Urban Males), but it is harder to compare the sizes of the categories. Take a look at ?spineplot. HTH, Henric On 13/10/06, laba diena [EMAIL PROTECTED] wrote: I`ve read all the manuals and still couln`t find what is the difference between the stacked and side-by-side barplots ? Could you explain me ? [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Compact presentation of multiple figures
Thanks Richard, That is roughly what I did. However, (1) the distance between the rows is still quite large compared to the distance between the columns; (2) (less important) I still have two ylab, one for each row. On Oct 15, 2006, at 5:12 PM, Richard M. Heiberger wrote: This will get you started x - rnorm(20) y - rnorm(20) par(pty = 's', xaxt = 'n', yaxt = 'n', mfrow = c(2, 3), mgp= c(1, 0, 0), mar = c(1.1, 2.1, 1.1, 0.1)) plot(y ~ x, main=1, xlab=, ylab=y1) plot(y ~ x, main=2, xlab=, ylab=) plot(y ~ x, main=3, xlab=, ylab=) plot(y ~ x, main=4, xlab=, ylab=y4) plot(y ~ x, main=5, xlab=, ylab=, xaxt=s) plot(y ~ x, main=6, xlab=, ylab=) _ Professor Michael Kubovy University of Virginia Department of Psychology USPS: P.O.Box 400400Charlottesville, VA 22904-4400 Parcels:Room 102Gilmer Hall McCormick RoadCharlottesville, VA 22903 Office:B011+1-434-982-4729 Lab:B019+1-434-982-4751 Fax:+1-434-982-4766 WWW:http://www.people.virginia.edu/~mk9y/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Mute script
Hi, I tried to run the following script with R 2.4.0. The data stuff is commented out because data are already in memory. #db - read.csv(normais.csv, sep=;, quote=, header=T) sink(normais-chi.txt, append=T, type = output, split=T) #sink(normais-chi.txt, append=T) table(db$agua, db$mBerg) chisq.test(db$agua, db$mBerg, correct=F) print(table(db$agua, db$mBerg)) print(chisq.test(db$agua, db$mBerg, correct=F)) # here there are several other table/chisq.test pairs sink() #rm(db) By introducing syntax errors in the script error messages will appear but normal output won't neither at screen nor at the sink file. On the other hand, by inserting print commands I get both screen and file outputs. Is it necessary to use the print command in every line? In interactive mode the simple command will produce the proper output. Or do I have a configuration problem? Haven't found the answer elsewhere in teh Net. Thanks in advance. -- Alexandre Aguiar Independent consultant for medical research SPS Consultoria Voice: +55-11-9320-2046 Fax: +55-11-5549-8760 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] overlapping intervals
Hello everybody, I have two series of intervals, and I'd like to output the shared regions. For example: series1-cbind(Start=c(10,21,40,300),End=c(20,26,70,350)) series2-cbind(Start=c(25,60,210,500),End=c(40,100,400,1000)) series1 Start End [1,]10 20 [2,]21 26 [3,]40 70 [4,] 300 350 series2 Start End [1,]25 40 [2,]60 100 [3,] 210 400 [4,] 500 1000 I'd like to have something like this as result: shared Start End [1,]25 26 [2,]60 70 [3,] 300 350 I found this post, but the solution finds the regions shared across all the intervals. http://finzi.psych.upenn.edu/R/Rhelp02a/archive/59594.html Can anybody help me with this? Thanks Giovanni __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Compact presentation of multiple figures
I would do something like this in lattice. It defaults to something prettier and gives more control if you don't want the defaults. tmp - data.frame(x=rnorm(120), y=rnorm(120), a=factor(rep(1:6, 20))) xyplot(y ~ x | a, data=tmp) xyplot(y ~ x | a, data=tmp, aspect=1, layout=c(3,2), scales=list(x=list(alternating=c(0,1,0)), y=list(alternating=FALSE))) Back to regular graphics. I misunderstood your request for a single y label. This comes closer. I used oma to squeeeze the two rows closer together. I used two axis(1) statements. The first puts the ticks on the border of the box. The second one puts the labels far enough away that they aren't cluttered. Similar adjustment is called for on the main title of each panel. x - rnorm(20) y - rnorm(20) old.par - par(pty = 's', xaxt = 'n', yaxt = 'n', mfrow = c(2, 3), mar = c(1.1, 2.1, 1.1, 0.1), oma=c(12,2,0,2)) plot(y ~ x, main=1, xlab=, ylab=) plot(y ~ x, main=2, xlab=, ylab=) plot(y ~ x, main=3, xlab=, ylab=) plot(y ~ x, main=4, xlab=, ylab=) mtext(Y label, side=2, at=2, line=2) plot(y ~ x, main=5, xlab=, ylab=) axis(1, xaxt=s, labels=FALSE) axis(1, xaxt=s, tick=FALSE, line=1) plot(y ~ x, main=6, xlab=, ylab=) par(old.par) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] extracting nested variances from lme4 model
On 10/15/06, Douglas Bates [EMAIL PROTECTED] wrote: On 10/14/06, Spencer Graves [EMAIL PROTECTED] wrote: You want to estimate a different 'cs' variance for each level of 'trth', as indicated by the following summary from your 'fake data set': tapply(dtf$x, dtf$trth, sd) FALSE TRUE 1.532251 8.378206 Since var(x[trth]) var(x[!trth]), I thought that the following should produce this: mod1.-lmer( x ~ (1|rtr)+ (trth|cs) , data=dtf) Unfortunately, this generates an error for me: Error in lmer(x ~ (1 | rtr) + (trth | cs), data = dtf) : the leading minor of order 1 is not positive definite I do not understand this error, and I don't have other ideas about how to get around it using 'lmer'. Hmm. It's an interesting problem. I can tell you where the error message originates but I can't yet tell you why. The error message originates when Cholesky decompositions of the variance-covariance matrices for the random effects are generated at the initial estimates. It looks as if one of the model matrices is not being formed correctly for some reason. I will need to do more debugging. OK, I have worked out why the error occurs and will upload lme4_0.9975-4, which fixes this problem and also has some speedups when fitting generalized linear mixed models. Just for the record, I had assumed that the cholmod_aat function in the CHOLMOD sparse matrix library (the function cholmod_aat calculates AA' from a sparse matrix A) returned a matrix in which the columns are sorted in increasing order of the rows. That is not always correct but the situation is easily remedied. (The typical way to sort the columns in sparse matrices is to take the transpose of the transpose, which had me confused when I first saw it in some of Tim Davis's code. It turns out that a side effect of this operation is to sort the columns in increasing order of the rows.) With this patched lme4 package the model Spencer proposed can be fit but the variance-covariance matrix of the two-dimensional random effects for the cs factor is singular because of the design (the level of trth is constant within each level of cs). I'm not sure how to specify the desired model in lmer. (fm2 - lmer(x ~(1|rtr) + (trth|cs), dtf)) Linear mixed-effects model fit by REML Formula: x ~ (1 | rtr) + (trth | cs) Data: dtf AIC BIC logLik MLdeviance REMLdeviance 252.5 260.9 -121.2 244.8242.5 Random effects: Groups NameVariance Std.Dev. Corr cs (Intercept) 1.9127 1.3830 trthTRUE24.1627 4.9156 1.000 rtr (Intercept) 1.9128 1.3830 Residual 18.5278 4.3044 number of obs: 40, groups: cs, 10; rtr, 4 Fixed effects: Estimate Std. Error t value (Intercept) -0.2128 1.2723 -0.1673 Warning message: nlminb returned message false convergence (8) in: `LMEoptimize-`(`*tmp*`, value = list(maxIter = 200, tolerance = 1.49011611938477e-08, It seems to me that 'lme' in library(nlme) should also be able to handle this problem, but 'lme' is designed for nested factors, and your 'rtr' and 'cs' are crossed. If it were my problem, I think I'd study ch. 4 and possibly ch. 5 of Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer) on 'pdMat' and 'corrStruct' classes, because I believe those may support models with crossed random effects like in your example AND might also support estimating different variance components for different levels of a fixed effect like 'trth' in your example. If we are lucky, someone else might educate us both. I'm sorry I couldn't answer your question, but I hope these comments might help in some way. Spencer Graves Frank Samuelson wrote: I have a model: mod1-lmer( x ~ (1|rtr)+ trth/(1|cs) , data=dtf) # Here, cs and rtr are crossed random effects. cs 1-5 are of type TRUE, cs 6-10 are of type FALSE, so cs is nested in trth, which is fixed. So for cs I should get a fit for 1-5 and 6-10. This appears to be the case from the random effects: mean( ranef(mod1)$cs[[1]][1:5] ) [1] -2.498002e-16 var( ranef(mod1)$cs[[1]][1:5] ) [1] 23.53083 mean( ranef(mod1)$cs[[1]][6:10] ) [1] 2.706169e-16 var( ranef(mod1)$cs[[1]][6:10] ) [1] 1.001065 However VarCorr(mod1) gives me only a single variance estimate for the cases. How do I get the other variance estimate? Thanks for any help. -Frank My fake data set: --- data: $frame x rtr trth cs 1 -4.4964808 1 TRUE 1 2 -0.6297254 1 TRUE 2 31.8062857 1 TRUE 3 42.7273275 1 TRUE 4 5 -0.6297254 1 TRUE 5 6 -1.3331767 1 FALSE 6 7 -0.3055796 1 FALSE 7 81.3331767 1 FALSE 8 90.1574314 1 FALSE 9 10 -0.1574314 1 FALSE 10 11 -3.0958803 2
Re: [R] overlapping intervals
Not the most efficient and requires integer values (maybe less than 1M). My results show an additional overlap at 40 - start end were the same -- does this count? If not, just delete rows that are the same in both columns. series1-cbind(Start=c(10,21,40,300),End=c(20,26,70,350)) series2-cbind(Start=c(25,60,210,500),End=c(40,100,400,1000)) x1 - x2 - logical(max(series1, series2)) # vector FALSE x1[unlist(mapply(seq, series1[,1], series1[,2]))] - TRUE x2[unlist(mapply(seq, series2[,1], series2[,2]))] - TRUE r - rle(x1 x2) # determine overlaps offset - cumsum(r$lengths) (z - cbind(offset[r$values] - r$lengths[r$values] + 1, offset[r$values])) [,1] [,2] [1,] 25 26 [2,] 40 40 [3,] 60 70 [4,] 300 350 # if you don't like dups for overlaps (@40) z[z[,1] != z[,2],] [,1] [,2] [1,] 25 26 [2,] 60 70 [3,] 300 350 On 10/15/06, Giovanni Coppola [EMAIL PROTECTED] wrote: Hello everybody, I have two series of intervals, and I'd like to output the shared regions. For example: series1-cbind(Start=c(10,21,40,300),End=c(20,26,70,350)) series2-cbind(Start=c(25,60,210,500),End=c(40,100,400,1000)) series1 Start End [1,]10 20 [2,]21 26 [3,]40 70 [4,] 300 350 series2 Start End [1,]25 40 [2,]60 100 [3,] 210 400 [4,] 500 1000 I'd like to have something like this as result: shared Start End [1,]25 26 [2,]60 70 [3,] 300 350 I found this post, but the solution finds the regions shared across all the intervals. http://finzi.psych.upenn.edu/R/Rhelp02a/archive/59594.html Can anybody help me with this? Thanks Giovanni __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Jim Holtman Cincinnati, OH +1 513 646 9390 What is the problem you are trying to solve? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] MASS4 page 360 cumulative hazard plot
I cannot reproduce Fig. 13.2 in MASS4. plot(gehan.surv,fun='cloglog') Warning message: 2 x values = 0 omitted from logarithmic plot in: xy.coords(x, y, xlabel, ylabel, log) and the x-axis is badly scaled. I was wondering if someone can help Regards Ross Darnell __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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] Contour plot of bivariate log-likelihood function
Thanks to prior help from the list, I've made progress in adapting a
univariate log-likelihood function to a bivariate one for Box-Cox
transformations. However, I'm having trouble plotting values of the
log-likelihood function (z) for given range of lambda(1) and lambda(2).
The current result of my function (below) is z. How can I adapt it to
output the values lambda(1), lambda(2) (for each row of all_lambda) and
z and then plot using contour()?
I'm also stuck on how to allow the variable xL to == log(x) when
lambda(1)==0 and/or lambda(2)==0.
Thanks.
--Dale
x-read.table(http://www.stat.lsu.edu/faculty/moser/exst7037/jwdata/T4-1.txt,col.names=x1;)
x-read.table(http://www.stat.lsu.edu/faculty/moser/exst7037/jwdata/T4-5.txt,col.names=x2;)
X - as.matrix(cbind(x1,x2))
bvboxcox - function(X, min, max, step)
{
seq - seq(min,max,step)
all_lambda - expand.grid(seq,seq)
all_lambda - as.matrix(all_lambda)
c - nrow(all_lambda)
res - numeric(c)
for (i in seq(1:c)){
lambda - all_lambda[i,]
xL - (X^lambda - 1)/lambda
n - nrow(X)
t1 - (-n/2)*log(det(cov(xL)))
gsum - apply(X,2,function(x) sum(log(x)))
t2 - sum((lambda - 1) * gsum)
res[i] - t1+t2
}
res
}
## For example ...
bvboxcox(X, 0.10, 0.17, 0.01)
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] Input buffer overflow
In gsubfn I replace matches with strings that represent calls to a function and then perform paste(eval(parse(text= ...)), collapse = ) on the result. One user of gsubfn is using it with very long strings (over 20,000 characters) and the parse is giving an input buffer overflow. Here is an artificial example: s - paste(rep(X, 25000), collapse = ) out - parse(text = shQuote(s)) Error in parse(text = shQuote(s)) : input buffer overflow Is there a way to increase the limit? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/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-pkgs] New package Ryacas
Ryacas is an R interface to the free yacas computer algebra system. Ryacas allows one to send R expressions, unprocessed yacas strings and certain other R objects to a separate yacas process from R and get back the result. It also has facilities for manipulating yacas strings and R expressions destined for yacas processing. It can be used for exact arithmetic, symbolic math, ASCII pretty printing and translating R to TeX. Online info. For overview, pointers to additional information, installation instructions and a sample session see: http://code.google.com/p/ryacas/ The vignettes can be viewed online here: http://ryacas.googlecode.com/svn/trunk/inst/doc/Ryacas.pdf http://ryacas.googlecode.com/svn/trunk/inst/doc/Ryacas-Sym.pdf More. Once Ryacas is installed, pointers to additional information can be found with these R commands: library(Ryacas) package?Ryacas --- Rob Goedman, goedman at mac dot com Gabor Grothendieck, ggrothendieck at gmail dot com Søren Højsgaard, Soren.Hojsgaard at agrsci dot dk Ayal Pinkus, apinkus at xs4all dot nl ___ R-packages mailing list R-packages@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help 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 Correcting Codes, Simplex
On 10/8/06, Egert, Bjoern [EMAIL PROTECTED] wrote: Hello, Is there a way in R to construct an (error correcting) binary code e.g. for an source alphabet containing integers from 1 to say 255 with the property that each pair of distinct codewords of length m is at Hamming distance exactly m/2 ? I was suggested to use so called simplex codes, which should be fairly standard, but I haven't found a direct way via R packages to do so, that's why I ask whether there might be in indirect way to solve this problem. Example: v1 =c(1,2,3,4) v2 =c(1,2,5,6) similarity(v1,v2)=0.5, (because 2 out of 4 elements are equal). Obviously, a binary representation of would yield a different similarity of: binary(v1) =001 010 011 100 binary(v1) =001 010 101 110 similarity(binary(v1),binary(v2))= 9/12 Remark: The focus here is not on error correction, but rather the binary encoding retaining similarity of the elements of vectors. Many thanks, Bjoern Bjoern, NB: I'm an R newbie and I only know a bit about error correcting codes. I haven't seen any responses to your questions and I don't know if you still have a need, but it is certainly possible to construct forward error correction codes with all the great math capability in R. It seems you want to generate code words that still have the original bits present. These are systematic codes and there are lots of them available to use. Many codes are specified by the code word length (n), number of original data bits in each code word (k), and the minimum Hamming distance of the code words (d) as a [n,k,d] code. Simplex Codes have these parameters: [2^k - 1, k, 2^(k - 1)]. These codes could be generated as a simple matrix multiply in R, but are you sure that's what you want? The code words will be quite long. Regards, Richard Graham __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
