[R] clusterCrit package produces Nan
I'm trying to compute the Silhouette value of a clustering partition done with SOM method using the function intCriteria of the clusterCrit package but it returns me a Nan value: somebody knows why? To simplify my case, here there is an example: test [,1] [,2] [,3] [,4] [,5] [,6] [1,] 37.28577 8.902218 17.93830 38.80381 6.825286 18.44061 [2,] 37.25598 8.693962 18.09256 38.81784 7.005092 18.17899 [3,] 37.54612 8.262074 18.07639 38.87837 6.592799 18.31604 [4,] 37.56661 8.651182 17.98653 38.76980 6.596529 18.46937 [5,] 37.59454 8.546921 17.93558 39.00177 6.508707 18.40102 [6,] 37.62195 8.422909 18.03132 38.86634 7.024104 18.30201 [7,] 37.19365 8.662608 18.01295 38.37456 7.173594 18.24273 [8,] 37.63161 8.495688 18.05909 38.92736 6.884675 18.44871 [9,] 37.30077 8.488438 18.02636 38.60844 7.004214 18.45600 [10,] 37.52518 8.610817 18.00498 38.57547 6.877532 18.36552 datissimi = som(test, grid=somgrid(xdim=2, ydim=4, topo='rectangular'), rlen=500, keep.data=T) somissimo = datissimi$unit.classif intCriteria(test, somissimo, Silhouette) $silhouette [1] NaN Thanks, Paola __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] reference for variance in svykm
Dear all, I am using the svykm function (in the survey package) to estimate survival in a two-phase study. I would like to know if there is any published reference for the estimator of the variance of survival curve under general sampling in order to cite it? The code uses the linearization strategy by Williams (Product-Limit Survival Functions with Correlated Survival Times,Lifetime Data Analysis 1995), but I found the formulas for general complex sampling in a technical note written by Thomas Lumley (http://staff.washington.edu/tlumley/survey/survcurve.pdf) and I was wondering if it was ever fomally validated and published. Thanks, Paola [[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] error in comparisons using glht function
hi everybody this is my last chance, I dont know what else to do I was trying to find my multiple comparisons using the code in R glht next is my function and the error I get vol.interallaovsol-aov(scoresol~ffjob+ffage+ffstudies+ffjob:ffstudies+ffjob:ffage+ffstudies:ffage+ffjob:ffage:ffstudies,data=vol18.df) vol.mcsol-glht(vol.interallaovsol, linfct=mcp(ffage=Tukey)) Warning messages: 1: In mcp2matrix(model, linfct = linfct) : covariate interactions found -- default contrast might be inappropriate 2: In glht.matrix(model = list(coefficients = c(4.83, -0.8846668, : 157 out of 210 coefficients not estimable in ‘model’ I dont have a clue what is this. Any idea??thanks -- View this message in context: http://r.789695.n4.nabble.com/error-in-comparisons-using-glht-function-tp4648467.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Contrasts in manova
Hi everybody I am trying to find contrast in MANOVA. I used next code contrasts(ffage)-ctr contrasts(ffage) MANOVA.agec-manova(Y1~ffage,data=vol18.df) summary(MANOVA.agec, split =list (ffage=list(0-17 v over 18=0, 18-25 v over 26=1, 26-31 v over 32=2, 32-42 v over 43=3, 43-65 v 66+=4))) But the output was only the overall fit contrasts(ffage)-ctr contrasts(ffage) [,1] [,2] [,3] [,4] [,5] 050000 1 -14000 2 -1 -1300 3 -1 -1 -120 4 -1 -1 -1 -11 5 -1 -1 -1 -1 -1 MANOVA.agec-manova(Y1~ffage,data=vol18.df) summary(MANOVA.agec, split =list (ffage=list(0-17 v over 18=0, 18-25 v over 26=1, 26-31 v over 32=2, 32-42 v over 43=3, 43-65 v 66+=4))) Df Pillai approx F num Df den Df Pr(F) ffage 5 0.30607 2.1543 20520 0.002681 ** Residuals 130 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 I try to use linearHypothesis, but i couldnt find anything, maybe I didnt use in a corect way, *any advice is welcome* I always get the error Error in solve.default(wcrossprod(model.matrix(model), w = wts)) : Lapack routine dgesv: system is exactly singular Thanks in advance Paola -- View this message in context: http://r.789695.n4.nabble.com/Contrasts-in-manova-tp4648306.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] error in lm
Hi everybody I am trying to run the next code but I have the next problem Y1-cbind(score.sol, score.com.ext, score.pur) vol.lm-lm(Y1~1, data=vol14.df) library(MASS) stepAIC(vol.lm,~fsex+fjob+fage+fstudies,data=vol14.df) Start: AIC=504.83 Y1 ~ 1 Error in addterm.mlm(fit, scope$add, scale = scale, trace = max(0, trace - : no addterm method implemented for mlm models Does anybody knows why I get this error? Thanks in advance Paola -- View this message in context: http://r.789695.n4.nabble.com/error-in-lm-tp4647840.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Change color in forest.rma (metafor)
Thank you so much!!! Could you tell me also how to change the size of the chart? There is not enough space below the chart to add the arrows! 2011/8/28 Uwe Ligges-3 [via R] ml-node+3774557-1567708350-262...@n4.nabble.com On 26.08.2011 15:50, Paola Tellaroli wrote: I lied, that was not my last question: how can I add two arrows at the bottom with the words in favor of A / B? This is not specified in the pdf and with text I have the impression that I can't add text below the x-axis. You can, see ?par and its xpd argument. Uwe Ligges 2011/8/26 Paola Tellaroli[hidden email]http://user/SendEmail.jtp?type=nodenode=3774557i=0 Dear Prof. Viechtbauer, thank you so much for your help and kindness. Clearly graphs are the minor problem in our work, and the parameters and options that can vary in R are so many that it is obvious that you can not expect to change everything you want! Your suggestions are very helpuf, but I have one last question. I'm trying to copy the style of a forest plot that I've seen and I like (the one in the attached file, page 1034): can I do this in R? Best wishes, *Paola* 2011/8/25 Viechtbauer Wolfgang (STAT)-2 [via R] [hidden email] http://user/SendEmail.jtp?type=nodenode=3774557i=1 The color of the squares is also currently hard coded. The thing is, there are so many different elements to a forest plot (squares, lines, polygons, text, axes, axis labels, etc.), if I would add arguments to set the color of each element, things would really get out of hand (as far as I am concerned, there are already too many arguments to begin with). I can think of one possibility: I could allow the col argument to accept a vector of colors and then apply the different elements of that vector to the various elements in the plot. Of course, there is also a limit to how far that can be taken. For example, what if somebody wants to have a different color for *one* of the squares and a different color for the other squares? Another possibility is to do some post-processing with other software. One can create the forest plot in R, save it for example as a postscript file, and the edit the plot in other software. Yes, I prefer it if I can create the plot in R and have it exactly the way I want it (without having to do any post-processing), but sometimes that may not be possible. Note that you can always add whatever you want to a plot created by the forest() function after it has been drawn. You can add text, lines, squares, polygons, whatever in any color you desire (e.g., with the text(), segments(), points(), polygon() functions). So, you could also just plot over the squares with: points(yi, 4:1, pch=15, col=red) To get rid of the black squares that are drawn by the forest function, add psize=0 as an argument in forest() (this will make the size of squares equal to 0, so essentially, they are invisible). If you want to make the size of the points inversely proportional to some function of the precision of the estimates, use points() together with the cex argument. For example: wi- 1/sqrt(vi) psize- wi/sum(wi) psize- (psize - min(psize)) / (max(psize) - min(psize)) psize- (psize * 1.0) + 0.5 points(yi, 4:1, pch=15, col=red, cex=psize) Best, Wolfgang -Original Message- From: Paola Tellaroli [mailto:[hidden email] http://user/SendEmail.jtp?type=nodenode=3768683i=0] Sent: Thursday, August 25, 2011 10:57 To: Viechtbauer Wolfgang (STAT) Cc: [hidden email] http://user/SendEmail.jtp?type=nodenode=3768683i=1; Bernd Weiss Subject: Re: [R] Change color in forest.rma (metafor) Thank you for your attention and help! In this way I get the diamond coloured, but actually I would have the squares representing the values of the individual studies coloured. Is it somehow possible? Paola 2011/8/24 Viechtbauer Wolfgang (STAT) [hidden email]http://user/SendEmail.jtp?type=nodenode=3768683i=2 Thank you, Bernd, for looking into this. Yes, at the moment, the color of the summary estimate for models without moderators is hard-coded (as black). I didn't think people may want to change that. I guess I was wrong =) A dirty solution for the moment is to add: addpoly(dfs, efac=6, row=-1, col=red, border=red, annotate=F, mlab=) after the call to forest(). You will get a warning message (since the border argument gets passed to the text() function inside addpoly() and that's not a par for text), but you can just ignore that. Best, -- Wolfgang Viechtbauer Department of Psychiatry and Neuropsychology School for Mental Health and Neuroscience Maastricht University, P.O. Box 616 6200 MD Maastricht, The Netherlands Tel: +31 (43) 368-5248 Fax: +31 (43) 368-8689 Web: http://www.wvbauer.com -Original Message- From: Bernd Weiss [mailto
Re: [R] Change color in forest.rma (metafor)
Dear Prof. Viechtbauer, thank you so much for your help and kindness. Clearly graphs are the minor problem in our work, and the parameters and options that can vary in R are so many that it is obvious that you can not expect to change everything you want! Your suggestions are very helpuf, but I have one last question. I'm trying to copy the style of a forest plot that I've seen and I like (the one in the attached file, page 1034): can I do this in R? Best wishes, *Paola* 2011/8/25 Viechtbauer Wolfgang (STAT)-2 [via R] ml-node+3768683-1225159815-262...@n4.nabble.com The color of the squares is also currently hard coded. The thing is, there are so many different elements to a forest plot (squares, lines, polygons, text, axes, axis labels, etc.), if I would add arguments to set the color of each element, things would really get out of hand (as far as I am concerned, there are already too many arguments to begin with). I can think of one possibility: I could allow the col argument to accept a vector of colors and then apply the different elements of that vector to the various elements in the plot. Of course, there is also a limit to how far that can be taken. For example, what if somebody wants to have a different color for *one* of the squares and a different color for the other squares? Another possibility is to do some post-processing with other software. One can create the forest plot in R, save it for example as a postscript file, and the edit the plot in other software. Yes, I prefer it if I can create the plot in R and have it exactly the way I want it (without having to do any post-processing), but sometimes that may not be possible. Note that you can always add whatever you want to a plot created by the forest() function after it has been drawn. You can add text, lines, squares, polygons, whatever in any color you desire (e.g., with the text(), segments(), points(), polygon() functions). So, you could also just plot over the squares with: points(yi, 4:1, pch=15, col=red) To get rid of the black squares that are drawn by the forest function, add psize=0 as an argument in forest() (this will make the size of squares equal to 0, so essentially, they are invisible). If you want to make the size of the points inversely proportional to some function of the precision of the estimates, use points() together with the cex argument. For example: wi - 1/sqrt(vi) psize - wi/sum(wi) psize - (psize - min(psize)) / (max(psize) - min(psize)) psize - (psize * 1.0) + 0.5 points(yi, 4:1, pch=15, col=red, cex=psize) Best, Wolfgang -Original Message- From: Paola Tellaroli [mailto:[hidden email]http://user/SendEmail.jtp?type=nodenode=3768683i=0] Sent: Thursday, August 25, 2011 10:57 To: Viechtbauer Wolfgang (STAT) Cc: [hidden email]http://user/SendEmail.jtp?type=nodenode=3768683i=1; Bernd Weiss Subject: Re: [R] Change color in forest.rma (metafor) Thank you for your attention and help! In this way I get the diamond coloured, but actually I would have the squares representing the values of the individual studies coloured. Is it somehow possible? Paola 2011/8/24 Viechtbauer Wolfgang (STAT) [hidden email] http://user/SendEmail.jtp?type=nodenode=3768683i=2 Thank you, Bernd, for looking into this. Yes, at the moment, the color of the summary estimate for models without moderators is hard-coded (as black). I didn't think people may want to change that. I guess I was wrong =) A dirty solution for the moment is to add: addpoly(dfs, efac=6, row=-1, col=red, border=red, annotate=F, mlab=) after the call to forest(). You will get a warning message (since the border argument gets passed to the text() function inside addpoly() and that's not a par for text), but you can just ignore that. Best, -- Wolfgang Viechtbauer Department of Psychiatry and Neuropsychology School for Mental Health and Neuroscience Maastricht University, P.O. Box 616 6200 MD Maastricht, The Netherlands Tel: +31 (43) 368-5248 Fax: +31 (43) 368-8689 Web: http://www.wvbauer.com -Original Message- From: Bernd Weiss [mailto:[hidden email]http://user/SendEmail.jtp?type=nodenode=3768683i=3] Sent: Wednesday, August 24, 2011 16:22 To: Paola Tellaroli Cc: [hidden email]http://user/SendEmail.jtp?type=nodenode=3768683i=4; [hidden email] http://user/SendEmail.jtp?type=nodenode=3768683i=5 Subject: Re: [R] Change color in forest.rma (metafor) Am 24.08.2011 07:50, schrieb Paola Tellaroli: My script is the following: library(metafor) yi-c(-0.1, 0.2, 0.3, 0.4) sei-c(0.4, 0.2, 0.6, 0.1) vi-sei^2 studi-c(A, B, C, D) eventi.c-c(10, 5, 7, 6) n.c-c(11, 34, 25, 20) eventi.a-c(2, 7, 6, 5) n.a-c(11, 35, 25, 15) dfs-rma(yi, vi, method=DL) dfs windows(height=6, width=10, pointsize=10) windowsFonts(B=windowsFont(Bookman Old Style
Re: [R] Change color in forest.rma (metafor)
I lied, that was not my last question: how can I add two arrows at the bottom with the words in favor of A / B? This is not specified in the pdf and with text I have the impression that I can't add text below the x-axis. 2011/8/26 Paola Tellaroli paola.tellar...@gmail.com Dear Prof. Viechtbauer, thank you so much for your help and kindness. Clearly graphs are the minor problem in our work, and the parameters and options that can vary in R are so many that it is obvious that you can not expect to change everything you want! Your suggestions are very helpuf, but I have one last question. I'm trying to copy the style of a forest plot that I've seen and I like (the one in the attached file, page 1034): can I do this in R? Best wishes, *Paola* 2011/8/25 Viechtbauer Wolfgang (STAT)-2 [via R] ml-node+3768683-1225159815-262...@n4.nabble.com The color of the squares is also currently hard coded. The thing is, there are so many different elements to a forest plot (squares, lines, polygons, text, axes, axis labels, etc.), if I would add arguments to set the color of each element, things would really get out of hand (as far as I am concerned, there are already too many arguments to begin with). I can think of one possibility: I could allow the col argument to accept a vector of colors and then apply the different elements of that vector to the various elements in the plot. Of course, there is also a limit to how far that can be taken. For example, what if somebody wants to have a different color for *one* of the squares and a different color for the other squares? Another possibility is to do some post-processing with other software. One can create the forest plot in R, save it for example as a postscript file, and the edit the plot in other software. Yes, I prefer it if I can create the plot in R and have it exactly the way I want it (without having to do any post-processing), but sometimes that may not be possible. Note that you can always add whatever you want to a plot created by the forest() function after it has been drawn. You can add text, lines, squares, polygons, whatever in any color you desire (e.g., with the text(), segments(), points(), polygon() functions). So, you could also just plot over the squares with: points(yi, 4:1, pch=15, col=red) To get rid of the black squares that are drawn by the forest function, add psize=0 as an argument in forest() (this will make the size of squares equal to 0, so essentially, they are invisible). If you want to make the size of the points inversely proportional to some function of the precision of the estimates, use points() together with the cex argument. For example: wi - 1/sqrt(vi) psize - wi/sum(wi) psize - (psize - min(psize)) / (max(psize) - min(psize)) psize - (psize * 1.0) + 0.5 points(yi, 4:1, pch=15, col=red, cex=psize) Best, Wolfgang -Original Message- From: Paola Tellaroli [mailto:[hidden email]http://user/SendEmail.jtp?type=nodenode=3768683i=0] Sent: Thursday, August 25, 2011 10:57 To: Viechtbauer Wolfgang (STAT) Cc: [hidden email]http://user/SendEmail.jtp?type=nodenode=3768683i=1; Bernd Weiss Subject: Re: [R] Change color in forest.rma (metafor) Thank you for your attention and help! In this way I get the diamond coloured, but actually I would have the squares representing the values of the individual studies coloured. Is it somehow possible? Paola 2011/8/24 Viechtbauer Wolfgang (STAT) [hidden email] http://user/SendEmail.jtp?type=nodenode=3768683i=2 Thank you, Bernd, for looking into this. Yes, at the moment, the color of the summary estimate for models without moderators is hard-coded (as black). I didn't think people may want to change that. I guess I was wrong =) A dirty solution for the moment is to add: addpoly(dfs, efac=6, row=-1, col=red, border=red, annotate=F, mlab=) after the call to forest(). You will get a warning message (since the border argument gets passed to the text() function inside addpoly() and that's not a par for text), but you can just ignore that. Best, -- Wolfgang Viechtbauer Department of Psychiatry and Neuropsychology School for Mental Health and Neuroscience Maastricht University, P.O. Box 616 6200 MD Maastricht, The Netherlands Tel: +31 (43) 368-5248 Fax: +31 (43) 368-8689 Web: http://www.wvbauer.com -Original Message- From: Bernd Weiss [mailto:[hidden email]http://user/SendEmail.jtp?type=nodenode=3768683i=3] Sent: Wednesday, August 24, 2011 16:22 To: Paola Tellaroli Cc: [hidden email]http://user/SendEmail.jtp?type=nodenode=3768683i=4; [hidden email] http://user/SendEmail.jtp?type=nodenode=3768683i=5 Subject: Re: [R] Change color in forest.rma (metafor) Am 24.08.2011 07:50, schrieb Paola Tellaroli: My script is the following: library(metafor) yi-c(-0.1, 0.2, 0.3, 0.4) sei-c(0.4, 0.2, 0.6
Re: [R] Change color in forest.rma (metafor)
Thank you for your attention and help!
In this way I get the diamond coloured, but actually I would have the
squares representing the values of the individual studies coloured. Is it
somehow possible?
*Paola*
2011/8/24 Viechtbauer Wolfgang (STAT)
wolfgang.viechtba...@maastrichtuniversity.nl
Thank you, Bernd, for looking into this.
Yes, at the moment, the color of the summary estimate for models without
moderators is hard-coded (as black). I didn't think people may want to
change that. I guess I was wrong =)
A dirty solution for the moment is to add:
addpoly(dfs, efac=6, row=-1, col=red, border=red, annotate=F, mlab=)
after the call to forest(). You will get a warning message (since the
border argument gets passed to the text() function inside addpoly() and
that's not a par for text), but you can just ignore that.
Best,
--
Wolfgang Viechtbauer
Department of Psychiatry and Neuropsychology
School for Mental Health and Neuroscience
Maastricht University, P.O. Box 616
6200 MD Maastricht, The Netherlands
Tel: +31 (43) 368-5248
Fax: +31 (43) 368-8689
Web: http://www.wvbauer.com
-Original Message-
From: Bernd Weiss [mailto:bernd.we...@uni-koeln.de]
Sent: Wednesday, August 24, 2011 16:22
To: Paola Tellaroli
Cc: w...@metafor-project.org; r-help@r-project.org
Subject: Re: [R] Change color in forest.rma (metafor)
Am 24.08.2011 07:50, schrieb Paola Tellaroli:
My script is the following:
library(metafor)
yi-c(-0.1, 0.2, 0.3, 0.4)
sei-c(0.4, 0.2, 0.6, 0.1)
vi-sei^2
studi-c(A, B, C, D)
eventi.c-c(10, 5, 7, 6)
n.c-c(11, 34, 25, 20)
eventi.a-c(2, 7, 6, 5)
n.a-c(11, 35, 25, 15)
dfs-rma(yi, vi, method=DL)
dfs
windows(height=6, width=10, pointsize=10)
windowsFonts(B=windowsFont(Bookman Old Style))
forest.rma(dfs, slab=studi, xlim=c(-15, 10), ilab=cbind(eventi.c, n.c,
eventi.a, n.a), ilab.xpos=c(-9.5, -8, -6, -4.5), cex=1.2, at=c(-2, -1,
0, 1,
2), family=B, xlab=Hazard Ratio (log scale), mlab=Random Effects
Model, efac=5, col=red, border=red)
text(-10, -1.3, paste(Heterogeneity: I-squared=,
paste(paste(round(dfs$I2,
2), %, sep=), paste(p, round(dfs$QEp, 4), sep==), sep=, ),
sep=), font=4, cex=1.2, family=B)
op-par(cex=1.2, font=2, family=B, oma=c(0.5, 0.5, 0.5, 0.5),
mar=c(0.5,
0.5, 0.5, 0.5))
text(x=c(-9.5, -8, -6, -4.5), 6, c(Events, N, Events, N),
cex=1.2 )
text(c(-8.7, -5.5, 8), 6.5, c(S, A, Log))
text(-15, 6, Trials, pos=4)
text(10, 6, Hazard Ratio [95% CI], pos=2)
par(op)
Even if I have specified col=red, border=red, color of squares
and
diamond rests black! Why?
As far as I know, col and border do only affect the fitted values
(diamonds), i.e. the FEM/REM estimators (see ?forest.rma: col:
character string specifying the name of a color to use for _the fitted_
values ('darkgray' by default).)
Furthermore, I had a quick look at the source code and it might be a
bug. If I replace in line 2770 the line
cex * efac), col = black, ...)
with
cex * efac), col = col, ...)
you can at least specify your own colour. Changing the border color
seems a bit more tricky...
However, Wolfgang Viechbauer (the package author) is always a very
responsive and helpful person and I suggest you better wait for his
answer.
Bernd
[[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] Change color in forest.rma (metafor)
My script is the following: library(metafor) yi-c(-0.1, 0.2, 0.3, 0.4) sei-c(0.4, 0.2, 0.6, 0.1) vi-sei^2 studi-c(A, B, C, D) eventi.c-c(10, 5, 7, 6) n.c-c(11, 34, 25, 20) eventi.a-c(2, 7, 6, 5) n.a-c(11, 35, 25, 15) dfs-rma(yi, vi, method=DL) dfs windows(height=6, width=10, pointsize=10) windowsFonts(B=windowsFont(Bookman Old Style)) forest.rma(dfs, slab=studi, xlim=c(-15, 10), ilab=cbind(eventi.c, n.c, eventi.a, n.a), ilab.xpos=c(-9.5, -8, -6, -4.5), cex=1.2, at=c(-2, -1, 0, 1, 2), family=B, xlab=Hazard Ratio (log scale), mlab=Random Effects Model, efac=5, col=red, border=red) text(-10, -1.3, paste(Heterogeneity: I-squared=, paste(paste(round(dfs$I2, 2), %, sep=), paste(p, round(dfs$QEp, 4), sep==), sep=, ), sep=), font=4, cex=1.2, family=B) op-par(cex=1.2, font=2, family=B, oma=c(0.5, 0.5, 0.5, 0.5), mar=c(0.5, 0.5, 0.5, 0.5)) text(x=c(-9.5, -8, -6, -4.5), 6, c(Events, N, Events, N), cex=1.2 ) text(c(-8.7, -5.5, 8), 6.5, c(S, A, Log)) text(-15, 6, Trials, pos=4) text(10, 6, Hazard Ratio [95% CI], pos=2) par(op) Even if I have specified col=red, border=red, color of squares and diamond rests black! Why? Thanks, Paola -- View this message in context: http://r.789695.n4.nabble.com/Change-color-in-forest-rma-metafor-tp3765090p3765090.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] help with nls fitting
Dear all, I'm trying to fit the following function slope_pp3_mrna = ( (k3 * v3_K_d *p1^v3_h) / ( (v3_Kd^v3_h) + p2^v3_h ) ) * ( 1/(1 + (p2/v4_Kd)^v4_h) ) - pp3_mrna to this experimental data in the datafraeme Data_pp3_mrna (see it at the end of this e-mail) I'm using the nls function in the following code. IN the last step of the fit fm_pp3_mrna_4, when I add to the funziont the paramter v4_Kd something goes wrong, and I reeive this message Error in numericDeriv(form[[3L]], names(ind), env) : Missing value or an infinity produced when evaluating the model What could be the error? I tried with differnet intial values of v4_Kd, but I did not fix the problem. Here below the code. Thanks in advance, Paola. Data_pp3_mrna - data.frame( p1 = protein_1, p2 = protein_2, pp3_mrna = protein_3_mrna, slope_pp3_mrna = stinemanSlopes(times, protein_3_mrna) ) one_par - nls(slope_pp3_mrna ~ ( (k3 * p1) / ( (1) + p1 ) ) * ( 1/(1 + (p2)) ) - pp3_mrna, data=Data_pp3_mrna, start=list( k3=0.1 )) summary(one_par) fm_pp3_mrna_1 - nls(slope_pp3_mrna ~ ((k3 * v3_Kd *p1) / ( (v3_Kd) + p1 ) ) * ( 1/(1 + (p2)) ) - pp3_mrna, data=Data_pp3_mrna, start=list(k3=24, v3_Kd=1 )) summary(fm_pp3_mrna_1) fm_pp3_mrna_2 - nls(slope_pp3_mrna ~ ((k3 * v3_Kd *p1^v3_h) / ( (v3_Kd)^v3_h + p1^v3_h ) ) * ( 1/(1 + (p2)) ) - pp3_mrna, control = list(maxiter = 500), data=Data_pp3_mrna, start=list(k3=69, v3_Kd=0.3238, v3_h=1 )) summary(fm_pp3_mrna_2) fm_pp3_mrna_3 - nls(slope_pp3_mrna ~ ((k3 * v3_Kd *p1^v3_h) / ( (v3_Kd)^v3_h + p1^v3_h ) ) * ( 1/(1 + (p2)^v4_h) ) - pp3_mrna, control = list(maxiter = 500), data=Data_pp3_mrna, start=list(k3=37.451, v3_Kd=0.59, v3_h=2.013, v4_h=0.01 )) summary(fm_pp3_mrna_3) fm_pp3_mrna_4 - nls(slope_pp3_mrna ~ ((k3 * v3_Kd *p1^v3_h) / ( (v3_Kd)^v3_h + p1^v3_h ) ) * ( 1/(1 + (p2/v4_Kd)^v4_h) ) - pp3_mrna, control = list(maxiter = 500), data=Data_pp3_mrna, start=list(k3=56.2823, v3_Kd=0.3366, v3_h=1.8040, v4_Kd=0.03, v4_h=0.7693 )) Here the data. Data_pp3_mrna p1 p2 pp3_mrna slope_pp3_mrna 1 1.006 0.921 0.041 8.63741887 2 2.235 2.047 2.9069031 2.82619343 3 3.744 3.937 4.052 0.84354113 4 4.222 9.340 4.3237353 0.47577213 5 9.022 14.609 4.531-0.03940131 6 11.326 22.765 3.510-2.0420 7 6.899 17.852 2.489-1.86822481 8 10.709 27.777 1.6222048-1.55625973 9 14.084 27.785 0.911-0.48800514 10 14.922 23.613 0.826-0.1700 11 14.340 18.422 0.741-0.22560156 12 13.066 24.085 0.599-0.2840 13 17.553 18.594 0.457-0.13847372 14 14.803 16.831 0.4550965 0.03588624 15 11.945 14.495 0.493 0.09536674 16 11.427 12.458 0.5505361 0.15549062 17 11.556 9.082 0.649 0.49638596 18 20.107 9.987 1.2486525 1.36871828 19 15.999 10.305 2.059 1.87868197 20 16.094 5.793 3.2285000 2.3390 21 11.752 6.944 4.398 0.80395869 22 15.841 5.575 4.651 0.5060 23 12.601 5.221 4.904 0.80270128 24 13.598 2.872 5.819 1.8300 25 13.879 2.883 6.734-0.03571884 26 16.270 2.213 6.4135000-0.6410 27 17.176 3.381 6.093-0.33913760 28 12.332 2.781 6.032-0.1220 29 12.373 3.073 5.971-0.41224459 30 14.781 2.948 5.489-0.9640 31 17.578 3.953 5.007-0.68258311 32 18.865 2.279 4.7568901-0.39425557 33 14.735 2.806 4.606 0.13264481 34 16.160 1.676 4.987 0.7620 35 13.416 2.478 5.368 0.63914248 36 13.864 1.394 5.6374239 0.45419096 37 16.219 2.299 5.827-0.12540453 38 13.249 1.457 5.2505000-1.1530 39 14.445 2.325 4.674-0.09880849 40 13.210 2.230 5.5576966 2.17200142 41 12.358 2.116 8.94711.38521228 -- *Paola Lecca, PhD* *The Microsoft Research - University of Trento* *Centre for Computational and Systems Biology* *Piazza Manci 17 38123 Povo/Trento, Italy* *Phome: +39 0461282843* *Fax: +39 0461282814* [[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] fit a 2-variables function to data
Dearl all, I have to fit a function y = f(x1, x2) to data experiemntal data describing the measured behavior of y. x1 and x2 are the independent variables. Could you suggest me wich R package can I use for this purpose? Thanks, Paola. -- *Paola Lecca, PhD* *The Microsoft Research - University of Trento* *Centre for Computational and Systems Biology* *Piazza Manci 17 38123 Povo/Trento, Italy* *Phome: +39 0461282843* *Fax: +39 0461282814* [[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] problemsn in using nls
Dear all, I tried to use nls, but I got the following error Error in numericDeriv(form[[3L]], names(ind), env) : Missing value or an infinity produced when evaluating the model Any suggestion? Thanks, Paola. The code I wrote is Data_pp2_mrna - data.frame( p1 = protein_1, p6 = protein_6, pp2_mrna = protein_2_mrna, slope_pp2_mrna = stinemanSlopes(times, protein_2_mrna) ) fm_pp2_mrna - nls(slope_pp2_mrna ~ ( (k1 * v2_Kd *p1^v2_h) / ( (v2_Kd^v2_h) + p1^v2_h ) ) * ( 1/(1 + (p6/v5_Kd)^v5_h) ) - pp2_mrna, data = Data_pp2_mrna, start = list(k1 = 1, v2_Kd = 1, v2_h = 1, v5_Kd = 1, v5_h = 1 )) The data are as follows Data_pp2_mrna p1 p6 pp2_mrna slope_pp2_mrna 1 1.006 1.234 0.000 1.0183976996 2 2.235 0.693 0.5565718 1.2167043185 3 3.744 0.451 1.230 1.6962541888 4 4.222 0.441 2.620 2.78 5 9.022 0.523 4.010 0.1298238635 6 11.326 0.845 3.9179535 -0.2048280861 7 6.899 0.674 3.805 -0.4595669210 8 10.709 1.369 3.386 -0.838000 9 14.084 1.646 2.967 -0.5032137310 10 14.922 2.561 2.822 -0.29 11 14.340 2.265 2.677 -0.300152 12 13.066 4.832 2.4634792 -0.4859794757 13 17.553 5.651 2.188 -0.7700435122 14 14.803 6.418 1.6045000 -1.167000 15 11.945 10.343 1.021 -0.3160241819 16 11.427 9.415 1.0435000 0.045000 17 11.556 12.610 1.066 -0.1700195282 18 20.107 16.171 0.8545000 -0.423000 19 15.999 13.189 0.643 -0.1721421868 20 16.094 17.022 0.6635000 0.041000 21 11.752 14.723 0.684 -0.0911220974 22 15.841 15.887 0.569 -0.23 23 12.601 15.604 0.454 -0.0022539547 24 13.598 18.202 0.5665000 0.225000 25 13.879 20.821 0.679 -0.0155566571 26 16.270 18.040 0.549 -0.26 27 17.176 16.015 0.419 -0.1195632728 28 12.332 18.429 0.425 0.012000 29 12.373 17.723 0.431 -0.1118527520 30 14.781 23.100 0.3095000 -0.243000 31 17.578 15.396 0.188 0.0004440266 32 18.865 21.408 0.310 0.244000 33 14.735 14.518 0.432 0.0473462057 34 16.160 21.351 0.361 -0.142000 35 13.416 16.689 0.290 0.0424156026 36 13.864 17.015 0.4065000 0.233000 37 16.219 21.776 0.523 0.1567965535 38 13.249 17.371 0.5650035 0.0498847737 39 14.445 19.218 0.573 -0.0435396634 40 13.210 21.309 0.5211514 -0.1713939796 41 12.358 24.966 0.400 -0.3132118164 -- *Paola Lecca, PhD* *The Microsoft Research - University of Trento* *Centre for Computational and Systems Biology* *Piazza Manci 17 38123 Povo/Trento, Italy* *Phome: +39 0461282843* *Fax: +39 0461282814* [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] use of modMCMC
Dear all, I used modFit of the package FME to fit a set of ODE to a ste of eperiemntal data. The summary of this fit give me the following error summary(Fit) Residual standard error: 984.1 on 452 degrees of freedom Error in cov2cor(x$cov.unscaled) : 'V' is not a square numeric matrix In addition: Warning message: In summary.modFit(Fit) : Cannot estimate covariance; system is singular This is due becasue the Hessian matrix has all the entries equal to 0. In these cases, on the help page of modFit, it is suggested to use modMCMC to generate new sets of parameters. modMCMC performs a Markov Chain Monte Carlo simulation. I do not understand very well how modMCMC can be used in a context of parameter estimation. Could someone help me in understanding the use of this function and its utility for parameter fiting? Thank you very much in advance, Paola. -- *Paola Lecca, PhD* *The Microsoft Research - University of Trento* *Centre for Computational and Systems Biology* *Piazza Manci 17 38123 Povo/Trento, Italy* *Phome: +39 0461282843* *Fax: +39 0461282814* [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] use of modMCMC
Dear all, I used modFit of the package FME to fit a set of ODE to a ste of experimental data. The summary of this fit give me the following error summary(Fit) Residual standard error: 984.1 on 452 degrees of freedom Error in cov2cor(x$cov.unscaled) : 'V' is not a square numeric matrix In addition: Warning message: In summary.modFit(Fit) : Cannot estimate covariance; system is singular This is due becasue the Hessian matrix has all the entries equal to 0 in my system. In these cases, on the help page of modFit, it is suggested to use modMCMC to generate new sets of parameters. modMCMC performs a Markov Chain Monte Carlo simulation. I do not understand very well how modMCMC can be used in a context of parameter estimation. Could someone help me in understanding the use of this function and its utility for parameter fitting? Thank you very much in advance, Paola. -- *Paola Lecca, PhD* *The Microsoft Research - University of Trento* *Centre for Computational and Systems Biology* *Piazza Manci 17 38123 Povo/Trento, Italy* *Phome: +39 0461282843* *Fax: +39 0461282814* [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Help with modFit of FME package
Dear R users, I'm trying to fit a set an ODE to an experimental time series. In the attachment you find the R code I wrote using modFit and modCost of FME package and the file of the time series. When I run summary(Fit) I obtain this error message, and the values of the parameters are equal to the initial guesses I gave to them. The problem is not due to the fact that I have only one equation (I tried also with more equations, but I still obtain this error). I would appreciate if someone could help me in understanding the reason of the error and in fixing it. Thanks for your attention, Paola Lecca. Here the error: summary(Fit) Parameters: Estimate Std. Error t value Pr(|t|) pro1_strength1 NA NA NA Residual standard error: 2.124 on 10 degrees of freedom Error in cov2cor(x$cov.unscaled) : 'V' is not a square numeric matrix In addition: Warning message: In summary.modFit(Fit) : Cannot estimate covariance; system is singular __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- *Paola Lecca, PhD* *The Microsoft Research - University of Trento* *Centre for Computational and Systems Biology* *Piazza Manci 17 38123 Povo/Trento, Italy* *Phome: +39 0461282843* *Fax: +39 0461282814* timepp1_mrna 0 0 2 2.754 4 2.958 6 4.058 8 3.41 10 3.459 12 2.453 14 1.234 16 2.385 18 3.691 20 3.252 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 with modFit of FME package 2
* Apologies for multiple posting *
I attached to my previous e-mail a .r file, and it was not permitted by the
rules of the mailing lis. Again, please receive my sincere apologies for
this.
I re-send again the e-mail with .txt attachemnt in the hope someone an help
me to solve my problem.
I'm trying to fit a set an ODE to an experimental time series. In the
attachment you find the R code I wrote using modFit and modCost of FME
package and the file of the time series.
When I run summary(Fit) I obtain this error message, and the values of the
parameters are equal to the initial guesses I gave to them.
The problem is not due to the fact that I have only one equation (I tried
also with more equations, but I still obtain this error).
I would appreciate if someone could help me in understanding the reason of
the error and in fixing it.
Thanks for your attention,
Paola Lecca.
Here the error:
summary(Fit)
Parameters:
Estimate Std. Error t value Pr(|t|)
pro1_strength1 NA NA NA
Residual standard error: 2.124 on 10 degrees of freedom
Error in cov2cor(x$cov.unscaled) : 'V' is not a square numeric matrix
In addition: Warning message:
In summary.modFit(Fit) : Cannot estimate covariance; system is singular
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
*Paola Lecca, PhD*
*The Microsoft Research - University of Trento*
*Centre for Computational and Systems Biology*
*Piazza Manci 17 38123 Povo/Trento, Italy*
*Phome: +39 0461282843*
*Fax: +39 0461282814*
timepp1_mrna
0 0
2 2.754
4 2.958
6 4.058
8 3.41
10 3.459
12 2.453
14 1.234
16 2.385
18 3.691
20 3.252
require(deSolve)
require(FME)
##
# PART 1
#
##
# Differential equations
model_1_part_1 - function(t, S, parameters)
{
with(as.list(parameters), {
#
cod1 = pro1_strength
#
pp1_mrna_degradation_rate - 1
###
#
v1 = cod1
v2 = pp1_mrna_degradation_rate * S[1]
#
#
#
dS1 = v1 - v2
#
#
list(c(dS1))
})
}
# Parameters
parms_part_1 - c(pro1_strength = 1.0)
# Initial values of the species concentration
S - c(pp1_mrna = 0)
times - seq(0, 20, by = 2)
# Solve the system
ode_solutions_part_1 - ode(S, times, model_1_part_1, parms = parms_part_1)
ode_solutions_part_1
summary(ode_solutions_part_1)
## Default plot method
plot(ode_solutions_part_1)
# Estimate of the parameters
experiment - read.table(./wild_pp1_mrna.txt, header=TRUE)
rw - dim(experiment)[1]
names - array(, rw)
for (i in 1:rw)
{
names[i] - pp1_mrna
}
names
observed_data_part_1 - data.frame(name = names,
time = experiment[,1], val = experiment[,2])
observed_data_part_1
ode_solutions_part_1
Cost_function - function (pars) {
out - ode_solutions_part_1
cost - modCost(model = out, obs = observed_data_part_1, y = val)
cost
}
Cost_function(parms)
# Fit the model to the observed data
Fit - modFit(f = Cost_function, p = parms_part_1)
Fit
# Summary of the fit
summary(Fit)
# Model coefficients
coef(Fit)
# Deviance of the fit
deviance(Fit)__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 with modFit of FME package
Dear R users, I’m trying to fit a set an ODE to an experimental time series. In the attachment you find the R code I wrote using modFit and modCost of FME package and the file of the time series. When I run summary(Fit) I obtain this error message, and the values of the parameters are equal to the initial guesses I gave to them. The problem is not due to the fact that I have only one equation (I tried also with more equations, but I still obtain this error). I would appreciate if someone could help me in understanding the reason of the error and in fixing it. Thanks for your attention, Paola. Here the error: summary(Fit) Parameters: Estimate Std. Error t value Pr(|t|) pro1_strength1 NA NA NA Residual standard error: 2.124 on 10 degrees of freedom Error in cov2cor(x$cov.unscaled) : 'V' is not a square numeric matrix In addition: Warning message: In summary.modFit(Fit) : Cannot estimate covariance; system is singular -- *Paola Lecca, PhD* *The Microsoft Research - University of Trento* *Centre for Computational and Systems Biology* *Piazza Manci 17 38123 Povo/Trento, Italy* *Phome: +39 0461282843* *Fax: +39 0461282814* timepp1_mrna 0 0 2 2.754 4 2.958 6 4.058 8 3.41 10 3.459 12 2.453 14 1.234 16 2.385 18 3.691 20 3.252 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 with modFit of FME package
Dear R users, I'm trying to fit a set an ODE to an experimental time series. In the attachment you find the R code I wrote using modFit and modCost of FME package and the file of the time series. When I run summary(Fit) I obtain this error message, and the values of the parameters are equal to the initial guesses I gave to them. The problem is not due to the fact that I have only one equation (I tried also with more equations, but I still obtain this error). I would appreciate if someone could help me in understanding the reason of the error and in fixing it. Thanks for your attention, Paola Lecca. Here the error: summary(Fit) Parameters: Estimate Std. Error t value Pr(|t|) pro1_strength1 NA NA NA Residual standard error: 2.124 on 10 degrees of freedom Error in cov2cor(x$cov.unscaled) : 'V' is not a square numeric matrix In addition: Warning message: In summary.modFit(Fit) : Cannot estimate covariance; system is singular timepp1_mrna 0 0 2 2.754 4 2.958 6 4.058 8 3.41 10 3.459 12 2.453 14 1.234 16 2.385 18 3.691 20 3.252 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] frequency table help
Hi, I have a dataset like this: Specieslength (cm) A 12.4 B 45 A 34.6 C 73 C 24.5 D 4.5 .. I'm trying to obtain a barplot with the class length in x (fixed classes, 5 cm) and the number of species in y, but using just a barplot is not the case (I think). So, maybe the best way is to first obtain a frequency table and then to plot the graph, but I managed only to obtain a frequency table with 5 cm of class length but the frequency are the numbers of individuals by class and this is not what I need. Do you know how I can do? Thank you in advance. Paola [[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.
