Re: [R] Computing time for matrix addition or subtraction
On Sun, 12 Nov 2006, YONGWAN CHUN wrote: I wonder by chance if there is a way to reduce computing time for matrix addition or subtraction. With a lot of iterations, it would be helpful to reduce a little amount time. Yes, by making use of an optimized BLAS: see the R-admin manual. On my (dual CPU) system this reduced your example from 36 to 6 seconds. BTW, it is the calculation of PP that is taking the most of time, not as in your subject line. Simple example is as below n - 2000 P - matrix(rnorm(n*n),n,n) PP - P %*% P M - diag(n) - P R - M + t(M) - diag(n) + PP I would like to reduce time in calculating R. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ 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] Computing time for matrix addition or subtraction
not concerning your subject line, but function crossprod may be useful, too Matthias - original message Subject: Re: [R] Computing time for matrix addition or subtraction Sent: Mon, 13 Nov 2006 From: Prof Brian Ripley[EMAIL PROTECTED] On Sun, 12 Nov 2006, YONGWAN CHUN wrote: I wonder by chance if there is a way to reduce computing time for matrix addition or subtraction. With a lot of iterations, it would be helpful to reduce a little amount time. Yes, by making use of an optimized BLAS: see the R-admin manual. On my (dual CPU) system this reduced your example from 36 to 6 seconds. BTW, it is the calculation of PP that is taking the most of time, not as in your subject line. Simple example is as below n - 2000 P - matrix(rnorm(n*n),n,n) PP - P %*% P M - diag(n) - P R - M + t(M) - diag(n) + PP I would like to reduce time in calculating R. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ 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. --- original message end __ 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] Need help in waveslim package: imodwt and universal.thresh.modwt
Hi: I have encountered problems with imodwt and universal.thresh.modwt and cannot find any reference in R Search. Hope someone can give me some ideas: Starting with modwt.la8 - modwt(xdata, la8, n.level=6) -- this seems to work fine (1) ydata - imodwt(modwt.la8) will always give ydata as numeric(0) (no values) instead of being a time series data with the same lenght as xdata. (2) thred.la8 - universal.thresh.modwt(modwt.la8, max.level=4, hard=F) will give me error at: abs(wc.fine) - cannot contain non-numeric field. Any help is much appreciated. Regards ___ YM - Â÷½u°T®§ ´Nºâ§A¨S¦³¤Wºô¡A§AªºªB¤Í¤´¥i¥H¯d¤U°T®§µ¹§A¡A·í§A¤Wºô®É´N¯à¥ß§Y¬Ý¨ì¡A¥ô¦ó»¡¸Ü³£ÉN¨«¥¢¡C [[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] doubt on Tukey HSD
dear all effect of a, b and c on d, total 48 comparisons, got one anova result in model1=aov(d~A*B*C). can we get all the result in one command? or can we interpret the whole comparisons from this result? how it work in Tukey HSD? jose Access over 1 million songs. __ 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] Problems with metaMDS from vegan
Hello Spencer, 1. Have you tried to develop an example so small and simple that it can be provided with a few lines of code in an email like this (as suggested in the trailer to every email distributed by r-help)? I did, but due to the simplicity of the example in my former posting you may have overlooked it. You should not criticize the maintainer for not replying to your earlier email: I did not see in your email any information that would allow anyone to solve your problem unless s/he just happened to be unbelievably lucky. Well, the mentioning of my post to the maintainer was not intended as criticism, but as justification for the public posting itself. Please include 'sessionInfo()' This does not reveal any additional information compared to that given in my former post, but anyhow (I noticed a new version of Vegan on Jari's HP and updated accordingly in the hope that it might cure the effect - which it didn't): On the machine where the IDENTICAL CODE AND DATA is defunct: sessionInfo() R version 2.2.1, 2005-12-20, i486-pc-linux-gnu attached base packages: [1] methods stats graphics grDevices utils datasets [7] base other attached packages: MASSvegan 7.2-24 1.9-3 On the machine where the code and data runs smoothly: sessionInfo() R version 2.2.1, 2005-12-20, i486-pc-linux-gnu attached base packages: [1] methods stats graphics grDevices utils datasets [7] base other attached packages: MASSvegan 7.2-24 1.6-10 (The third machine [Win XP with R 2.4.0 and Vegan 1.8.3-x], show the defunct symptoms, is not available for testing right now) 2. Have you tried 'traceback()' after your failed call to 'metaMDS'? traceback() on the defunct installation yields: traceback() 6: stop('x' must be an array of at least two dimensions) 5: rowSums(x, na.rm = TRUE) 4: any(rowSums(x, na.rm = TRUE) == 0) 3: vegdist(comm, method = distance, ...) 2: metaMDSdist(comm, distance = distance, autotransform = autotransform, noshare = noshare, trace = trace, commname = commname, ...) 1: metaMDS(distab.dist, distance = bray, k, trymax = 50) With luck, this will identify the function actually generating the error. Then you can list that function and try find rowSums(x, a.rm = TRUE). Then hunt for x in that code. That rowSums() code is identical on both installations, but the higher layers of invocation differ in differing amounts. I checked with mgdiff. The function that includes this call to 'rowSums' expects 'x' to be an array of at least two dimensions. Either you pass it something that does NOT meet that description or someplace x might be created from selecting rows or columns of a larger array and either 0 or 1 rows were selected, thereby reducing the dimension of the array to 0 or 1. Yes, I guessed such a dependency already, but do not see the exact place where to dig into. The major code differences are in the vegdist() and the metaMDSdist() functions but I cannot point my finger to a certain code part where I estimate it to break - sorry! 3. If 1 and 2 fail to provide enlightenment, I would then try the 'debug' functions (provided it was important enough for me to do so, of course). ... This is my next step to try as the first two do not seem to yield further insights for me, but I just wanted to return a timely reply to the other points you mentioned in your post and with a certain hope that it might give some longer term R users some ideas to communicate some ideas. Thanks for your involvement, Peter __ 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] A printing macro
I am exploring the result of clustering a large multivariate data set into a number of groups, represented, say, by a factor G. I wrote a function to see how categorical variables vary between groups: ddisp - function(dvar) { + csqt - chisq.test(G,dvar) + print(csqt$statistic) + print(csqt$observed) + print(round(csqt$expected)) + round(csqt$residuals) + } x - ceiling(4*runif(100)) G - gl(4,1,100) ddisp(x) X-squared 6.235645 dvar G1 2 3 4 1 10 5 5 5 2 6 9 5 5 3 8 6 5 6 4 7 4 4 10 dvar G 1 2 3 4 1 8 6 5 6 2 8 6 5 6 3 8 6 5 6 4 8 6 5 6 dvar G1 2 3 4 1 1 0 0 -1 2 -1 1 0 -1 3 0 0 0 0 4 0 -1 0 1 Warning message: Chi-squared approximation may be incorrect in: chisq.test(G, dvar) As I need to apply this function to a large number of variables x it would be helpful if the function printed x rather than the formal argument dvar. I have a vague idea that things like deparse() and substitute() will come into the solution but I have not yet come up with the right incantation. Any help appreciated! Murray Jorgensen -- Dr Murray Jorgensen http://www.stats.waikato.ac.nz/Staff/maj.html Department of Statistics, University of Waikato, Hamilton, New Zealand Email: [EMAIL PROTECTED]Fax 7 838 4155 Phone +64 7 838 4773 wkHome +64 7 825 0441Mobile 021 1395 862 __ 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] stepAIC for overdispersed Poisson
I am wondering if stepAIC in the MASS library may be used for model selection in an overdispersed Poisson situation. What I thought of doing was to get an estimate of the overdispersion parameter phi from fitting a model with all or most of the available predictors (we have a large number of observations so this should not be problematical) and then use stepAIC with scale = phi. Should this be OK? Murray Jorgensen -- Dr Murray Jorgensen http://www.stats.waikato.ac.nz/Staff/maj.html Department of Statistics, University of Waikato, Hamilton, New Zealand Email: [EMAIL PROTECTED]Fax 7 838 4155 Phone +64 7 838 4773 wkHome +64 7 825 0441Mobile 021 1395 862 __ 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] Problems with metaMDS: Resolution!
The problem is solved: With the help of the Vegan library maintainer the error was found immediately. It was a misinterpretation/overlooking of the provided documentation. The problem lies here: meta - metaMDS(distab.dist, distance=bray, k, trymax=50) The provided primary data matrix - distab.dist - was in my case a distance matrix, not a original data set matrix as demanded in the documentation. As Jari Oksanen told me, several other users seem to have been trapped by the same pragmatic error also, which are reported in his newer versions of the Vegan library only. At least version 1.6-10 ran smoothly without reporting an error, thus leaving the user with the faulty notion of a correctly running program. It seems that the pragmatic error of not only me stemmed from the fact that the frequently used methods cmdscale and isoMDS from the original R distribution scope require indeed *distance matrices* as their first input value. So, switching from them to the Vegan metaMDS, makes one easily overlook the important change in required data structure. The error condition that the newer versions of Vegan now returns helps to detect this faulty data provision. Kind regards and many thanks to the list and the maintainer, Peter __ 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] legend positioning problems
Hello, I have several Barplots I want to plot. I also want to include a legend. Since my Barplots are very different I ve dicided to put the legend in the top left corner. Unfortunately, sometimes there is a part of a bar just below the legend. This makes it difficult to see the legend itself or the barplot. All I found out so far, is to make the backgroundcolor of the legend white or so, but I would couldn t see the bar anymore. Is there a way, to position the legend excel-like on the outside the plot on the side? or Is there a possibility to use sort of an alpha-factor to make the background white but slightly invisible so that there is a contrast? I hope you can help me out and thank you in advance! #Example: barplot(VADeaths) legend(topleft, c(1x, 2x, 3x, 4x, 5x), bg=white)) Markus Schweitzer - http://www.pokertips.tk __ 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] stepAIC for overdispersed Poisson
On Mon, 13 Nov 2006, Murray Jorgensen wrote: I am wondering if stepAIC in the MASS library may be used for model selection in an overdispersed Poisson situation. What I thought of doing was to get an estimate of the overdispersion parameter phi from fitting a model with all or most of the available predictors (we have a large number of observations so this should not be problematical) and then use stepAIC with scale = phi. Should this be OK? Well no, as that quasi-Poisson model does not have an AIC. Remember AIC assumes maximum likelihood fitting, and you don't have a likelihood here (even for fixed phi). The problem is that an 'overdispersed Poisson situation' is a situation, not a model (in the usual sense of a probability measure over possible outcomes). You could use a negative binomial GLM (with fixed theta). -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ 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] change \baselinestretch in Soutput ?
Hi there, Is it possible to change \baselinestretch in Soutput ? What should \DefineVerbatimEnvironment{Soutput}{Verbatim}{fontsize=\footnotesize} look like, if it is possible to effect the change here? Google did not help here. best Christian -- Dr. Christian W. Hoffmann, Swiss Federal Research Institute WSL Zuercherstrasse 111, CH-8903 Birmensdorf, Switzerland Tel +41-44-7392-277 (office), -111(exchange), -215 (fax) [EMAIL PROTECTED], www.wsl.ch/staff/christian.hoffmann __ 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 in waveslim package: imodwt and universal.thresh.modwt
The first two commands worked fine to me. I used xdata = rnorm(128). To let us help you more, specify your R version, some data, OS etc, as usually asked. Rogerio Porto. -- Cabeçalho original --- De: [EMAIL PROTECTED] Para: r-help@stat.math.ethz.ch Cópia: Data: Sun, 12 Nov 2006 16:02:45 +0800 (CST) Assunto: [R] Need help in waveslim package: imodwt and universal.thresh.modwt Hi: I have encountered problems with imodwt and universal.thresh.modwt and cannot find any reference in R Search. Hope someone can give me some ideas: Starting with modwt.la8 - modwt(xdata, la8, n.level=6) -- this seems to work fine (1) ydata - imodwt(modwt.la8) will always give ydata as numeric(0) (no values) instead of being a time series data with the same lenght as xdata. (2) thred.la8 - universal.thresh.modwt(modwt.la8, max.level=4, hard=F) will give me error at: abs(wc.fine) - cannot contain non-numeric field. Any help is much appreciated. Regards ___ YM - Â÷½u°T®§ ´Nºâ§A¨S¦³¤Wºô¡A§AªºªB¤Í¤´¥i¥H¯d¤U°T®§µ¹§A¡A·í§A¤Wºô®É´N¯à¥ß§Y¬Ý¨ì¡A¥ô¦ó»¡¸Ü³£ÉN¨«¥¢¡C [[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] Heatmap.2 in gplots
I used the function heatmap.2 some times ago, and it worked as intended (for me...). I tried to reuse my old R program and I encounter some trouble. I suppose that something was changed in newer versions. One of the examples given in the help page shows the same problem : -- this part wa ## For variable clustering, rather use distance based on cor(): data(USJudgeRatings) symnum( cU - cor(USJudgeRatings) ) hU - heatmap.2(cU, Rowv=FALSE, symm=TRUE, col=topo.colors(16), distfun=function(c) as.dist(1 - c), trace=none) ## The Correlation matrix with same reordering: hM - format(round(cU[hU[[1]], hU[[2]]], 2)) hM # now with the correlation matrix on the plot itself heatmap.2(cU, Rowv=FALSE, symm=TRUE, col=rev(heat.colors(16)), distfun=function(c) as.dist(1 - c), trace=none, cellnote=hM) When running this correlation values are supposed to appear on the plot. But, as one can see in this example, the same value is not associated with the same color ( sometimes 1.00 is red or pale yellow). Looks like some reordering done on colors was'nt done on cellnotes. [[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] GESCHAEFTSVORCHLAG
Sehr geehrter Geschaefts-Vorschlag! Sicher sind Sie verwundert, diese Nachricht von jemandem zu erhalten, den sie nicht personlich kennen. Der Grund weshalb ich mich an Sie wende ist das ich Ehret Olds bin, der erstgeborene Sohn von Martin Olds, einem der bekanntesten schwarzen Farmer in Zimbabwe. Er wurde vor kurzem im Laufe des Landstreites in meinem Land ermordet. Ich erhielt Ihre Adresse durch das Internet und habe mich entschieden mich an Sie zu wenden. Bevor mein Vater starb, nahm er mich mit nach Johannesburg um dort einen Betrag von 8,5 Millionen US$ (Acht Millionen funfhunderttausend US Dollar) bei einer privaten Sicherheitsfirma zu deponieren. Er hat die drohende Gefahr geahnt und daher diesen Betrag in einer Kassette dort hinterlegt, die er als normales Wertgut deklarierte um keinen Verdacht zu erwecken und ausserdem die Kosten fuer die Hinterlegung gering zu halten. Niemand weis also von dem Bargeld in dieser Kassette. Ursprunglich war dieser Betrag fuer neue Maschinen und Chemikalien für seine Farm gedacht, ausserdem wollte er neues Farmland in Schwerz erwerben und erschliesen. Das Landproblem entstand als der Präsident von Zimbabwe, Mr. Robert Mugabe, seine neue Landreform durchsetzte bei der nur die reichen weisen Farmer und einige wenige schwarze Farmer Vorteile verschafft wurden, fuer alle anderen war die Existenz nun bedroht. Daraus resultierten dann schlielich die Unruhen die sich ausbreiteten und denen viele Menschen zum Opfer fielen und ihr Leben lassen mussten. Leider ist mein Vater ebenfalls eines dieser Opfer das sein Leben verlor. Vor diesem Hintergrund bin ich mit der Familie unter lebensgefährlichen Umständen aus Zimbabwe nach Holland geflohen, hier haben wir politisches Asyl beantragt. Das in Südafrika zurückgelassene Geld wollen wir nun der Sicherheit wegen auf ein Konto hier übertragen. Das Gesetz in Holland verbietet es Menschen die politisches Asyl beantragt haben während dieses Verfahrens ein Konto zu eröffnen oder irgendwelche Transfers/Transaktionen die Über die hollandische Grenze hinausgehen abzuwickeln. Als der erstgeborene Sohn bin ich nun für meine ganze Familie veranwortlich und habe die Rolle meines Vaters als Bewacher und Beschitzer Übernommen. Nun stehe ich vor dem Problem dieses Geld ohne Wissen der afrikanischen Regierung hierher zu transferieren, ansonsten wurde man uns den gesamten Betrag enteignen, dieses Geld ist aber alles das uns noch geblieben ist. Die Südafrikanische Regierung unterstützt wohl die Regierung Zimbabwes, so dass auch dort nichts von meinem Vorhaben bekannt werden darf. Als Geschäftsmann suche ich nun nach einem Partner dem ich voll und ganz vertrauen kann, dem ich meine und auch die Zukunft meiner ganzen Familie anvertrauen kann. Ich möchte Sie wissen lassen, dass ein Plan völlig ohne Risiko ist, falls Sie mir und meiner Familie helfen wollen. Alles worum ich Sie bitte, ist eine Vereinbarung mit der hollandischen Niederlassung der Sicherheitsfirma zu treffen, das hinterlegte Wertgut auszuhändigen. Ich meinerseits habe bereits Anweisung gegeben die Kassette von Südafrika nach Holland zu senden. Vorher allerdings müssen noch einige wichtige Formalitäten erledigt werden, wie die Änderung des Begunstigten an diesem Wertgut. Weiterhin habe ich vor diese Geld nach Erhalt möglichst gut und gewinnbringend zu investieren. Ich möchte Ihnen zwei Vorschläge unterbreiten.Erstens biete ich Ihnen einen Teil der Gesamtsumme an, wenn Sie bereit wären, Ihr eigenes Konto für diese Transaktion Verfügung zu stellen.Oder aber Sie sind an einer Partnerschaft mit mir interessiert um diese Summe profitabel in Ihrem Land zu investieren.Egal welcher Vorschlag Ihnen mehr zusagt, zögern Sie nicht mir Ihre Entscheidung mitzuteilen. Fuer alle entstehenden Unkosten habe ich 5% des Gesamtbetrages ingeplant. Sollten Sie nicht an einer Partnerschaft mit Mir interessieren sein, biete ich Ihnen 25% des Betrages an und werde die verbleibenden 70% fuer meine Investitionen in Ihrem Land nutzen. Sie können mich jederzeit unter dieser e-Mail Adresse erreichen: ( [EMAIL PROTECTED] ) Ich bitte Sie jedoch um Ihre absolute Verschwiegenheit in dieser Angelegenheit Dritten gegenüber. Gott schutze Sie! Hochachtungsvoll, Ihr ergebener, Ehret Olds. __ 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] Fetching Intraday data from Bloomberg
Hi Everyone. I am downloading intraday Bloomberg data from R. The code I give is: library(zoo) library(chron) library(RBloomberg) conn-blpConnect(show.days=trading,na.action=previous.days,periodici ty=daily) dat-blpGetData(conn, VG1 Index, c(LAST_PRICE), start=as.chron(as.Date(2006-9-01, %Y-%m-%d)),barfields=OPEN, barsize=10,retval=data.frame) blpDisconnect(conn) write.table(dat,Z:\h1.csv,sep=,,row.names=F,col.names=T) Now, though I give show.days=trading, in the above code, my intraday downloads also have Saturdays and Sunday data (non-active days). Why is this so? Is it not updated in R? Note that this happens only for intraday downloads and not for the usual historical downloads where time is not considered. Please answer me ASAP... Thanks in advance. [[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] Confidence intervals for relative risk
Wolfgang, It is common to handle relative risk problems using Poisson regression. In your example you have 8 events out of 508 tries, and 0/500 in the second data set. tdata - data.frame(y=c(8,0), n=c(508,500), group=1:0) fit - glm(y ~ group + offset(log(n)), data=tdata, family=poisson) Because of the zero, the standard beta/se(beta) confidence intervals don't work. In fact, the actual MLE for the ratio is infinity, which the above glm decides is exp(33.85) = 5* 10^14 --- which is close enough to infinity for me. What you want is the lower limit of the interval, which you can find with a profile likelihood. That is, look at the curve of deviance vs beta, draw a horizontal line 3.84 units down from the top (chisq on 1 df), and where it itersects the curve is your confidence limit. For the above, since it is 2 groups with 2 coefficients, the deviance of the full model happens to be 0. Your approximation gave a lower limit of .97 which is about exp(0), so I'll guess a solution between -2 and 2. xx - seq(-2, 2, length=21) for (i in 1:21) { fit - glm(y ~ offset(group*xx[i] + log(n)), poisson, tdata) print(c(xx[i], fit$deviance)) } [1] -2.0 33.80736 [1] -1.8 31.02994 [1] -1.6 28.33138 [1] -1.4 25.72327 [1] -1.2 23.21769 [1] -1.0 20.82691 [1] -0.8 18.56281 [1] -0.6 16.43629 [1] -0.4 14.45668 [1] -0.2 12.63108 [1] 0.0 10.96387 [1] 0.2 9.45639 [1] 0.40 8.106805 [1] 0.60 6.910274 [1] 0.80 5.859303 [1] 1.00 4.944278 [1] 1.20 4.154088 [1] 1.40 3.476757 [1] 1.6 2.90003 [1] 1.8 2.41186 [1] 2.00 2.000785 So the exact lower limit is somewhere between exp(1.4) and exp(1.2), which is around 3.6. One can easily refine this by tucking the fit into an iterative search, or just resetting xx to a new range. The offset(log(n)) part of the model, BTW, is a well known trick in poisson models. It has to do with the fact that the likelihood is in terms of the number of events, but we want coefficients in terms of rates, and E(y) = rate*n. Adding an offset of xx[i] * group fits a model with the group effect fixed at xx, but allowing the intercept to vary. For more conceptual depth, you could look up 'rate regression' or 'standardized mortality ratio' in the Encyclopedia of Biostatistics --- it is in the latter computation that this approach is quite common. The idea of adding a fraction is found in Anscombe(1949), Transformations of Poisson, binomial and negative-binomial data. Biometrika, vol 35, p246-254, but for each single estimate not for the ratio. Terry Therneau Mayo Clinic __ 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] segfault in RMySQL dbConnect error handling
Hi there. I see in a post from 2002 that you got the following problem with RMySQL: con - dbConnect(m) Process R segmentation fault at Wed Aug 28 08:21:11 2002 I have the same problem today: drv=dbDriver(MySQL) dbConnect(drv) # or with pretty much any other failing options Program received signal SIGSEGV, Segmentation fault. 0x2b28d60c9fa0 in strlen () from /lib/libc.so.6 (gdb) bt #0 0x2b28d60c9fa0 in strlen () from /lib/libc.so.6 #1 0x2b28d609c7f4 in vfprintf () from /lib/libc.so.6 #2 0x2b28d60b6a69 in vsprintf () from /lib/libc.so.6 #3 0x2b28d60a1f76 in sprintf () from /lib/libc.so.6 #4 0x2b28d7f66b1b in ?? () #5 0x00fbda01 in ?? () #6 0xa285 in ?? () #7 0x011eda68 in ?? () #8 0x00d982a0 in ?? () #9 0x in ?? () (gdb) info regi rax0x8d 141 rbx0x1 1 rcx0x3 3 rdx0x7f9ccb08 140737481853704 rsi0x2b28d7f67f01 47454421942017 (argument 2) nrecognised R/S type for group rdi0x8d 141 (argument 1) rbp0x7f9cc990 0x7f9cc990 rsp0x7f9cc2f8 0x7f9cc2f8 (gdb) c *** caught segfault *** address 0x8d, cause 'memory not mapped' Traceback: 1: .Call(RS_MySQL_newConnection, drvId, con.params, groups, default.file, PACKAGE = .MySQLPkgName) 2: mysqlNewConnection(drv, ...) 3: .class1(object) 4: .class1(object) 5: is(object, Cl) 6: .valueClassTest(standardGeneric(dbConnect), DBIConnection, dbConnect) 7: dbConnect(drv) Possible actions: 1: abort (with core dump) 2: normal R exit 3: exit R without saving workspace 4: exit R saving workspace Selection: It appears that it is crashing while handling an error relating to the database. If I give it parameters which allow it to connect, it does not segfault. System information: XEON x86_64 with VT Ubuntu 6.06.1 LTS \n \lRN) MySQL 5.0.22-Debian_0ubuntu6.06.2-log R 2.4.0 R version 2.4.0 (2006-10-03) from deb http://cran.R- project.org/bin/linux/ubuntu dapper/ It's a pristine 2.4.0 R install with the following libraries: DBI_0.1-11.tar.gz RMySQL_0.5-10.tar.gz installed in a custom library location /export/downloads/R/packages (normal location seems to behave the same) note: it also failed on 2.3.0 on this machine --Alex Brown __ 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] random forest regression
Dear all, I am doing a regression in ramdomForest, using the option sampsize reduce the number of records used to produce the randomForest object. The manual says For classification, if sampsize is a vector of the length the number of strata, then sampling is stratified by strata, and the elements of sampsize indicate the numbers to be drawn from the strata. I need my sampling to be done with factors, but I am doing a regression. Does anyone know a way to do that? Thanks, Naiara. [[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] R and Fortran 9x -- advice
Tamas K Papp [EMAIL PROTECTED] wrote: I found some bottlenecks in my R code with Rprof. First I wanted to rewrite them in C, but a colleague keeps suggesting that I learn Fortran, so maybe this is the time to do it... 1) I hear bad things about Fortran. Sure, F77 looks archaic, but F90/95 seems nicer. Is it worth learning, especially when I plan to use it mainly from R? Dusting off my C knowledge would take a bit of time too, and I never liked C that much. I'll answer this from the perspective of someone who uses Fortran 95 regularly. It is a modern language, far more reliable and flexible than Fortran 77 and quite well suited to most scientific problems. I do think it's worth learning, particularly if C is not to your taste. Two free compilers for Fortran 95 are available. It seems that g95 is complete, while gfortran is nearing completion. There are also several high-quality commercial compilers, some of which are free under certain operating systems and/or conditions and others of which (Lahey) are quite inexpensive if one is willing to work from the command line or one's own editor. I can't address questions of R interoperability -- not something I've done. -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ 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] Can I execute the content of a character vector?
Hi I want to know if there is any possibility of executing the content of a vector, for example: example=c(Test,1,0,0,0,seq(14,42,by=2),0,0,1) i want to know if there is anything like execute(example[6]) i really need this because this object example is created from a parameter file with names, flags and this seq somethimes will be something like c(1,99,3) or even c() executing the content would be an easier step for me thanks! -- Luiz Rodrigo Lins Tozzi [EMAIL PROTECTED] (21)91318150 http://luizrodrigotozzi.multiply.com/ http://www.flogao.com.br/luizrodrigotozzi [[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] [R-SIG-Finance] looking for functions that can test/estimate CAPM, APT, Fama's factor model, etc.
Hi Brian, thanks a lot for your pointers! I've taken a look at the book and the R example website. That's super! Some of the examples there are very good. Yet I am still looking for Fama 3 factor model and Ross' APT implementation. The concept is not hard per se, however I am not sure how to classify some companies as H, and some companies as L and others as Value companies, and the others as Growth companies. There are a lot of implementation details that I am not sure of. Thus I guess learning from other people's implementation is a better approach for me as a green-hand. (The Rmetrics are not very complete I guess.) Thanks a lot Brian and thanks everybody... On 11/13/06, Brian G. Peterson [EMAIL PROTECTED] wrote: On Sunday 12 November 2006 22:41, Michael wrote: I am also looking for interesting statistical experiments about testing and estimating CAPM, APT, Fama models, etc. using R using financial series data... please give me some pointers... I have been searching the R archives for the past a few hours and I vaguely got to know that there are programs do these interesting statistical things, but I just could not find where are they... Look at at the 'portfolio' and 'RMetrics' packages. RMetrics is actually a group of packages, not just a single module. These should give you CAPM and more. I can also recommend this excellent introductory book: http://www.amazon.com/Statistics-Finance-Introduction-David-Ruppert/dp/0387202706 Statistics and Finance: An Introduction by David Ruppert with R examples available here: http://www.stat.tamu.edu/~ljin/Finance/stat689-R.htm Once you've worked on this for a while, if you need pointers on a specific problem, this list may be able to be of greater assistance. Regards, - Brian ___ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-finance [[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] Problem with fitted.value function
Who can help me in finding the right fitted-value-function?Based upon a linear model (2-dimensional, 1 result-variable and 1 input-variable). I found the perfect regression model (not always simply linear) by the function lm(). Based on the given input and output tuples. For these exact input values I can get the exact fitted values by the function fitted.values. The problem is, once I want to extrapolate with the automaticly found function (not always the same) and new input-variables I don't know how to calculate automaticly the fitted values. I don't know how to use the fitted values function with a given function and given input-variables but yet unknown result-values. Who can help me get the right fitted value function and tell me how to use it. Thanks in advance for any help. __ 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] Can I execute the content of a character vector?
On Mon, 13 Nov 2006, Luiz Rodrigo Tozzi wrote: Hi I want to know if there is any possibility of executing the content of a vector, for example: example=c(Test,1,0,0,0,seq(14,42,by=2),0,0,1) i want to know if there is anything like execute(example[6]) You can say: eval(parse(text=example[6])) but it is difficult to avoid shooting people in the feet with it, unless that is what you want ... eval(parse(text=example[1])) i really need this because this object example is created from a parameter file with names, flags and this seq somethimes will be something like c(1,99,3) or even c() executing the content would be an easier step for me thanks! -- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: [EMAIL PROTECTED] __ 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] Can I execute the content of a character vector?
On 11/13/06, Roger Bivand [EMAIL PROTECTED] wrote: On Mon, 13 Nov 2006, Luiz Rodrigo Tozzi wrote: Hi I want to know if there is any possibility of executing the content of a vector, for example: example=c(Test,1,0,0,0,seq(14,42,by=2),0,0,1) i want to know if there is anything like execute(example[6]) You can say: eval(parse(text=example[6])) but it is difficult to avoid shooting people in the feet with it, unless that is what you want ... eval(parse(text=example[1])) We could enclose strings to be evaluated in backticks: example - c(Test,1,0,0,0,`seq(14,42,by=2)`,0,0,1) pat - ^`(.*)`$ f - function(x) # strip backticks in 1st and last posn and evaluate if (regexpr(pat, x) 0) eval.parent(parse(text = sub(pat, \\1, x))) else x example - c(Test,1,0,0,0,`seq(14,42,by=2)`,0,0,1) lapply(example, f) i really need this because this object example is created from a parameter file with names, flags and this seq somethimes will be something like c(1,99,3) or even c() executing the content would be an easier step for me thanks! -- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: [EMAIL PROTECTED] __ 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] R and Fortran 9x -- advice
Mike, Can you recommend any good books on Fortran 90/95? I had been an old user of Fortran 77 but haven't followed the developments in the last 15 years or so... Thanks. Christos Hatzis, Ph.D. Nuvera Biosciences, Inc. 400 West Cummings Park Suite 5350 Woburn, MA 01801 Tel: 781-938-3830 www.nuverabio.com -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Mike Prager Sent: Monday, November 13, 2006 1:00 PM To: r-help@stat.math.ethz.ch Subject: Re: [R] R and Fortran 9x -- advice Tamas K Papp [EMAIL PROTECTED] wrote: I found some bottlenecks in my R code with Rprof. First I wanted to rewrite them in C, but a colleague keeps suggesting that I learn Fortran, so maybe this is the time to do it... 1) I hear bad things about Fortran. Sure, F77 looks archaic, but F90/95 seems nicer. Is it worth learning, especially when I plan to use it mainly from R? Dusting off my C knowledge would take a bit of time too, and I never liked C that much. I'll answer this from the perspective of someone who uses Fortran 95 regularly. It is a modern language, far more reliable and flexible than Fortran 77 and quite well suited to most scientific problems. I do think it's worth learning, particularly if C is not to your taste. Two free compilers for Fortran 95 are available. It seems that g95 is complete, while gfortran is nearing completion. There are also several high-quality commercial compilers, some of which are free under certain operating systems and/or conditions and others of which (Lahey) are quite inexpensive if one is willing to work from the command line or one's own editor. I can't address questions of R interoperability -- not something I've done. -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ 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] R and Fortran 9x -- advice
Metcalf and Reid - FORTRAN 90/95 Explained Ravi. --- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Christos Hatzis Sent: Monday, November 13, 2006 2:05 PM To: 'Mike Prager'; r-help@stat.math.ethz.ch Subject: Re: [R] R and Fortran 9x -- advice Mike, Can you recommend any good books on Fortran 90/95? I had been an old user of Fortran 77 but haven't followed the developments in the last 15 years or so... Thanks. Christos Hatzis, Ph.D. Nuvera Biosciences, Inc. 400 West Cummings Park Suite 5350 Woburn, MA 01801 Tel: 781-938-3830 www.nuverabio.com -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Mike Prager Sent: Monday, November 13, 2006 1:00 PM To: r-help@stat.math.ethz.ch Subject: Re: [R] R and Fortran 9x -- advice Tamas K Papp [EMAIL PROTECTED] wrote: I found some bottlenecks in my R code with Rprof. First I wanted to rewrite them in C, but a colleague keeps suggesting that I learn Fortran, so maybe this is the time to do it... 1) I hear bad things about Fortran. Sure, F77 looks archaic, but F90/95 seems nicer. Is it worth learning, especially when I plan to use it mainly from R? Dusting off my C knowledge would take a bit of time too, and I never liked C that much. I'll answer this from the perspective of someone who uses Fortran 95 regularly. It is a modern language, far more reliable and flexible than Fortran 77 and quite well suited to most scientific problems. I do think it's worth learning, particularly if C is not to your taste. Two free compilers for Fortran 95 are available. It seems that g95 is complete, while gfortran is nearing completion. There are also several high-quality commercial compilers, some of which are free under certain operating systems and/or conditions and others of which (Lahey) are quite inexpensive if one is willing to work from the command line or one's own editor. I can't address questions of R interoperability -- not something I've done. -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ 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-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] 95% CI of spectrum
Dear R-list Members, Sorry if I have double posted this message. I am not sure if my previous message got through. Can someone explain me what exactly is the vertical bar of the cross in the upper right corner of the spectrum plot produced by spec.plot? I know it is the 95% confidence of the spectrum, but do not understand what it is exactly. Based on the section 9.5 of Bloomfield 2000, I understand how the pointwise 95% CI of spectrum at a certain frequency is computed. Is the vertical bar of the cross for the mean spectrum or total spectrum? What is the exact interpretation of the vertical bar? Many thanks in advance! Sincerely, Guiming Wang __ 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 with syntax of nlme call.
I am getting an error message in a call to nlme and cannot understand what is happening. I explain the steps below in the hope that someone can explain the error and how to correct it. STEP 1: Data set: name: marouane.data. This is a data frame whose first few lines are as follows: marouane.data[1:13,] species plant leaf irradiance photosynthesis chlorophyll 1 Asclepias.incarnata 11 0 -2.091359 0.02619 2 Asclepias.incarnata 11 50 1.153241 0.02619 3 Asclepias.incarnata 11100 2.241963 0.02619 4 Asclepias.incarnata 11200 3.541258 0.02619 5 Asclepias.incarnata 11300 5.012547 0.02619 6 Asclepias.incarnata 11400 5.710689 0.02619 7 Asclepias.incarnata 11600 6.632571 0.02619 8 Asclepias.incarnata 11800 7.314284 0.02619 9 Asclepias.incarnata 11 1000 8.092402 0.02619 10 Asclepias.incarnata 11 1310 8.110368 0.02619 11 Asclepias.incarnata 12 0 -2.051934 0.02707 12 Asclepias.incarnata 12 50 1.213854 0.02707 13 Asclepias.incarnata 12100 2.271453 0.02707 absorbance nitrogen sla 1 0.8816 18.35913 77.53 2 0.8816 18.35913 77.53 3 0.8816 18.35913 77.53 4 0.8816 18.35913 77.53 5 0.8816 18.35913 77.53 6 0.8816 18.35913 77.53 7 0.8816 18.35913 77.53 8 0.8816 18.35913 77.53 9 0.8816 18.35913 77.53 10 0.8816 18.35913 77.53 11 0.8813 18.35913 77.53 12 0.8813 18.35913 77.53 13 0.8813 18.35913 77.53 The nesting structure is species (25), plants within species (4 each species) and leaves within plants (2 each plant). There is only 1 missing value in the data set. STEP 2: I constructed a self starting function (photo) that gives the nonlinear function and provides initial estimates of the three parameters. This self starting function works. I use this to call the nlsList function, which fits separate nonlinear functions to each species (I am only working on a 2 level model so far: between and within species): fit-nlsList(photosynthesis~photo(irradiance,Am,Q,LCP)|species, + data=marouane.data,na.action=na.omit) This works and I get the three parameter estimates per species. STEP 3: use nlsList to call nlme. I use the nlsList object (fit) to fit a variance components model: fit.nlme-nlme(fit) This works (at least it runs without an error message). To see what syntax is used, I extract the call: fit.nlme$call nlme.formula(model = photosynthesis ~ photo(irradiance, Am, Q, LCP), data = marouane.data, fixed = list(Am ~ 1, Q ~ 1, LCP ~ 1), random = list(species = c(2.41887429868478, -1.53351210977838, 2.47282942435836, 0, 0, 0)), groups = ~species, start = list( fixed = c(11.7384877502532, 0.109091641284836, 9.91905527125462 )), na.action = na.omit) Question 1: the syntax to the random argument seems wrong! It should be something like: random=(list(Am~1,Q~1,LCP~1)). I have no idea where the 6 numbers (2.41887...) are coming from in the random argument. Can someone explain how the nlme function obtains such an argument for random? The values in fixed are the estimated fixed effects from nlsList. If I follow the instructions in Pinheiro Bates, the call should be (with a self-starting function): nlme(model = photosynthesis ~ photo(irradiance, Am, Q, LCP), data = marouane.data, fixed = list(Am ~ 1, Q ~ 1, LCP ~ 1), random = list(Am~1,Q~1,LCP~1),groups=~species, na.action=na.omit). When I use this, I get an error message: nlme(model=photosynthesis~photo(irradiance,Am,Q,LCP),data=marouane.data, + fixed=list(Am~1,Q~1,LCP~1),groups=~species, + random=list(Am~1,Q~1,LCP~1), + na.action = na.omit) Error in nlmeCall[[i]] - NULL : subscript out of bounds Bill Shipley [[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] indexing question
I have the following set of indices, call it idx, that correspond to the indices of a vector say temp. [1] 31 36 41 61 66 71 91 96 101 121 126 131 151 156 161 181 186 191 211 216 221 241 246 251 271 276 281 301 306 311 331 336 341 361 366 [36] 371 391 396 401 421 426 431 451 456 461 481 486 491 511 516 521 541 546 551 571 576 581 601 606 611 631 636 641 661 666 671 691 696 701 721 [71] 726 731 751 756 761 781 786 791 811 816 821 841 846 851 871 876 881 901 906 911 931 936 941 961 966 971 991 996 1001 1021 1026 1031 1051 1056 1061 [106] 1081 1086 1091 1116 1121 1141 1146 1151 1171 1176 1181 1201 1206 1211 1231 1236 1241 1261 1266 1271 1291 1296 1301 1321 1326 1331 1351 1356 1361 1381 1386 1391 1411 1416 [141] 1421 I want to calculate temp[36] - temp[31] and temp[41] - temp[36] Similarly, temp[66] - temp[61] and temp[71] - temp[66] . . . . Similarly temp[1416]-temp[1411] temp[1421] - temp[1416] I'm doing this because the above subractions represent pairs of returns and the correlations between them wil be calculated eventually. In other words, eventually I will have X_36_31 ( i.e : temp[36] - temp[31] ) X_66-61 X_96-91 . . . . . . . X_1411-1416 as one vector and Y_41-36 Y_71-66 Y_101-96 . . . . . Y_1416_1421 as another vector. and will calculate the correlation between the two vectors in order to get one number. The point is I am really only using the indices 31, 61, 91 etc as anchor's so a regular diff(temp[idx]) won't work because it will diff all the elements that are next to each other ? This is a weird problem. I'm still thinking about it. I'm hoping to figure it out before someone sends me something but I won't mind so much if I get an external solution first. I have no pride. This is not an offer (or solicitation of an offer) to buy/se...{{dropped}} __ 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] hybrid in fisher.test broken?
The hybrid feature in fisher.test looks to me like an excellent way to analyze my two-way tables. The only problem is that it does not seem to be implemented. Am I right about this? An example is pasted below. I note that I get the warning message only when I shouldn't: for a 2x2 table hybrid seems to be ignored without warning. In no case does fisher.test seem to be checking Cochran criteria as promised. Any help will be appreciated. JD -- R : Copyright 2006, The R Foundation for Statistical Computing Version 2.3.1 (2006-06-01) . . . m = matrix(10*1:9, nc=3) fisher.test(m, hybrid=TRUE) Error in fisher.test(m, hybrid = TRUE) : FEXACT error 6. LDKEY is too small for this problem. Try increasing the size of the workspace. In addition: Warning message: 'hybrid' is ignored for a 2 x 2 table in: fisher.test(m, hybrid = TRUE) __ 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] indexing question
diff(tmp[idx]) cheers, b On Nov 13, 2006, at 3:06 PM, Leeds, Mark ((IED)) wrote: I have the following set of indices, call it idx, that correspond to the indices of a vector say temp. [1] 31 36 41 61 66 71 91 96 101 121 126 131 151 156 161 181 186 191 211 216 221 241 246 251 271 276 281 301 306 311 331 336 341 361 366 [36] 371 391 396 401 421 426 431 451 456 461 481 486 491 511 516 521 541 546 551 571 576 581 601 606 611 631 636 641 661 666 671 691 696 701 721 [71] 726 731 751 756 761 781 786 791 811 816 821 841 846 851 871 876 881 901 906 911 931 936 941 961 966 971 991 996 1001 1021 1026 1031 1051 1056 1061 [106] 1081 1086 1091 1116 1121 1141 1146 1151 1171 1176 1181 1201 1206 1211 1231 1236 1241 1261 1266 1271 1291 1296 1301 1321 1326 1331 1351 1356 1361 1381 1386 1391 1411 1416 [141] 1421 I want to calculate temp[36] - temp[31] and temp[41] - temp[36] Similarly, temp[66] - temp[61] and temp[71] - temp[66] . . . . Similarly temp[1416]-temp[1411] temp[1421] - temp[1416] I'm doing this because the above subractions represent pairs of returns and the correlations between them wil be calculated eventually. In other words, eventually I will have X_36_31 ( i.e : temp[36] - temp[31] ) X_66-61 X_96-91 . . . . . . . X_1411-1416 as one vector and Y_41-36 Y_71-66 Y_101-96 . . . . . Y_1416_1421 as another vector. and will calculate the correlation between the two vectors in order to get one number. The point is I am really only using the indices 31, 61, 91 etc as anchor's so a regular diff(temp[idx]) won't work because it will diff all the elements that are next to each other ? This is a weird problem. I'm still thinking about it. I'm hoping to figure it out before someone sends me something but I won't mind so much if I get an external solution first. I have no pride. This is not an offer (or solicitation of an offer) to buy/se... {{dropped}} __ 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] indexing question
thanks beilton but that won't work. A diff will also include 61-41 etc and I don't want to include those. I'm working on using lapply or sapply with a seq along 31, 61, etc. I'll let you know if it works. -Original Message- From: Benilton Carvalho [mailto:[EMAIL PROTECTED] Sent: Monday, November 13, 2006 3:18 PM To: Leeds, Mark (IED) Cc: r-help@stat.math.ethz.ch Subject: Re: [R] indexing question diff(tmp[idx]) cheers, b On Nov 13, 2006, at 3:06 PM, Leeds, Mark ((IED)) wrote: I have the following set of indices, call it idx, that correspond to the indices of a vector say temp. [1] 31 36 41 61 66 71 91 96 101 121 126 131 151 156 161 181 186 191 211 216 221 241 246 251 271 276 281 301 306 311 331 336 341 361 366 [36] 371 391 396 401 421 426 431 451 456 461 481 486 491 511 516 521 541 546 551 571 576 581 601 606 611 631 636 641 661 666 671 691 696 701 721 [71] 726 731 751 756 761 781 786 791 811 816 821 841 846 851 871 876 881 901 906 911 931 936 941 961 966 971 991 996 1001 1021 1026 1031 1051 1056 1061 [106] 1081 1086 1091 1116 1121 1141 1146 1151 1171 1176 1181 1201 1206 1211 1231 1236 1241 1261 1266 1271 1291 1296 1301 1321 1326 1331 1351 1356 1361 1381 1386 1391 1411 1416 [141] 1421 I want to calculate temp[36] - temp[31] and temp[41] - temp[36] Similarly, temp[66] - temp[61] and temp[71] - temp[66] . . . . Similarly temp[1416]-temp[1411] temp[1421] - temp[1416] I'm doing this because the above subractions represent pairs of returns and the correlations between them wil be calculated eventually. In other words, eventually I will have X_36_31 ( i.e : temp[36] - temp[31] ) X_66-61 X_96-91 . . . . . . . X_1411-1416 as one vector and Y_41-36 Y_71-66 Y_101-96 . . . . . Y_1416_1421 as another vector. and will calculate the correlation between the two vectors in order to get one number. The point is I am really only using the indices 31, 61, 91 etc as anchor's so a regular diff(temp[idx]) won't work because it will diff all the elements that are next to each other ? This is a weird problem. I'm still thinking about it. I'm hoping to figure it out before someone sends me something but I won't mind so much if I get an external solution first. I have no pride. This is not an offer (or solicitation of an offer) to buy/se... {{dropped}} __ 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. This is not an offer (or solicitation of an offer) to buy/se...{{dropped}} __ 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] handling time units--hours, minutes, days--from file times
Dear R-helpers, I am trying to generalize my function for recording measurement times from file times mtime--my intervals are minutes to hours over the course of several days. I want to use hours as my units, and I have had trouble dealing with time units in a general way. I have a simple solution for the dealing with regular intervals, but I have not been able to handle irregular intervals. I tried using the zoo package, and had similar problems with the conversion of time units. Code and sample data is below-- for your info--after this, I am feeding the times and measurements to nls. Thanks in advance, -- Warren #Here is the failed code: myfilenames-row.names(myfileinfo) fileno-length(myfilenames) filetimes-numeric() for (i in 2:fileno){rel.read.time-myfileinfo$mtime[i]-myfileinfo$mtime[1] rel.read.time-(as.numeric(rel.read.time)) ##this next part is problematic ##i would love to say time.in.hours-as.hours (myfileinfo$mtime[i]-myfileinfo$mtime[1]) ifelse(rel.read.time =5, hour.read.time-(rel.read.time /60), hour.read.time-rel.read.time) filetimes-c(filetimes,rel.read.time)} #For regular intervals of readings, the code below generates measurement times from the first interval: myfilenames-row.names(myfileinfo) fileno-length(myfilenames) intervalue- myfileinfo$mtime[2] - myfileinfo$mtime[1] unitless.interval-(as.numeric(intervalue)) ifelse(unitless.interval=5, hour.interval-(unitless.interval/60), hour.interval-unitless.interval) hours-0 for (i in 1:(fileno-1)){hours-c(hours,(i*hour.interval))} myfileinfo - size isdir mode mtime ctime G2659310 2006-307-21-09-57.txt 10220 FALSE 666 2006-11-03 21:04:00 2006-11-13 10:56:39 G2659310 2006-307-21-45-55.txt 10230 FALSE 666 2006-11-03 21:40:00 2006-11-13 10:56:39 G2659310 2006-307-22-23-00.txt 10236 FALSE 666 2006-11-03 22:17:00 2006-11-13 10:56:39 G2659310 2006-307-23-00-33.txt 10236 FALSE 666 2006-11-03 22:54:00 2006-11-13 10:56:39 G2659310 2006-307-23-39-39.txt 10234 FALSE 666 2006-11-03 23:33:00 2006-11-13 10:56:39 G2659310 2006-308-00-17-25.txt 10234 FALSE 666 2006-11-04 00:11:00 2006-11-13 10:56:39 G2659310 2006-308-00-53-45.txt 10228 FALSE 666 2006-11-04 00:48:00 2006-11-13 10:56:39 G2659310 2006-308-01-30-03.txt 10208 FALSE 666 2006-11-04 01:25:00 2006-11-13 10:56:39 G2659310 2006-308-02-06-00.txt 10188 FALSE 666 2006-11-04 02:00:00 2006-11-13 10:56:39 G2659310 2006-308-02-41-57.txt 10168 FALSE 666 2006-11-04 02:36:00 2006-11-13 10:56:39 G2659310 2006-308-03-17-54.txt 10158 FALSE 666 2006-11-04 03:12:00 2006-11-13 10:56:39 G2659310 2006-308-03-53-51.txt 10148 FALSE 666 2006-11-04 03:48:00 2006-11-13 10:56:39 G2659310 2006-308-04-29-48.txt 10144 FALSE 666 2006-11-04 04:24:00 2006-11-13 10:56:39 G2659310 2006-308-05-05-44.txt 10132 FALSE 666 2006-11-04 05:00:00 2006-11-13 10:56:39 G2659310 2006-308-05-41-41.txt 10136 FALSE 666 2006-11-04 05:35:00 2006-11-13 10:56:39 G2659310 2006-308-06-17-39.txt 10128 FALSE 666 2006-11-04 06:11:00 2006-11-13 10:56:39 G2659310 2006-308-06-52-27.txt 10128 FALSE 666 2006-11-04 06:47:00 2006-11-13 10:56:39 G2659310 2006-308-07-28-25.txt 10130 FALSE 666 2006-11-04 07:23:00 2006-11-13 10:56:39 G2659310 2006-308-08-04-22.txt 10130 FALSE 666 2006-11-04 07:59:00 2006-11-13 10:56:39 G2659310 2006-308-08-40-20.txt 10134 FALSE 666 2006-11-04 08:35:00 2006-11-13 10:56:39 G2659310 2006-308-09-16-18.txt 10166 FALSE 666 2006-11-04 09:11:00 2006-11-13 10:56:39 G2659310 2006-308-09-52-15.txt 10180 FALSE 666 2006-11-04 09:46:00 2006-11-13 10:56:39 G2659310 2006-308-10-28-13.txt 10192 FALSE 666 2006-11-04 10:22:00 2006-11-13 10:56:39 G2659310 2006-308-11-04-10.txt 10210 FALSE 666 2006-11-04 10:58:00 2006-11-13 10:56:39 G2659310 2006-308-11-40-08.txt 10210 FALSE 666 2006-11-04 11:34:00 2006-11-13 10:56:39 G2659310 2006-308-12-16-06.txt 10222 FALSE 666 2006-11-04 12:10:00 2006-11-13 10:56:39 G2659310 2006-308-12-50-55.txt 10222 FALSE 666 2006-11-04 12:46:00 2006-11-13 10:56:39 G2659310 2006-308-13-26-54.txt 10224 FALSE 666 2006-11-04 13:21:00 2006-11-13 10:56:39 G2659310 2006-308-14-02-52.txt 10220 FALSE 666 2006-11-04 13:57:00 2006-11-13 10:56:39 G2659310 2006-308-14-38-51.txt 10224 FALSE 666 2006-11-04 14:33:00 2006-11-13 10:56:39 G2659310 2006-308-15-14-50.txt 10216 FALSE 666 2006-11-04 15:09:00 2006-11-13 10:56:39 G2659310 2006-308-15-50-47.txt 10202 FALSE 666 2006-11-04 15:45:00 2006-11-13 10:56:39 G2659310 2006-308-16-26-45.txt 10202 FALSE 666 2006-11-04 16:21:00 2006-11-13 10:56:39 G2659310 2006-308-17-02-43.txt 10190 FALSE 666 2006-11-04 16:57:00 2006-11-13 10:56:39 G2659310 2006-308-17-38-41.txt 10188 FALSE 666 2006-11-04 17:32:00 2006-11-13 10:56:39 G2659310 2006-308-18-14-40.txt 10174 FALSE 666 2006-11-04 18:08:00 2006-11-13 10:56:39 G2659310 2006-308-18-49-29.txt 10170 FALSE 666 2006-11-04 18:44:00 2006-11-13 10:56:39 G2659310 2006-308-19-25-28.txt 10158 FALSE 666 2006-11-04 19:20:00 2006-11-13 10:56:39 G2659310
Re: [R] indexing question
On Tuesday 14 November 2006 09:28, Leeds, Mark (IED) wrote: thanks beilton but that won't work. A diff will also include 61-41 etc and I don't want to include those. I'm working on using lapply or sapply with a seq along 31, 61, etc. I'll let you know if it works. Try looking at: dim(idx) - c(3, length(idx)/3) then: tmp[idx[2, ]] - temp[idx[1, ]] temp[idx[3, ]] - temp[idx[2, ]] HTH Ray Brownrigg __ 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] R and Fortran 9x -- advice
Christos Hatzis [EMAIL PROTECTED] wrote: Can you recommend any good books on Fortran 90/95? I had been an old user of Fortran 77 but haven't followed the developments in the last 15 years or so... Ravi Varadhan has already recommended Metcalf and Reid - Fortran 90/95 Explained. I own and like that book, so I'll second that, with a few caveats-- 1. The same authors + M. Cohen have released Fortran 95/2003 Explained. If you buy, why not buy the latest? There are as yet no Fortran 2003 compilers, but F95 compilers are starting to introduce F2003 features. 2. The style of MR(C) is quite condensed. Think R help pages. If you respond to terse explanations, this is your book. If not, check Amazon.com. 3. The book is NOT specifically written to help those moving from F77 to F95. For that, you might look at a Web tutorial, such as http://www.nsc.liu.se/~boein/f77to90/f77to90.html --and a couple of comments: Because standard F95, unlike F77, meets most programming needs, many F95 programmers try to write strictly standard-conforming code, rather than use vendor extensions. (I have ported a 10,000+ line program from Windows to Linux with NO changes.) A book with a detailed explanation of the Fortran 95 standard is Fortran 95 Handbook by Adams et al. I have a copy and use it rarely, but when I just can't seem to find something in MRC, I'm glad to have it. The Usenet group comp.lang.fortran has a high S/N ratio and courteous tone, with several members of the Fortran standards committee being frequent contributors. It's a great resource. HTH, ...Mike -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ 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] Nominal Respose Model in R
Hi: I have been working in Item Response Theory, exactly, with Nominal Response Model (NRM). Exist in R a function for estimate parameter and ability from database for this Model?. Thank you, Xavier G. Ordóñez [[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] question on MSM warning message
Hello After (simple random cluster) resampling with replacement I ran MSM function and I'm getting the following warning message ,which I'm not sure why. I don't have any absorbing stage set up in my MSM model. I have a 4 stage uni-directional MSM model. The only thing I see might be a problem is that when the last stage (stage 4 in my data) gets repeated it seems to give me a warning message. although, it picks only one observation whereas there might be other such cases in the data (printing only the first case?) out1.msm2 -msm(stage~years,subject=patient,data=mydata2,qmatrix=new.twoway9.q,con trol=list(trace=1,report=1),covariates=~black+tport11,constraint=list(bl ack=c(1,1,1),tport11=c(1,1,1))) ... Warning message: Absorbing - absorbing transition at observation 30 in: msm.check.model(msmdata$fromstate, msmdata$tostate, msmdata$obs, My data at observation around 30 look like ('years' is the time variable) mydata2[27:35,] count patient black years stage tport11 pcount 11029 1 0 1 1 17 1 1029 1 0 1 1 17 2 1029 1194 1 17 2 1029 1194 1 17 1 1033 0 0 1 1 20 1 1033 0 0 1 1 20 2 1033 0184 1 20 2 1033 0184 1 20 I'm wondering if I can ignore this warning message or should I have not done resampling with replacement. Any input would be greatly appreciated. Thank you so much for your time, kelly - KyungAh (Kelly) Im Epidemiology Data Center Graduate School of Public Health Department of Epidemiology University of Pittsburgh PA 15261 PHONE:(412) 624- 4612 FAX: (412) 624- 3775 Email: [EMAIL PROTECTED] http://www.edc.gsph.pitt.edu __ 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] For MacBook, best way to do R?
I'm getting a MacBook and I'd like to stick with OS X rather than convert it to Linux just yet. However, my main concern is having decent performance. What's my best option: *use the existing binary for R? *compile R fresh under OS X? * install Linux and run R under that? Does anyone have any recent experience they can share? Thanks! -- I can answer any question. I don't know is an answer. I don't know yet is a better answer. __ 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] A printing macro
Thanks for these suggestions, Professor Ripley. It's interesting that the function parameters in R are not truly dummy as they can effect the result of a function. Murray Prof Brian Ripley wrote: ddisp - function(dvar) { yn - substitute(dvar) csqt - eval.parent(substitute(chisq.test(G,dvar), list(dvar=yn))) } There are other ways, such as forming the cross-classification table, setting its dimnames and passing that to chisq.test. On Mon, 13 Nov 2006, Murray Jorgensen wrote: I am exploring the result of clustering a large multivariate data set into a number of groups, represented, say, by a factor G. I wrote a function to see how categorical variables vary between groups: ddisp - function(dvar) { + csqt - chisq.test(G,dvar) + print(csqt$statistic) + print(csqt$observed) + print(round(csqt$expected)) + round(csqt$residuals) + } x - ceiling(4*runif(100)) G - gl(4,1,100) ddisp(x) X-squared 6.235645 dvar G1 2 3 4 1 10 5 5 5 2 6 9 5 5 3 8 6 5 6 4 7 4 4 10 dvar G 1 2 3 4 1 8 6 5 6 2 8 6 5 6 3 8 6 5 6 4 8 6 5 6 dvar G1 2 3 4 1 1 0 0 -1 2 -1 1 0 -1 3 0 0 0 0 4 0 -1 0 1 Warning message: Chi-squared approximation may be incorrect in: chisq.test(G, dvar) As I need to apply this function to a large number of variables x it would be helpful if the function printed x rather than the formal argument dvar. I have a vague idea that things like deparse() and substitute() will come into the solution but I have not yet come up with the right incantation. Any help appreciated! Murray Jorgensen -- Dr Murray Jorgensen http://www.stats.waikato.ac.nz/Staff/maj.html Department of Statistics, University of Waikato, Hamilton, New Zealand Email: [EMAIL PROTECTED]Fax 7 838 4155 Phone +64 7 838 4773 wkHome +64 7 825 0441Mobile 021 1395 862 -- Dr Murray Jorgensen http://www.stats.waikato.ac.nz/Staff/maj.html Department of Statistics, University of Waikato, Hamilton, New Zealand Email: [EMAIL PROTECTED]Fax 7 838 4155 Phone +64 7 838 4773 wkHome +64 7 825 0441Mobile 021 1395 862 __ 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] Converting Character Strings to Numerical Values
Does R have anything like strtod, atoi, atod, etc. in C? I would like to covert 3 to 3 and 3.5 to 3.5 . Thanks, Peter Lauren. __ 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] Nominal Respose Model in R
AFAIK there exist functions to fit IRT models only for ordinal data. In particular, look at function grm() that fits the Graded Response Model in the 'ltm' package (using Marginal Maximum Likelihood), function PCM() in the 'eRm' package (using Condtional Maximum Likelihood), and function MCMCordfactanal() in the 'MCMCpack' package (using a Bayesian approach). I hope it helps. Best, Dimitris Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/(0)16/336899 Fax: +32/(0)16/337015 Web: http://med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm Quoting Xavier Giovanni Ordóñez Camacho [EMAIL PROTECTED]: Hi: I have been working in Item Response Theory, exactly, with Nominal Response Model (NRM). Exist in R a function for estimate parameter and ability from database for this Model?. Thank you, Xavier G. Ordóñez [[alternative HTML version deleted]] Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm __ 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] For MacBook, best way to do R?
Running R using the precomplied binary works really well and you should have no issues. Unless you really know what you are doing compiling R to work on a mac from sctrach can be at the very least tedious, at most, daunting, and there is really no point in doing so, since the binary already takes care of all the issues you may run into. In addition, you will need in install additional software (free though, from the OS X dvd) before you can compile. Lanre On 11/13/06, Mitchell Maltenfort [EMAIL PROTECTED] wrote: I'm getting a MacBook and I'd like to stick with OS X rather than convert it to Linux just yet. However, my main concern is having decent performance. What's my best option: *use the existing binary for R? *compile R fresh under OS X? * install Linux and run R under that? Does anyone have any recent experience they can share? Thanks! -- I can answer any question. I don't know is an answer. I don't know yet is a better answer. __ 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] Converting Character Strings to Numerical Values
?as.numeric x - 3 y - 3.5 x [1] 3 as.numeric(x) [1] 3 y [1] 3.5 as.numeric(y) [1] 3.5 Sarah On 11/13/06, Peter Lauren [EMAIL PROTECTED] wrote: Does R have anything like strtod, atoi, atod, etc. in C? I would like to covert 3 to 3 and 3.5 to 3.5 . Thanks, Peter Lauren. __ 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. -- Sarah Goslee http://www.functionaldiversity.org __ 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] Multivariate time-series
Hi all, I'm looking for R packages that estimate multivariate time-series models or vector-autoregression (VAR) time-series models. Thanks David -- === David Kaplan, Ph.D. Professor Department of Educational Psychology University of Wisconsin - Madison Educational Sciences, Room, 1061 1025 W. Johnson Street Madison, WI 53706 email: [EMAIL PROTECTED] homepage: http://www.education.wisc.edu/edpsych/facstaff/kaplan/kaplan.htm Phone: 608-262-0836 __ 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] Multivariate time-series
I think I remember seeing something called sem ? It's listed as on the packages on www.r-project.org and the explanation is what you want. I just can't be sure of the name. There are probably more than just sem. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of David Kaplan Sent: Monday, November 13, 2006 4:42 PM To: r-help@stat.math.ethz.ch Subject: [R] Multivariate time-series Hi all, I'm looking for R packages that estimate multivariate time-series models or vector-autoregression (VAR) time-series models. Thanks David -- === David Kaplan, Ph.D. Professor Department of Educational Psychology University of Wisconsin - Madison Educational Sciences, Room, 1061 1025 W. Johnson Street Madison, WI 53706 email: [EMAIL PROTECTED] homepage: http://www.education.wisc.edu/edpsych/facstaff/kaplan/kaplan.htm Phone: 608-262-0836 __ 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. This is not an offer (or solicitation of an offer) to buy/se...{{dropped}} __ 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] hybrid in fisher.test broken?
Jonathan Dushoff [EMAIL PROTECTED] writes: The hybrid feature in fisher.test looks to me like an excellent way to analyze my two-way tables. The only problem is that it does not seem to be implemented. Am I right about this? An example is pasted below. I note that I get the warning message only when I shouldn't: for a 2x2 table hybrid seems to be ignored without warning. In no case does fisher.test seem to be checking Cochran criteria as promised. Any help will be appreciated. Well, something seems to have been mangled in connection with these changes: r36358 | ripley | 2005-11-15 17:50:12 +0100 (Tue, 15 Nov 2005) | 2 lines add simulation option to fisher.test r36266 | ripley | 2005-11-10 22:19:34 +0100 (Thu, 10 Nov 2005) | 2 lines move 2x2 case in fisher.test to similar code as that used for or != 1 so that the warning else if (hybrid) { warning('hybrid' is ignored for a 2 x 2 table) now appears inside if (nr != 2 || nc != 2) { which is clearly wrong. However, the call to fexact that follows does have parameters set for the hybrid algorithm, so it's just the warning that is spurious. Notice that the hybrid algorithm may save a bit of time on probability calculations, but still needs to generate a large number of tables, so you really do need a larger workspace: fisher.test(matrix(10*1:9, nc=3),hybrid=TRUE,workspace=1e6) Fisher's Exact Test for Count Data data: matrix(10 * 1:9, nc = 3) p-value = 0.3149 alternative hypothesis: two.sided Warning message: 'hybrid' is ignored for a 2 x 2 table in: fisher.test(matrix(10 * 1:9, nc = 3), hybrid = TRUE, workspace = 1e+06) whereas (notice the slightly different p value) fisher.test(matrix(10*1:9, nc=3),hybrid=FALSE,workspace=1e6) Fisher's Exact Test for Count Data data: matrix(10 * 1:9, nc = 3) p-value = 0.3145 alternative hypothesis: two.sided JD -- R : Copyright 2006, The R Foundation for Statistical Computing Version 2.3.1 (2006-06-01) . . . m = matrix(10*1:9, nc=3) fisher.test(m, hybrid=TRUE) Error in fisher.test(m, hybrid = TRUE) : FEXACT error 6. LDKEY is too small for this problem. Try increasing the size of the workspace. In addition: Warning message: 'hybrid' is ignored for a 2 x 2 table in: fisher.test(m, hybrid = TRUE) __ 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. -- 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.
[R] Creating data for logistic regression and Cox proportional hazards regression
I know that mvrnorm from MASS (generously provided by Profs. Venables and Ripley) can be used to generate multivariable normal data that can be used in a linear regression with certain desired characteristics (e.g. a given mean for each variable as well as a given variance-covariance pattern). Is there any similar facility that can be used to generate data for (1) a logistic regression and (2) a Cox proportional hazards regression? Thanks, John John Sorkin M.D., Ph.D. Chief, Biostatistics and Informatics Baltimore VA Medical Center GRECC, University of Maryland School of Medicine Claude D. Pepper OAIC, University of Maryland Clinical Nutrition Research Unit, and Baltimore VA Center Stroke of Excellence University of Maryland School of Medicine Division of Gerontology Baltimore VA Medical Center 10 North Greene Street GRECC (BT/18/GR) Baltimore, MD 21201-1524 (Phone) 410-605-7119 (Fax) 410-605-7913 (Please call phone number above prior to faxing) [EMAIL PROTECTED] Confidentiality Statement: This email message, including any attachments, is for the so...{{dropped}} __ 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] Multivariate time-series
On Mon, 13 Nov 2006 15:42:06 -0600 David Kaplan wrote: Hi all, I'm looking for R packages that estimate multivariate time-series models or vector-autoregression (VAR) time-series models. The Econometrics task view at http://CRAN.R-project.org/src/contrib/Views/Econometrics.html has in the Time-series modelling section: For estimating VAR models, several methods are available: simple models can be fitted by ar() in stats, more elaborate models are provided package vars, estVARXls() in dse and a Bayesian approach is available in MSBVAR. hth, Z Thanks David -- === David Kaplan, Ph.D. Professor Department of Educational Psychology University of Wisconsin - Madison Educational Sciences, Room, 1061 1025 W. Johnson Street Madison, WI 53706 email: [EMAIL PROTECTED] homepage: http://www.education.wisc.edu/edpsych/facstaff/kaplan/kaplan.htm Phone: 608-262-0836 __ 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] Embedded carriage returns in text document
Colleagues, I am using R 2.4.0 on both a Mac (10.4.8) and Linux (RedHat 9). To read data from an Excel spreadsheet, I do save as in Excel, then select the Text (tab-delimited) format. The resulting file uses a tab separator and I can usually read the file using read.delim. Sometimes, the header row contains embedded carriage returns. When I view the file, these carriage returns appear as ^M. Now the problem: When I read.delim these files, they do not read correctly. Sometimes I get error messages; sometimes only the first line is read. Interestingly, invoking the option skip=1 (or a larger N) does not appear to bypass the problem. I can solve the problem by manually deleting these carriage returns either in the original Excel file or the .txt version. However, this is not an ideal solution. Does anyone have a work-around within R? Dennis Dennis Fisher MD P (The P Less Than Company) Phone: 1-866-PLessThan (1-866-753-7784) Fax: 1-415-564-2220 www.PLessThan.com [[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] Creating data for logistic regression and Cox proportional hazards regression
John Sorkin wrote: I know that mvrnorm from MASS (generously provided by Profs. Venables and Ripley) can be used to generate multivariable normal data that can be used in a linear regression with certain desired characteristics (e.g. a given mean for each variable as well as a given variance-covariance pattern). Is there any similar facility that can be used to generate data for (1) a logistic regression and (2) a Cox proportional hazards regression? Thanks, John See the examples in the help file for the lrm and cph functions in the Design package as well as in the spower function in the Hmisc package. Frank John Sorkin M.D., Ph.D. __ 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] hybrid in fisher.test broken?
Thank you very much for this useful response. I did a search for fexact both inside R and on R-project.org before posting: can somebody tell me where I could have found this information? In particular, I am still not clear on whether the Cochran criteria are being tested. Thanks, JD -- Jonathan Dushoff [EMAIL PROTECTED] writes: The hybrid feature in fisher.test looks to me like an excellent way to analyze my two-way tables. The only problem is that it does not seem to be implemented. Am I right about this? An example is pasted below. I note that I get the warning message only when I shouldn't: for a 2x2 table hybrid seems to be ignored without warning. In no case does fisher.test seem to be checking Cochran criteria as promised. Any help will be appreciated. Well, something seems to have been mangled in connection with these changes: r36358 | ripley | 2005-11-15 17:50:12 +0100 (Tue, 15 Nov 2005) | 2 lines add simulation option to fisher.test r36266 | ripley | 2005-11-10 22:19:34 +0100 (Thu, 10 Nov 2005) | 2 lines move 2x2 case in fisher.test to similar code as that used for or != 1 so that the warning else if (hybrid) { warning('hybrid' is ignored for a 2 x 2 table) now appears inside if (nr != 2 || nc != 2) { which is clearly wrong. However, the call to fexact that follows does have parameters set for the hybrid algorithm, so it's just the warning that is spurious. Notice that the hybrid algorithm may save a bit of time on probability calculations, but still needs to generate a large number of tables, so you really do need a larger workspace: fisher.test(matrix(10*1:9, nc=3),hybrid=TRUE,workspace=1e6) Fisher's Exact Test for Count Data data: matrix(10 * 1:9, nc = 3) p-value = 0.3149 alternative hypothesis: two.sided Warning message: 'hybrid' is ignored for a 2 x 2 table in: fisher.test(matrix(10 * 1:9, nc = 3), hybrid = TRUE, workspace = 1e+06) whereas (notice the slightly different p value) . . . __ 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] Fitting mean and covariance of Multivariate normal with censored data
Nobody answered, but this is what I did. I used an iterative poor man's imputation and filled in the data with a single value at each iteration. My covariance matrices and variance will be a little underestimated, but it'll do for me at this stage of this project. I iteratively filled in the data using the current data to estimate the parameters of the normal distribution. Using the current estimates, I figured out how much of the current distribution is outside each cut-off. Then fill-in at the mid-point of the probability density that is outside the area. I limited my corrections to 3*sigma. Here is the code section inside the loop (I am not showing another part of the code that does outlier elimination, nor the rest of the loop that recomputes the covariance and the parameters of the normal distribution) Trim the data to only use lines with at least 60% real data. Trim the columns to those that have enough data to estimate the parameters of the normal distribution (e.g. ~ 20) Loop over iter iterations ( from 1 to iterlast) Compute mean, variance, and co-variance matrix. Code to detect outliers.. Only activated after iterlast iterations Sort and trim worst outliers (left and right, up to 10% of data) For Each line of data, find out which elements were missing in original data set. for(im in missings) { # Compute the quantile on the cutoff away from N(mode,variance) for both cutoff quantOver-pnorm(c(highR),mean=mu[im],sd=sigmaMu[im]) quantUnder-pnorm(c(lowR),mean=mu[im],sd=sigmaMu[im]) # Find out which end of the distribution is the most outside of the limit quantOverMid- (1.0-quantOver)/2 quantUnderMid- quantUnder/2 # correct the mean, but underestimate the variance a little bit in cases # where width of pattern width of cutoff-region fixUp-mu[im] fixLow-mu[im] # Don't fix if distribution is cutoff at edges. if(quantOverMid0.05) { quantOverMid-0.0 } else { fixUp-qnorm(c(1.0-quantOverMid),mean=mu[im],sd=sigmaMu[im]) } if(quantUnderMid0.05) { quantUnderMid-0.0 } else { fixLow-qnorm(c(quantUnderMid),mean=mu[im],sd=sigmaMu[im]) } # Censored variables pull from both sides. quantNormSum - (quantOverMid+quantUnderMid) if(quantNormSum=0.05) { quantNormSum- 1.0/quantNormSum fix - (fixUp*quantOverMid+fixLow*quantUnderMid)*quantNormSum } else { # If missing data is in the tails, correct by cutoff. if(quantNormSum1e-10) { fix- (highR*(1.0-quantOver)+lowR*quantUnder)/((1.0-quantOver)+quantUnder) } else { fix-mu[im] } } if(!is.finite(fix)) { if(1.0-quantOver quantUnder) { fix-mu[im]+3*sigmaMu[im] } else { fix- mu[im]-3*sigmaMu[im] } } else { if(fixmu[im]+3*sigmaMu[im]) { fix-mu[im]+3*sigmaMu[im] } else if (fix mu[im]-3*sigmaMu[im]) { fix- mu[im]-3*sigmaMu[im] } } _ From: Sicotte, Hugues Ph.D. Sent: Wednesday, November 01, 2006 11:38 AM To: 'r-help@stat.math.ethz.ch' Subject: Fitting mean and covariance of Multivariate normal with censored data Hello, I have been googling for 2 days and I cannot find the answer in previous posts. I have a set of d-dimensional data elements (d=11 .. 14), each data point can be censored at different values both Lower-limit and upper limit. N = 2000 sets of vectors of D=11 data points per vector. Each of the N*D points can have different upper and lower limits. I simply want to fit a Multivariate distribution to the data and get the mean and covariance matrix of the sample. The closest post I saw mentionned the function survreg in the survival5 package, but I can't figure out how To use it for a
Re: [R] Embedded carriage returns in text document
Dennis Fisher [EMAIL PROTECTED] writes: Colleagues, I am using R 2.4.0 on both a Mac (10.4.8) and Linux (RedHat 9). To read data from an Excel spreadsheet, I do save as in Excel, then select the Text (tab-delimited) format. The resulting file uses a tab separator and I can usually read the file using read.delim. Sometimes, the header row contains embedded carriage returns. When I view the file, these carriage returns appear as ^M. Now the problem: When I read.delim these files, they do not read correctly. Sometimes I get error messages; sometimes only the first line is read. Interestingly, invoking the option skip=1 (or a larger N) does not appear to bypass the problem. I can solve the problem by manually deleting these carriage returns either in the original Excel file or the .txt version. However, this is not an ideal solution. Does anyone have a work-around within R? Hmm,... I suppose that CR messes with what R or the system thinks is the line-end character in this particular file. My first idea would be to read from a pipe() which executed sed 's/\r//' myfile.dat or something in that vein. Beware of quoting and differences in sed versions -- 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] handling time units--hours, minutes, days--from file times
hardly the most efficient way to go, but consider using a substring function to extract the time bits from your data, then reading them as POSIX dates and using difftime. for example, mytime - c(G2659310 2006-310-10-55-32.txt 10134 FALSE 666 2006-11-06 10:49:00 2006-11-13 10:56:41, G2659310 2006-310-11-31-43.txt 10136 FALSE 666 2006-11-06 11:25:00 2006-11-13 10:56:41, G2659310 2006-310-12-08-04.txt 10140 FALSE 666 2006-11-06 12:02:00 2006-11-13 10:56:41, G2659310 2006-310-12-43-17.txt 10148 FALSE 666 2006-11-06 12:38:00 2006-11-13 10:56:41, G2659310 2006-310-13-20-21.txt 10138 FALSE 666 2006-11-06 13:15:00 2006-11-13 10:56:41, G2659310 2006-310-13-56-28.txt 10080 FALSE 666 2006-11-06 13:51:00 2006-11-13 10:56:38) time1 - substring(mytime,49,67) time2 - substring(mytime,69,87) interval - difftime(as.POSIXct(time1, tz = AET10ADT), as.POSIXct(time2, tz = AET10ADT), units = hours) interval Time differences of -168.1281, -167.5281, -166.9114, -166.3114, -165.6947, -165.0939 hours i'm not exactly sure what you're after, but here's a vector of time differences. in general the POSIX class holds the most detailed information about time -- note the specification of tz (time zone) -- and should be used when a great deal of precision is required. tc On 11/14/06, Warren [EMAIL PROTECTED] wrote: Dear R-helpers, I am trying to generalize my function for recording measurement times from file times mtime--my intervals are minutes to hours over the course of several days. I want to use hours as my units, and I have had trouble dealing with time units in a general way. I have a simple solution for the dealing with regular intervals, but I have not been able to handle irregular intervals. I tried using the zoo package, and had similar problems with the conversion of time units. Code and sample data is below-- for your info--after this, I am feeding the times and measurements to nls. Thanks in advance, -- Warren #Here is the failed code: myfilenames-row.names(myfileinfo) fileno-length(myfilenames) filetimes-numeric() for (i in 2:fileno){rel.read.time-myfileinfo$mtime[i]-myfileinfo$mtime[1] rel.read.time-(as.numeric(rel.read.time)) ##this next part is problematic ##i would love to say time.in.hours-as.hours (myfileinfo$mtime[i]-myfileinfo$mtime[1]) ifelse(rel.read.time =5, hour.read.time-(rel.read.time /60), hour.read.time-rel.read.time) filetimes-c(filetimes,rel.read.time)} #For regular intervals of readings, the code below generates measurement times from the first interval: myfilenames-row.names(myfileinfo) fileno-length(myfilenames) intervalue- myfileinfo$mtime[2] - myfileinfo$mtime[1] unitless.interval-(as.numeric(intervalue)) ifelse(unitless.interval=5, hour.interval-(unitless.interval/60), hour.interval-unitless.interval) hours-0 for (i in 1:(fileno-1)){hours-c(hours,(i*hour.interval))} myfileinfo - size isdir mode mtime ctime G2659310 2006-307-21-09-57.txt 10220 FALSE 666 2006-11-03 21:04:00 2006-11-13 10:56:39 G2659310 2006-307-21-45-55.txt 10230 FALSE 666 2006-11-03 21:40:00 2006-11-13 10:56:39 G2659310 2006-307-22-23-00.txt 10236 FALSE 666 2006-11-03 22:17:00 2006-11-13 10:56:39 G2659310 2006-307-23-00-33.txt 10236 FALSE 666 2006-11-03 22:54:00 2006-11-13 10:56:39 G2659310 2006-307-23-39-39.txt 10234 FALSE 666 2006-11-03 23:33:00 2006-11-13 10:56:39 G2659310 2006-308-00-17-25.txt 10234 FALSE 666 2006-11-04 00:11:00 2006-11-13 10:56:39 G2659310 2006-308-00-53-45.txt 10228 FALSE 666 2006-11-04 00:48:00 2006-11-13 10:56:39 G2659310 2006-308-01-30-03.txt 10208 FALSE 666 2006-11-04 01:25:00 2006-11-13 10:56:39 G2659310 2006-308-02-06-00.txt 10188 FALSE 666 2006-11-04 02:00:00 2006-11-13 10:56:39 G2659310 2006-308-02-41-57.txt 10168 FALSE 666 2006-11-04 02:36:00 2006-11-13 10:56:39 G2659310 2006-308-03-17-54.txt 10158 FALSE 666 2006-11-04 03:12:00 2006-11-13 10:56:39 G2659310 2006-308-03-53-51.txt 10148 FALSE 666 2006-11-04 03:48:00 2006-11-13 10:56:39 G2659310 2006-308-04-29-48.txt 10144 FALSE 666 2006-11-04 04:24:00 2006-11-13 10:56:39 G2659310 2006-308-05-05-44.txt 10132 FALSE 666 2006-11-04 05:00:00 2006-11-13 10:56:39 G2659310 2006-308-05-41-41.txt 10136 FALSE 666 2006-11-04 05:35:00 2006-11-13 10:56:39 G2659310 2006-308-06-17-39.txt 10128 FALSE 666 2006-11-04 06:11:00 2006-11-13 10:56:39 G2659310 2006-308-06-52-27.txt 10128 FALSE 666 2006-11-04 06:47:00 2006-11-13 10:56:39 G2659310 2006-308-07-28-25.txt 10130 FALSE 666 2006-11-04 07:23:00 2006-11-13 10:56:39 G2659310 2006-308-08-04-22.txt 10130 FALSE 666 2006-11-04 07:59:00 2006-11-13 10:56:39 G2659310 2006-308-08-40-20.txt 10134 FALSE 666 2006-11-04 08:35:00 2006-11-13 10:56:39 G2659310 2006-308-09-16-18.txt 10166 FALSE 666 2006-11-04 09:11:00 2006-11-13 10:56:39 G2659310 2006-308-09-52-15.txt 10180 FALSE 666 2006-11-04 09:46:00 2006-11-13 10:56:39 G2659310
Re: [R] hybrid in fisher.test broken?
Jonathan Dushoff [EMAIL PROTECTED] writes: Thank you very much for this useful response. I did a search for fexact both inside R and on R-project.org before posting: can somebody tell me where I could have found this information? In particular, I am still not clear on whether the Cochran criteria are being tested. You need the actual source for FEXACT for this, then look at the 5th argument: EXPECT - Expected value used in the hybrid algorithm for deciding when to use asymptotic theory probabilities. (Input) If EXPECT = 0.0 then asymptotic theory probabilities are not used and Fisher exact test probabilities are computed. Otherwise, if PERCNT or more of the cells in the remaining table have estimated expected values of EXPECT or more, with no remaining cell having expected value less than EMIN, then asymptotic chi-squared probabilities are used. See the algorithm section of the manual document for details. Use EXPECT = 5.0 to obtain the 'Cochran' condition. and compare with PVAL - .C(fexact, nr, nc, x, nr, as.double(5), as.double(80), as.double(1), double(1), p = double(1), as.integer(workspace), mult = as.integer(mult), PACKAGE = stats)$p The source file in question is src/library/stats/src/fexact.c available in the R source tarballs or via https://svn.r-project.org/R Thanks, JD -- Jonathan Dushoff [EMAIL PROTECTED] writes: The hybrid feature in fisher.test looks to me like an excellent way to analyze my two-way tables. The only problem is that it does not seem to be implemented. Am I right about this? An example is pasted below. I note that I get the warning message only when I shouldn't: for a 2x2 table hybrid seems to be ignored without warning. In no case does fisher.test seem to be checking Cochran criteria as promised. Any help will be appreciated. Well, something seems to have been mangled in connection with these changes: r36358 | ripley | 2005-11-15 17:50:12 +0100 (Tue, 15 Nov 2005) | 2 lines add simulation option to fisher.test r36266 | ripley | 2005-11-10 22:19:34 +0100 (Thu, 10 Nov 2005) | 2 lines move 2x2 case in fisher.test to similar code as that used for or != 1 so that the warning else if (hybrid) { warning('hybrid' is ignored for a 2 x 2 table) now appears inside if (nr != 2 || nc != 2) { which is clearly wrong. However, the call to fexact that follows does have parameters set for the hybrid algorithm, so it's just the warning that is spurious. Notice that the hybrid algorithm may save a bit of time on probability calculations, but still needs to generate a large number of tables, so you really do need a larger workspace: fisher.test(matrix(10*1:9, nc=3),hybrid=TRUE,workspace=1e6) Fisher's Exact Test for Count Data data: matrix(10 * 1:9, nc = 3) p-value = 0.3149 alternative hypothesis: two.sided Warning message: 'hybrid' is ignored for a 2 x 2 table in: fisher.test(matrix(10 * 1:9, nc = 3), hybrid = TRUE, workspace = 1e+06) whereas (notice the slightly different p value) . . . -- 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.
[R] Forcing the intercept
Dear R-users: I am doing multiple regressions using the lm function and would like to force the intercept to be equal to a specific value (such as 4.3). I was able to find out how to force it through the origin but this does not work for other values. I am also interested in forcing the regression parameters obtained from one regression in another regression with a subset of the data. Are either of these possible in R? I have been searching the help guide for hours and have been unsuccessful. Many thanks, Heather __ 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] Question on applying vectors to data.frames by row
# Newbie alert # I am wanting to multiply the rows in a dataframe by a vector. # However, the default behavior appears to be for the vector to be applied # column wise. For example: vct - 1:4 df - data.frame(c1 = 5:10, c2= 6:11, c3=7:12, c4=8:13) multTheTwo - vct * df multTheTwo # This results in the vector getting cycled columnwise # c1 c2 c3 c4 # 1 5 18 7 24 # 2 12 28 16 36 # 3 21 8 27 10 # 4 32 18 40 22 # 5 9 30 11 36 # 6 20 44 24 52 # But what I actually want is: # c1 c2 c3 c4 #1 5 12 21 32which is 5*1, 6*2, 7*3, 8*4 #2 6 14 24 36same pattern applied to the next row #3 7 16 27 40so on #4 8 18 30 44 #5 9 20 33 48 #6 10 22 36 52 # I am currently using a for statement to do this, but that just doesn't feel # right. Is there another option that is more R like. for(i in 1:length(vct)) multTheTwo[,i] - vct[i] * df[,i] multTheTwo # Thanks Richard -- __ 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] Profile confidence intervals and LR chi-square test
System: R 2.3.1 on Windows XP machine. I am building a logistic regression model for a sample of 100 cases in dataframe d, in which there are 3 binary covariates: x1, x2 and x3. summary(d) y x1 x2 x3 0:54 0:50 0:64 0:78 1:46 1:50 1:36 1:22 fit - glm(y ~ x1 + x2 + x3, data=d, family=binomial(link=logit)) summary(fit) Call: glm(formula = y ~ x1 + x2 + x3, family = binomial(link = logit), data = d) Deviance Residuals: Min 1Q Median 3Q Max -1.6503 -1.0220 -0.7284 0.9965 1.7069 Coefficients: Estimate Std. Error z value Pr(|z|) (Intercept) -0.3772 0.3721 -1.014 0.3107 x11 -0.8144 0.4422 -1.842 0.0655 . x21 0.9226 0.4609 2.002 0.0453 * x31 1.3347 0.5576 2.394 0.0167 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 137.99 on 99 degrees of freedom Residual deviance: 120.65 on 96 degrees of freedom AIC: 128.65 Number of Fisher Scoring iterations: 4 exp(fit$coef) (Intercept) x11 x21 x31 0.6858006 0.4429233 2.5157321 3.7989873 --- After reading the appropriate sections in MASS4 (7.2 and 8.4 in particular), I decided to estimate the 95% confidence intervals for the odds ratios using the profile method implemented in the confint function. I then used the anova function to perform the deviance chi-square tests for each covariate. --- ci - confint(fit); exp(ci) Waiting for profiling to be done... 2.5 %97.5 % (Intercept) 0.3246680 1.413684 x11 0.1834819 1.048154 x21 1.0256096 6.314473 x31 1.3221533 12.129210 anova(fit, test='Chisq') Analysis of Deviance Table Model: binomial, link: logit Response: y Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev P(|Chi|) NULL99137.989 x115.85698132.133 0.016 x215.27197126.862 0.022 x316.21296120.650 0.013 My question relates to the interpretation of the significance of variable x1. The OR for x1 is 0.443 and its profile confidence interval is 0.183-1.048. If a type I error rate of 5% is assumed, this result would tend to suggest that x1 is NOT a significant predictor of y. However, the deviance chi-square test has a P-value of 0.016, which suggests that x1 is indeed a significant predictor of y. How do I reconcile these two differing messages? I do recognize that the upper bound of the confidence interval is pretty close to 1, but I am certain that some journal reviewer will point out the problem as inconsistent. Brant Inman __ 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] Forcing the intercept
you can just subtract 4.3 from the independent variable and then do through zero. That will Give you a force through 4.3. I don't undersarand the second part of your statement. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Heather Maughan Sent: Monday, November 13, 2006 6:31 PM To: r-help@stat.math.ethz.ch Subject: [R] Forcing the intercept Dear R-users: I am doing multiple regressions using the lm function and would like to force the intercept to be equal to a specific value (such as 4.3). I was able to find out how to force it through the origin but this does not work for other values. I am also interested in forcing the regression parameters obtained from one regression in another regression with a subset of the data. Are either of these possible in R? I have been searching the help guide for hours and have been unsuccessful. Many thanks, Heather __ 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. This is not an offer (or solicitation of an offer) to buy/se...{{dropped}} __ 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] Forcing the intercept
Second part: you want to look at the residuals after specifying a model? Just use the dependent variables in your subsample as the dependent variables in your regression equation and subtract that from your outcome variable in your subsample. Might not be the answer to the question you're asking though. Jeff. -- Jeffrey R. Spies http://www.nd.edu/~jspies/ On Nov 13, 2006, at 7:00 PM, Leeds, Mark (IED) wrote: you can just subtract 4.3 from the independent variable and then do through zero. That will Give you a force through 4.3. I don't undersarand the second part of your statement. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Heather Maughan Sent: Monday, November 13, 2006 6:31 PM To: r-help@stat.math.ethz.ch Subject: [R] Forcing the intercept Dear R-users: I am doing multiple regressions using the lm function and would like to force the intercept to be equal to a specific value (such as 4.3). I was able to find out how to force it through the origin but this does not work for other values. I am also interested in forcing the regression parameters obtained from one regression in another regression with a subset of the data. Are either of these possible in R? I have been searching the help guide for hours and have been unsuccessful. Many thanks, Heather __ 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. This is not an offer (or solicitation of an offer) to buy/se... {{dropped}} __ 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. [[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] Question on applying vectors to data.frames by row
Here are a couple of possibilities: as.data.frame(t(t(df) * vct)) df * rep(vct, each = nrow(df)) On 11/13/06, RDM [EMAIL PROTECTED] wrote: # Newbie alert # I am wanting to multiply the rows in a dataframe by a vector. # However, the default behavior appears to be for the vector to be applied # column wise. For example: vct - 1:4 df - data.frame(c1 = 5:10, c2= 6:11, c3=7:12, c4=8:13) multTheTwo - vct * df multTheTwo # This results in the vector getting cycled columnwise # c1 c2 c3 c4 # 1 5 18 7 24 # 2 12 28 16 36 # 3 21 8 27 10 # 4 32 18 40 22 # 5 9 30 11 36 # 6 20 44 24 52 # But what I actually want is: # c1 c2 c3 c4 #1 5 12 21 32which is 5*1, 6*2, 7*3, 8*4 #2 6 14 24 36same pattern applied to the next row #3 7 16 27 40so on #4 8 18 30 44 #5 9 20 33 48 #6 10 22 36 52 # I am currently using a for statement to do this, but that just doesn't feel # right. Is there another option that is more R like. for(i in 1:length(vct)) multTheTwo[,i] - vct[i] * df[,i] multTheTwo # Thanks Richard -- __ 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] Confidence interval for relative risk
Michael Dewey [EMAIL PROTECTED] wrote Subject: [R] Confidence interval for relative risk The concrete problem is that I am refereeing a paper where a confidence interval is presented for the risk ratio and I do not find it credible. I show below my attempts to do this in R. The example is slightly changed from the authors'. I can obtain a confidence interval for the odds ratio from fisher.test of course === fisher.test example === outcome - matrix(c(500, 0, 500, 8), ncol = 2, byrow = TRUE) fisher.test(outcome) Fisher's Exact Test for Count Data data: outcome p-value = 0.00761 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 1.694792 Inf sample estimates: odds ratio Inf === end example === Since RR = p1 +1 p1=a/(a+b), p2=c/(c+d) -- p2 1 - + --- 1-p2OR you can backsolve for the RR CI. x - function(p1, p2, oddsr) p1 + 1/(p2/(1-p2) + 1/oddsr) 1/x(8/508, 0, 1/1.694792) [1] 1.650734 Another approach is the Cornfield method -- which is a profile likelihood based CI using the Gibbs Chi-square (LRTS), but IIRC originally used the Pearson chi-square. I don't think there are R implementations. David Duffy. -- | David Duffy (MBBS PhD) ,-_|\ | email: [EMAIL PROTECTED] ph: INT+61+7+3362-0217 fax: -0101 / * | Epidemiology Unit, Queensland Institute of Medical Research \_,-._/ | 300 Herston Rd, Brisbane, Queensland 4029, Australia GPG 4D0B994A v __ 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] Generating confusion matrix / Kappa stats
Hi there, I classified an image and checked out it on the field. Now I have a table with three fields like: Field_ID Field_Class Image_Class 1 1 1 2 3 5 3 4 1 4 1 1 5 2 1 ... And now I need gerating a confusion matrix to compute Kappa statistic. First of all, how can I generate the confusion matrix from my input table? What package is good for compute Kappa statistics? Regards, Miltinho Brazil - [[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] Profile confidence intervals and LR chi-square test
On 2006-11-14 00:41, Inman, Brant A. M.D. skrev: System: R 2.3.1 on Windows XP machine. Time to upgrade! I am building a logistic regression model for a sample of 100 cases in dataframe d, in which there are 3 binary covariates: x1, x2 and x3. Please provide a reproducible example (as suggested by the posting guide). summary(d) y x1 x2 x3 0:54 0:50 0:64 0:78 1:46 1:50 1:36 1:22 fit - glm(y ~ x1 + x2 + x3, data=d, family=binomial(link=logit)) summary(fit) Call: glm(formula = y ~ x1 + x2 + x3, family = binomial(link = logit), data = d) Deviance Residuals: Min 1Q Median 3Q Max -1.6503 -1.0220 -0.7284 0.9965 1.7069 Coefficients: Estimate Std. Error z value Pr(|z|) (Intercept) -0.3772 0.3721 -1.014 0.3107 x11 -0.8144 0.4422 -1.842 0.0655 . x21 0.9226 0.4609 2.002 0.0453 * x31 1.3347 0.5576 2.394 0.0167 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 137.99 on 99 degrees of freedom Residual deviance: 120.65 on 96 degrees of freedom AIC: 128.65 Number of Fisher Scoring iterations: 4 exp(fit$coef) (Intercept) x11 x21 x31 0.6858006 0.4429233 2.5157321 3.7989873 --- After reading the appropriate sections in MASS4 (7.2 and 8.4 in particular), I decided to estimate the 95% confidence intervals for the odds ratios using the profile method implemented in the confint function. I then used the anova function to perform the deviance chi-square tests for each covariate. --- ci - confint(fit); exp(ci) Waiting for profiling to be done... 2.5 %97.5 % (Intercept) 0.3246680 1.413684 x11 0.1834819 1.048154 x21 1.0256096 6.314473 x31 1.3221533 12.129210 anova(fit, test='Chisq') Analysis of Deviance Table Model: binomial, link: logit Response: y Terms added sequentially (first to last) Hence, your use of the `anova' function doesn't return tests corresponding to the CIs computed above. Df Deviance Resid. Df Resid. Dev P(|Chi|) NULL99137.989 x115.85698132.133 0.016 x215.27197126.862 0.022 x316.21296120.650 0.013 My question relates to the interpretation of the significance of variable x1. The OR for x1 is 0.443 and its profile confidence interval is 0.183-1.048. If a type I error rate of 5% is assumed, this result would tend to suggest that x1 is NOT a significant predictor of y. This is also suggested by the Wald test computed by the `summary' function. However, the deviance chi-square test has a P-value of 0.016, which suggests that x1 is indeed a significant predictor of y. How do I That p-value corresponds to adding x1 to a model containing only the intercept term. reconcile these two differing messages? I do recognize that the upper Generally, in order to compute the LR test for the null hypothesis of some subset of the parameters being equal to zero, you need to explicitly fit both the restricted and the unrestricted model and compare them using the `anova' function. Also, see FAQ 7.18. HTH, Henric bound of the confidence interval is pretty close to 1, but I am certain that some journal reviewer will point out the problem as inconsistent. Brant Inman __ 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] Correspondence Analysis
Dear All, I am in the process of teaching myself R and am getting the hang of it slowly, and so apologies for what may be a novice question. (Many thanks to those who have helped me so far). I have been performing normal Correspondence Analysis for a number of years using a variety of packages. CA is available in a number of R packages (e.g. vegan). Could anyone advise me which one(s) would provide me with the standard Greenacre diagnostic statistics (ie. quality, mass, contribution etc as discussed in Correspondence Analysis in Practice, esp. chapter 11) and how I get hold of them? Many thanks in advance, Kris Lockyear. _ Be the first to hear what's new at MSN - sign up to our free newsletters! __ 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] Building R from source
Hello, I was trying to build R from source on Windows XP. I installed software which are mentioned from the follow web page http://www.murdoch-sutherland.com/Rtools/ (Last accessed on Nov. 13th, 2006) . Unfortunately, I got error messages whenever I run 'make all recommended' without modifying 'MkRules' file. I have removed software and reinstalled them several times but I still failed to get it done. The below message is what I got. If anybody gives information, I would appreciate it very much. Thanks, Yongwan === D:\Rsource\R-2.4.0\src\gnuwin32make all recommended make[1]: `Rpwd.exe' is up to date. make[4]: Nothing to be done for `svnonly'. installing C headers make[2]: `all' is up to date. make[2]: `libRblas.dll.a' is up to date. make[5]: Nothing to be done for `svnonly'. installing C headers make --no-print-directory -C ../extra/intl -f Makefile.win make --no-print-directory -C ../appl OPTFLAGS='-O3 -Wall -pedantic -std=gnu99' FOPTFLAGS='-O3 -Wall' -f Makefile.win make --no-print-directory -C ../nmath OPTFLAGS='-O3 -Wall -pedantic -std=gnu99' -f Makefile.win make --no-print-directory -C ../main OPTFLAGS='-O3 -Wall -pedantic -std=gnu99 -DLEA_MALLOC' FFLAGS='-O3 -Wall' -f Makef ile.win make --no-print-directory -C ./graphapp OPTFLAGS='-O3 -Wall -pedantic -std=gnu99' make --no-print-directory -C ./getline OPTFLAGS='-O3 -Wall -pedantic -std=gnu99' make[4]: `gl.a' is up to date. make -f Makefile.win chartables.h make[5]: `chartables.h' is up to date. make -f Makefile.win makeMakedeps make -f Makefile.win libpcre.a make[5]: `libpcre.a' is up to date. make[4]: Nothing to be done for `all'. make[4]: Nothing to be done for `all'. gcc -shared -s -mwindows -o R.dll R.def console.o dataentry.o dynload.o edit.o editor.o embeddedR.o extra.o opt.o pager .o preferences.o psignal.o rhome.o rui.o run.o shext.o sys-win32.o system.o e_pow.o malloc.o ../main/libmain.a ../appl/l ibappl.a ../nmath/libnmath.a graphapp/ga.a getline/gl.a ../extra/xdr/libxdr.a ../extra/zlib/libz.a ../extra/pcre/libpcre .a ../extra/bzip2/libbz2.a ../extra/intl/libintl.a ../extra/trio/libtrio.a dllversion.o -L. -lg2c -lRblas -lcomctl32 -lv ersion ../extra/pcre/libpcre.a(pcre_dfa_exec.o):pcre_dfa_exec.c:(.text+0x15f9): undefined reference to `_pcre_ucp_findprop' ../extra/pcre/libpcre.a(pcre_dfa_exec.o):pcre_dfa_exec.c:(.text+0x1639): undefined reference to `_pcre_ucp_findprop' ../extra/pcre/libpcre.a(pcre_dfa_exec.o):pcre_dfa_exec.c:(.text+0x1736): undefined reference to `_pcre_ucp_findprop' ../extra/pcre/libpcre.a(pcre_dfa_exec.o):pcre_dfa_exec.c:(.text+0x1770): undefined reference to `_pcre_ucp_findprop' ../extra/pcre/libpcre.a(pcre_dfa_exec.o):pcre_dfa_exec.c:(.text+0x184d): undefined reference to `_pcre_ucp_findprop' ../extra/pcre/libpcre.a(pcre_dfa_exec.o):pcre_dfa_exec.c:(.text+0x1e0f): more undefined references to `_pcre_ucp_findpro p' follow collect2: ld returned 1 exit status make[3]: *** [R.dll] Error 1 make[2]: *** [../../bin/R.dll] Error 2 make[1]: *** [rbuild] Error 2 make: *** [all] Error 2 __ 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] Problem with file size
Hi everyone, I have 2 environments (2 different R sessions) as described below: Session 1: Name of the environment: CrlmmInfo Objects in the environment: index1: logical index - length 238304 index2: logical index - length 238304 priors: list of 4 - (matrix 6x6, 2 vectors of length 6, vector of length 2) - all num params: list of 4: centers [238304 x 3 x 2]: num scales [238304 x 3 x 2]: num N [238304 x 3]: num f0 [scalar]: num If I save this environment to a file, I get a file of 23MB. Great. Session 2: Analogous to Session 1, but replace 238304 by 262264. If I save the environment on Session 2, I get a file of 8.4GB. I applied object.size on each of the objects in each environment, and this is what I got: For Session 1: index1: 16204864 index2: 16204864 priors: 3336 params: 74353584 For Session 2: index1: 1049096 index2: 1049096 priors: 3336 params: 81829104 Is this increase from 23MB to 8.4GB expected to happen? Benilton sessionInfo() on both sessions is identical: sessionInfo() R version 2.4.0 (2006-10-03) x86_64-unknown-linux-gnu locale: LC_CTYPE=en_US.iso885915;LC_NUMERIC=C;LC_TIME=en_US.iso885915;LC_COLLATE =en_US.iso885915;LC_MONETARY=en_US.iso885915;LC_MESSAGES=en_US.iso885915 ;LC_PAPER=en_US.iso885915;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASU REMENT=en_US.iso885915;LC_IDENTIFICATION=C attached base packages: [1] methods stats graphics grDevices utils datasets [7] base version _ platform x86_64-unknown-linux-gnu arch x86_64 os linux-gnu system x86_64, linux-gnu status major 2 minor 4.0 year 2006 month 10 day03 svn rev39566 language R version.string R version 2.4.0 (2006-10-03) __ 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] Installing package rpvm under Windows
Hello, I'm trying to install the rpvm package under Windows, but I am having problems. I have pvm3.4 installed properly. I've defined the system variables PVM_ROOT = C:\PROGRA~1\pvm3.4\ PVM_ARCH = win32 When I try to install, I get this: C:\R\PackagesRcmd INSTALL rpvm_1.0.1.tar.gz -- Making package rpvm ** WARNING: this package has a configure script It probably needs manual configuration ** adding build stamp to DESCRIPTION making DLL ... gcc -shared -s -o rpvm.dll rpvm.def rpvm_core.o rpvm_ser.o utils.orpvm_res.o -Lc:/PROGRA~1/R/R-24~1.0/bin -lR rpvm_core.o:rpvm_core.c:(.text+0x44): undefined reference to `pvm_perror' rpvm_core.o:rpvm_core.c:(.text+0x49): undefined reference to `pvm_exit' rpvm_core.o:rpvm_core.c:(.text+0x197): undefined reference to `pvm_mytid' rpvm_core.o:rpvm_core.c:(.text+0x1c7): undefined reference to `pvm_parent' rpvm_core.o:rpvm_core.c:(.text+0x207): undefined reference to `pvm_exit' rpvm_core.o:rpvm_core.c:(.text+0x2a8): undefined reference to `pvm_pstat' rpvm_core.o:rpvm_core.c:(.text+0x352): undefined reference to `pvm_kill' rpvm_core.o:rpvm_core.c:(.text+0x3ba): undefined reference to `pvm_tasks' rpvm_core.o:rpvm_core.c:(.text+0x6bf): undefined reference to `pvm_spawn' rpvm_core.o:rpvm_core.c:(.text+0x7b7): undefined reference to `pvm_initsend' rpvm_core.o:rpvm_core.c:(.text+0x7f7): undefined reference to `pvm_mkbuf' rpvm_core.o:rpvm_core.c:(.text+0x837): undefined reference to `pvm_freebuf' rpvm_core.o:rpvm_core.c:(.text+0x867): undefined reference to `pvm_getsbuf' rpvm_core.o:rpvm_core.c:(.text+0x897): undefined reference to `pvm_getrbuf' rpvm_core.o:rpvm_core.c:(.text+0x8d7): undefined reference to `pvm_setsbuf' rpvm_core.o:rpvm_core.c:(.text+0x917): undefined reference to `pvm_setrbuf' rpvm_core.o:rpvm_core.c:(.text+0x97d): undefined reference to `pvm_pkint' rpvm_core.o:rpvm_core.c:(.text+0x9f2): undefined reference to `pvm_pkdouble' rpvm_core.o:rpvm_core.c:(.text+0xa43): undefined reference to `pvm_pkstr' rpvm_core.o:rpvm_core.c:(.text+0xa9e): undefined reference to `pvm_pkint' rpvm_core.o:rpvm_core.c:(.text+0xad4): undefined reference to `pvm_pkint' rpvm_core.o:rpvm_core.c:(.text+0xb2e): undefined reference to `pvm_pkint' rpvm_core.o:rpvm_core.c:(.text+0xb5f): undefined reference to `pvm_pkdouble' rpvm_core.o:rpvm_core.c:(.text+0xbc1): undefined reference to `pvm_pkint' rpvm_core.o:rpvm_core.c:(.text+0xbf8): undefined reference to `pvm_pkstr' rpvm_core.o:rpvm_core.c:(.text+0xcc5): undefined reference to `pvm_pkint' rpvm_core.o:rpvm_core.c:(.text+0xcf5): undefined reference to `pvm_pkint' rpvm_core.o:rpvm_core.c:(.text+0xdc5): undefined reference to `pvm_pkint' rpvm_core.o:rpvm_core.c:(.text+0xdf5): undefined reference to `pvm_pkdouble' rpvm_core.o:rpvm_core.c:(.text+0xe6e): undefined reference to `pvm_send' rpvm_core.o:rpvm_core.c:(.text+0xee4): undefined reference to `pvm_mcast' rpvm_core.o:rpvm_core.c:(.text+0xf43): undefined reference to `pvm_recv' rpvm_core.o:rpvm_core.c:(.text+0xf9e): undefined reference to `pvm_nrecv' rpvm_core.o:rpvm_core.c:(.text+0xffe): undefined reference to `pvm_probe' rpvm_core.o:rpvm_core.c:(.text+0x107c): undefined reference to `pvm_trecv' rpvm_core.o:rpvm_core.c:(.text+0x113d): undefined reference to `pvm_bufinfo' rpvm_core.o:rpvm_core.c:(.text+0x1217): undefined reference to `pvm_upkint' rpvm_core.o:rpvm_core.c:(.text+0x12ac): undefined reference to `pvm_upkdouble' rpvm_core.o:rpvm_core.c:(.text+0x130b): undefined reference to `pvm_upkstr' rpvm_core.o:rpvm_core.c:(.text+0x1374): undefined reference to `pvm_upkint' rpvm_core.o:rpvm_core.c:(.text+0x13b4): undefined reference to `pvm_upkstr' rpvm_core.o:rpvm_core.c:(.text+0x1421): undefined reference to `pvm_upkint' rpvm_core.o:rpvm_core.c:(.text+0x146f): undefined reference to `pvm_upkint' rpvm_core.o:rpvm_core.c:(.text+0x14c1): undefined reference to `pvm_upkint' rpvm_core.o:rpvm_core.c:(.text+0x150f): undefined reference to `pvm_upkdouble' rpvm_core.o:rpvm_core.c:(.text+0x1561): undefined reference to `pvm_upkint' rpvm_core.o:rpvm_core.c:(.text+0x15be): undefined reference to `pvm_upkint' rpvm_core.o:rpvm_core.c:(.text+0x1611): undefined reference to `pvm_upkint' rpvm_core.o:rpvm_core.c:(.text+0x166e): undefined reference to `pvm_upkdouble' rpvm_core.o:rpvm_core.c:(.text+0x16c0): undefined reference to `pvm_config' rpvm_core.o:rpvm_core.c:(.text+0x1894): undefined reference to `pvm_start_pvmd' rpvm_core.o:rpvm_core.c:(.text+0x191c): undefined reference to `pvm_addhosts' rpvm_core.o:rpvm_core.c:(.text+0x19c7): undefined reference to `pvm_delhosts' rpvm_core.o:rpvm_core.c:(.text+0x1a27): undefined reference to `pvm_halt' rpvm_core.o:rpvm_core.c:(.text+0x1acc): undefined reference to `pvm_mstat' rpvm_core.o:rpvm_core.c:(.text+0x1b53): undefined reference to `pvm_joingroup' rpvm_core.o:rpvm_core.c:(.text+0x1ba3): undefined reference to
Re: [R] Using lrm
Hello Nitin, if you examine the help information for lrm carefully, at the bottom you will find numerous examples that you can follow. ?lrm By the way, asking a queestion like this it's best to clarify if you mean the lrm from the Design package or some other one. I assume it's from Design, as the syntax you quote is a good match. It's also best to study the help files carefully. Andrew On Mon, Nov 13, 2006 at 11:01:59PM -0600, nitin jindal wrote: Hi, I have to build a logistic regression model on a data set that I have. I have three input variables (x1, x2, x3) and one output variable (y). The syntax of lrm function looks like this lrm(formula, data, subset, na.action=na.delete, method=lrm.fit, model=FALSE, x=FALSE, y=FALSE, linear.predictors=TRUE, se.fit=FALSE, penalty=0, penalty.matrix, tol=1e-7, strata.penalty=0, var.penalty=c('simple','sandwich'), weights, normwt, ...) Any logistic regression model would take y and x1,x2,x3 as parameters and output the model (probabilities). So, I dont know where to fit in these values in this function. Any help is appreciated. I am chasing a deadline in my project. Nitin [[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. -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/ __ 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] 回覆: Re: Need help in waveslim package: imodwt and universal.thresh.modwt
Hi Rogerio: I am using Waveslim 1.5 on R 2.3.0, running on Chinese WinXP (SP2). The data, attached in the CSV file, is a stock data (0001 from Hong Kong) downloaded from Yahoo!Finance. The time series data are the Adj. Close prices (last column) from 4-Jan-00 to 30-Nov-2005. The waveslim mra rountines are working fine on these data. Thks rdporto1 [EMAIL PROTECTED] 說: The first two commands worked fine to me. I used xdata = rnorm(128). To let us help you more, specify your R version, some data, OS etc, as usually asked. Rogerio Porto. -- Cabe蓷lho original --- De: [EMAIL PROTECTED] Para: r-help@stat.math.ethz.ch C鏕ia: Data: Sun, 12 Nov 2006 16:02:45 +0800 (CST) Assunto: [R] Need help in waveslim package: imodwt and universal.thresh.modwt Hi: I have encountered problems with imodwt and universal.thresh.modwt and cannot find any reference in R Search. Hope someone can give me some ideas: Starting with modwt.la8 - modwt(xdata, la8, n.level=6) -- this seems to work fine (1) ydata - imodwt(modwt.la8) will always give ydata as numeric(0) (no values) instead of being a time series data with the same lenght as xdata. (2) thred.la8 - universal.thresh.modwt(modwt.la8, max.level=4, hard=F) will give me error at: abs(wc.fine) - cannot contain non-numeric field. Any help is much appreciated. Regards Content-Type: text/html; charset=big5 Content-Transfer-Encoding: 8bit DIVHinbsp; Rogerio:/DIV DIVnbsp;/DIV DIVI am using Waveslim 1.5 on R 2.3.0, running on Chinese WinXP (SP2).nbsp; /DIV DIVnbsp;/DIV DIVThe data, attached in the CSV file, is a stock data (0001 from Hong Kong) downloaded from Yahoo!Finance.nbsp; The time series data are thenbsp;Adj. Close prices (last column) fromnbsp;4-Jan-00 to 30-Nov-2005./DIV DIVnbsp;/DIV DIVThe waveslimnbsp;mra rountines are working fine on these data./DIV DIVnbsp;/DIV DIVThksBRBRBIrdporto1 lt;[EMAIL PROTECTED]gt;/I/B 說:/DIV BLOCKQUOTE class=replbq style=PADDING-LEFT: 5px; MARGIN-LEFT: 5px; BORDER-LEFT: #1010ff 2px solidThe first two commands worked fine to me. I usedBRxdata = rnorm(128).BRBRTo let us help you more, specify your R version,BRsome data, OS etc, as usually asked.BRBRRogerio Porto.BRBR-- Cabe蓷lho original ---BRBRDe: [EMAIL PROTECTED]BRPara: r-help@stat.math.ethz.chBRC鏕ia:BRData: Sun, 12 Nov 2006 16:02:45 +0800 (CST)BRAssunto: [R] Need help in waveslim package: imodwt and universal.thresh.modwtBRBRgt; Hi:BRgt; I have encountered problems with imodwt and universal.thresh.modwt and cannot find any reference in R Search. Hope someone can give me some ideas:BRgt;BRgt; Starting withBRgt; modwt.la8 lt;- modwt(xdata, la8, n.level=6) lt;-- this seems to work fineBRgt;BRgt; (1) ydata lt;- imodwt(modwt.la8)BRgt; will always give ydata as numeric(0) (no values) instead of being a time series data with the same lenght as xdata.BRgt;BRgt; (2) thred.la8 lt;- universal.thresh.modwt(modwt.la8, max.level=4, hard=F) will give me error at:BRgt; abs(wc.fine) - cannot contain non-numeric field.BRgt;BRgt; Any help is much appreciated.BRgt;BRgt; RegardsBRgt;BRgt;BRBR/BLOCKQUOTEBRp#32;___br YM - 離線訊息br Content-Type: application/octet-stream; name=CKH-0001-ALL-data-only.csv Content-Transfer-Encoding: base64 Content-Description: 3176266068-CKH-0001-ALL-data-only.csv Content-Disposition: attachment; filename=CKH-0001-ALL-data-only.csv MDg2MTEsNzkuNjUNCg== __ 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] Problem with file size
On Mon, 13 Nov 2006, Benilton Carvalho wrote: Hi everyone, I have 2 environments (2 different R sessions) as described below: Session 1: Name of the environment: CrlmmInfo Objects in the environment: index1: logical index - length 238304 index2: logical index - length 238304 priors: list of 4 - (matrix 6x6, 2 vectors of length 6, vector of length 2) - all num params: list of 4: centers [238304 x 3 x 2]: num scales [238304 x 3 x 2]: num N [238304 x 3]: num f0 [scalar]: num If I save this environment to a file, I get a file of 23MB. Great. Session 2: Analogous to Session 1, but replace 238304 by 262264. If I save the environment on Session 2, I get a file of 8.4GB. I applied object.size on each of the objects in each environment, and this is what I got: For Session 1: index1: 16204864 index2: 16204864 priors: 3336 params: 74353584 For Session 2: index1: 1049096 index2: 1049096 priors: 3336 params: 81829104 Is this increase from 23MB to 8.4GB expected to happen? We don't have enough information to expect anything. Saving environments is a tricky business: read the description in the R Internals manual if this is news to you. But note that the index[12] objects are much smaller in Session 2, and those in Session 1 are about 16x larger than needed for the description given. So my guess is that they have attributes which are promises in Session 2, and save() there is pulling in another environment to enable those attributes to be evaluated later. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ 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] Embedded carriage returns in text document
For OS X, http://dos2unix.darwinports.com/ -- Jeffrey R. Spies http://www.nd.edu/~jspies/ On Nov 13, 2006, at 11:06 PM, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote: Dennis, You can get rid of the '^M' by reformating the file from DOS to UNIX. Withing a UNIX system send the command dos2unix filename1 filename2. Filename2 won't have the '^M'. Hope it helps, Augusto Augusto Sanabria. MSc, PhD. Mathematical Modeller Risk Research Group Geospatial Earth Monitoring Division Geoscience Australia (www.ga.gov.au) Cnr. Jerrabomberra Av. Hindmarsh Dr. Symonston ACT 2601 Ph. (02) 6249-9155 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Dennis Fisher Sent: Tuesday, 14 November 2006 9:06 To: r-help@stat.math.ethz.ch Subject: [R] Embedded carriage returns in text document Colleagues, I am using R 2.4.0 on both a Mac (10.4.8) and Linux (RedHat 9). To read data from an Excel spreadsheet, I do save as in Excel, then select the Text (tab-delimited) format. The resulting file uses a tab separator and I can usually read the file using read.delim. Sometimes, the header row contains embedded carriage returns. When I view the file, these carriage returns appear as ^M. Now the problem: When I read.delim these files, they do not read correctly. Sometimes I get error messages; sometimes only the first line is read. Interestingly, invoking the option skip=1 (or a larger N) does not appear to bypass the problem. I can solve the problem by manually deleting these carriage returns either in the original Excel file or the .txt version. However, this is not an ideal solution. Does anyone have a work-around within R? Dennis Dennis Fisher MD P (The P Less Than Company) Phone: 1-866-PLessThan (1-866-753-7784) Fax: 1-415-564-2220 www.PLessThan.com [[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. [[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] Variance of a complex estimator using survey package ...
Hi, I want to compute the variance of two complex statistics. The first statistic is the a ratio R1=Q(a1)/Q(a2) where Q denote de quantile at a1 and a2. The second is also a ratio but not a classic one this ratio is R2=sum(x_{i}|x_{i}Q(a1))/sum(x_{i} }|Q(a2)) Linearisation for those statistic is very complex. I want to use bootstrap and JKn How can I procede since I dont have a function like svrepanyfunction where you just specify the function to compute the statistic. Sincerly. Justin BEM Elève Ingénieur Statisticien Economiste BP 294 Yaoundé. Tél (00237)9597295. ___ Découvrez une nouvelle façon d'obtenir des réponses à toutes vos questions ! Profitez des connaissances, des opinions et des expériences des internautes sur Yahoo! Questions/Réponses [[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] reconstruction of data after PCA and interpolation
Dear all, I have a matrix of size N x M. I purchase a PCA analysis through prcomp function. Then I keep the H eigenvectors which explain 90% of the total variance and interpolate each vectors of the matrix H x M, to obtain a new matrix of size H x K (K M). My question is : from this last matrix (H x K), how to (re)construct a matrix of size N x K which have the same units of the inital matrix (N x M). Regards -- Cordialement Emmanuel Poizot Cnam/Intechmer B.P. 324 50103 Cherbourg Cedex Phone (Direct) : (00 33)(0)233887342 Fax : (00 33)(0)233887339 __ 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] Installing package rpvm under Windows
Plese heed the warning enclosed in ***. You haven't set up the pvm libraries in the link. If you don't know how to do the manual configuration, talk to the package maintainer (as the posting guide asked you to do before posting). On Mon, 13 Nov 2006, Adrian Dragulescu wrote: I'm trying to install the rpvm package under Windows, but I am having problems. I have pvm3.4 installed properly. I've defined the system variables PVM_ROOT = C:\PROGRA~1\pvm3.4\ Using \ is unlikely to work with the tools that we use, but I see no sign that this is being used. PVM_ARCH = win32 When I try to install, I get this: C:\R\PackagesRcmd INSTALL rpvm_1.0.1.tar.gz -- Making package rpvm ** WARNING: this package has a configure script It probably needs manual configuration ** adding build stamp to DESCRIPTION making DLL ... gcc -shared -s -o rpvm.dll rpvm.def rpvm_core.o rpvm_ser.o utils.orpvm_res.o -Lc:/PROGRA~1/R/R-24~1.0/bin -lR rpvm_core.o:rpvm_core.c:(.text+0x44): undefined reference to `pvm_perror' [...] -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ 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.