gmane.comp.lang.r.general

Yet another set of codes to optimize

Daren Tan — 2008-12-05T02:41:23

I have problems converting my dataset from long to wide format. Previous attempts using reshape package and aggregate function were unsuccessful as they took too long. Apparently, my simplified solution also lasted as long. My complete codes is given below. When sample.size = 10000, the execution takes about 20 seconds. But sample.size = 100000 seems to take eternity. My actual sample.size is 15000000 i.e. 15 million. sample.size <- 10000 m <- data.frame(Name=sample(1:100000, sample.size, T), Type=sample(1:1000, sample.size, T), Predictor=sample(LETTERS[1:10], sample.size, T)) res <- function(m) { m.12.unique <- unique(m[,1:2]) m.12.unique <- m.12.unique[order(m.12.unique[,1], m.12.unique[,2]),] v1 <- paste(m.12.unique[,1], m.12.unique[,2], sep=".") v2 <- c(sort(unique(m[,3]))) res <- matrix(0, nr=length(v1), nc=length(v2), dimnames=list(v1, v2)) m.ids <- paste(m[,1], m[,2], sep=".") for(i in 1:nrow(m)) { x <- m.ids[i] y <- m[i,3] res[x, y] <- res[x,

making sense of posterior statistics in the deal package

Aaron Tarone — 2008-12-05T00:35:05

Hello, I'm doing bayesian network analyses with the deal package. I am at a loss for how to interpret output from the analysis (i.e. what is a good score, what is a bad score, which stats tell me what about the network edges/nodes). Here is an example node with its posterior scores for all parent nodes. ------------------------------------------------------------ Conditional Posterior: Yp1| 3 4 5 6 9 11 12 15 18 [[1]] [[1]]$tau [,1] [,2] [,3] [,4] [,5] [,6] [1,] 138.0000000 -201.944190 -61.827901 -29.5419149 11.7780877 -56.1691436 [2,] -201.9441898 379.014299 101.336606 49.2886631 -9.5976678 99.0119458 [3,] -61.8279013 101.336606 55.301879 18.3175413 0.4718180 31.7741275 [4,] -29.5419149 49.288663 18.317541 18.5074653 0.7297184 14.7963722 [5,] 11.7780877 -9.597668 0.471818 0.7297184 11.9705940 -0.1152971 [6,] -56.1691436 99.011946 31.774127 14.7963722 -0.1152971 33.0750507 [7,] 11.8398168 -11.819652 2

Re: How to optimize this codes ?

Daren Tan — 2008-12-05T01:17:43

Hi William, Thanks for your suggestion. I learn alot from it, some concepts are new which need some digestion. > Subject: Re: [R] How to optimize this codes ?> Date: Thu, 4 Dec 2008 11:49:30 -0800> From: wdunlap< at >tibco.com> To: daren76< at >hotmail.com> CC: R-help< at >stat.math.ethz.ch> > > -----Original Message-----> > From: William Dunlap > > Sent: Thursday, December 04, 2008 9:59 AM> > To: 'daren76< at >hotmail.com'> > Cc: 'R help'> > Subject: Re: [R] How to optimize this codes ?> > > > [R] How to optimize this codes ?> > Daren Tan daren76 at hotmail.com> > Thu Dec 4 17:02:49 CET 2008> > > > How to optimize the for-loop to be reasonably fast for > > sample.size=100000000 ? > > You may want to change sample.size=1000 to have an idea > > what I am achieving. > > > > set.seed(143)> > A <- matrix(sample(0:1, sample.s ize, TRUE), ncol=10, > > dimnames=list(NULL, LETTERS[1:10]))> > > > B <- list()> > for(i in 1:10) { > > B[[i]] <- apply(combn(LETTERS[1:10], i), 2, function(x) > > { sum(apply(data.frame(A[! ,x]), 1, all)

Re: R] adding a new dataset to the default R distribution

John Maindonald — 2008-12-05T00:33:04

Making data, especially data that have been the subject of published papers, widely available, can be a useful spinoff from the R project, another gift to the scientific community beyond the provision of computing and analytic tools. Nowadays, in a complete publication of a scientific result, there is every reason for the data to be part of that publication. The Gentleman and Lang 2004 paper "Statistical Analyses and Reproducible Research" takes this further still, making a compelling case for opening the analysis to ready view. (http://www.bepress.com/bioconductor/paper2/ ) How else can critics know what analysis was done, and whether the data do really support the claimed conclusions? As I see it, the first recourse should be use of archives that individual communities may establish. Instructions on how to input the data into R would be a useful small item of ancillary information. Links to such archives (under Data Archives, maybe) might be included on CRAN. The Open Archae

Question about lrandom effects specification in lme4 :Correction

Bert Gunter — 2008-12-05T00:27:46

... Of course what I should have said is: lmer( y ~ ... + (-1+A|Grp) + (-1+B|Grp) + (-1+C|Grp)) Cheers, Bert Gunter Genentech ______________________________________________ R-help< at >r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help< at >r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.

Re: Formula parsing and updating

Charles C. Berry — 2008-12-05T00:21:41

You need to give us something to reproduce. I get an error when I try Error in model.matrix.default(newTerms, m) : model frame and formula mismatch in model.matrix() so I do not see how you obtained 'somemod'. When I protect with I(quant>0), I have no problems obtaining 'fit' nor update()ing it: Call: coxph(formula = Surv(time, status) ~ I(quant > 0) + strata(sex), data = test1) coef exp(coef) se(coef) z p I(quant > 0)TRUE -0.746 0.474 1.27 -0.589 0.56 Likelihood ratio test=0.37 on 1 df, p=0.545 n=6 (1 observation deleted due to missingness) HTH, Chuck Charles C. Berry (858) 534-2098 Dept of Family/Preventive Medicine E mailto:cberry< at >tajo.ucsd.edu UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 ______________________________________________ R-help< at >r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE

Question about lrandom effects specification in lme4

Bert Gunter — 2008-12-05T00:21:17

Folks: Suppose I have 3 random effects, A,B, and C. Using the older lme() function (in nlme) it was possible (using the pdMat classes) to specify that they are uncorrelated with identical variances. Is it possible to do this with lmer? My understanding is that if I specify them as lmer( y ~ ... + (A|Grp) + (B|Grp) + (C|Grp)) then they are uncorrelated but have different variances. Motivation: I'd like to use lmer instead of lme for fitting smoothing splines to longitudinal data. Cheers, Bert Gunter Genentech ______________________________________________ R-help< at >r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.

integration within maximum likelihood

Cristian Angelo Guevara — 2008-12-04T23:29:14

Hi: I'm trying to estimate a latent variable model in mnl discrete choice framework using R. I need to do first a uni dimensional integral within each observation (row) in the database and then sum over observations. I'm stacked in the point shown below. Apparently I have a dimensionality problem in the definition of the integral. Maybe it does not identify that what I need is only one row at a time. I'm no sure. I would appreciate any kind of help. Below is my code and the error report. Thanks Angelo + th1<-theta[1] + th2<-theta[2] + tha<-theta[3] + thb<-theta[4] + thc<-theta[5] + thp<-theta[6] + thmu<-theta[7] + alfz<-theta[8] + alfp<-theta[9] + mu1<- alfz*y$z_a1 + alfp*y$p1 + mu2<- alfz*y$z_a2 + alfp*y$p2 + mu3<- alfz*y$z_a3 + alfp*y$p3 + U1<- th1 +tha*y$a1 +thb*y$b1 +thc*y$c1 +thp*y$p1 +thmu*mu1 + U2<- th2 +tha*y$a2 +thb*y$b2 +thc*y$c2 +thp*y$p2 +thmu*mu2 + U3<- +tha*y$a3 +thb*y$b3 +thc*y$c3 +thp*y$p3 +thmu*mu3 + Usum<- (exp(U1)+exp(U2)+exp(U3)) + arg1<- function(mu1) {(exp(U1)/Usum)*dnorm(m

Re: How to get Greenhouse-Geisser epsilons from anova?

Peter Dalgaard — 2008-12-04T22:23:32

Why? It is supposed to work. You just need to work out the X and M specification for the relevant error strata and set test="Spherical" for anova.mlm, or work out the T contrast matrix explicitly if that suits your temper better.

Re: Running R under Sun Grid Engine with OpenMPI tightintegration

Sean Davis — 2008-12-04T22:11:13

No, I haven't gotten things to work. Even under openmpi without SGE, I have R crashes on a regular basis when using Rmpi-0.5.5, so I haven't really pursued it further. Sean ______________________________________________ R-help< at >r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.

Plotting a kriging on a map

Rodrigo Aluizio — 2008-12-04T21:27:35

Hi list, Well, this time Ive a doubt with mapping generation. I was already able to read and plot shapefiles, plot point on this map. All this without any problems, but now I want to do something that I think, only Golden Software Surfer is capable of. I would like to plot a kriging result on the existing map (map script below). Well I looked for on the net, but I was not able to be sure of which packages Ill need to perform the kriging (maybe geoR?) and then plot it in the area map. Is R able to do that? What packages I need to know to do so? I just need a starting point to go forward. I need to generate something like this: http://i288.photobucket.com/albums/ll162/Godrigos/Kriging.jpg (but without the top right area, just the different grayscale kriging area). library(xlsReadWrite) library(maps) library(mapdata) library(maptools) Pontos<-read.xls('Pontos.xls',sheet=1,rowNames=T) Prof<-read.xls('Pontos.xls',sheet=2,rowNames=T) Batimetria<-readShapeSpatial(Batimetria.shp') Munici

comparing SAS and R survival analysis with time-dependentcovariates

Svetlana Eden — 2008-12-04T21:05:33

Dear R-help, I was comparing SAS (I do not know what version it is) and R (version 2.6.0 (2007-10-03) on Linux) survival analyses with time-dependent covariates. The results differed significantly so I tried to understand on a short example where I went wrong. The following example shows that even when argument 'method' in R function coxph and argument 'ties' in SAS procedure phreg are the same, the results of Cox regr. are different. This seems to happen when there are ties in the events/covariates times. My question is what software, R or SAS, is more reliable for the survival analysis with time-dependent covariates or if you could point out a problem in the following example. Example. SAS gives HR=3.236: data trythis; input id days timedeli stat; datalines; 1 3 .5 1 2 1.5 1 1 3 6 1000 0 4 8 1000 1 5 8 1 0 6 21 1000 1 7 11 3 1 run; proc phreg data=trythis; model days*stat(0)=de

Formula parsing and updating

Johannes Huesing — 2008-12-04T21:00:07

Hi all, I can't come over a problem with formula. Suppose I have a coxmod model with the following formula: Surv(lebzeit, tot == 1) ~ sex + (alter >= 65) + diff3k + zelltyp_k_c + q_nuc_3k + kar_k80_g80 + stadium and I want to drop the stadium explanatory variable from the model with I get the following messages: Fehler in coxph(formula = Surv(lebzeit, tot == 1) ~ sex + alter >= 65 + : No (non-missing) observations Zusätzlich: Warning messages: 1: In Ops.factor(sex, alter) : + not meaningful for factors 2: In Ops.factor(65, diff3k) : + not meaningful for factors 3: In Ops.factor(65 + diff3k, zelltyp_k_c) : + not meaningful for factors 4: In Ops.factor(65 + diff3k + zelltyp_k_c, q_nuc_3k) : + not meaningful for factors 5: In Ops.factor(65 + diff3k + zelltyp_k_c + q_nuc_3k, kar_k80_g80) : + not meaningful for factors What seems to happen is that "alter >= 65" is not parsed as an entity. It does not look as a problem with the environment to me, as is outlined in http://developer.r-project.o

Re: How to optimize this codes ?

William Dunlap — 2008-12-04T19:49:30

A hastily written recursive version is: computeB3 <- function(A) { B <- rep(list(integer(0)), ncol(A)) storage.mode(A) <- "logical" recurse <- function(thisCol = 1, includedCols = rep(NA, ncol(A)), allSoFar = rep(TRUE, nrow(A))) { if (thisCol < ncol(A)) # skip this column Recall(thisCol+1, replace(includedCols, thisCol, FALSE), allSoFar = allSoFar) else if (thisCol > ncol(A)) return() includedCols[thisCol] <- TRUE allSoFar <- allSoFar & A[,thisCol] nIncludedCols <- sum(includedCols[1:thisCol]) # cat(thisCol, ":", nIncludedCols, ":", includedCols[1:thisCol], sum(allSoFar), "\n") # Note B<<- in next line uses lexical scoping to change computeB3:B. Does not work in S+. B[[nIncludedCols]][length(B[[nIncludedCols]])+1] <<- sum(allSoFar) Recall(thisCol+1, includedCols, allSoFar = allSoFar) } recurse() B } The elements of B[[i]], for each i, are in a different order than they are in the other functions, but the histogram

Approximate matching: agrep and TRE?

Stavros Macrakis — 2008-12-04T19:03:36

agrep{base} is very useful, but it would be very useful to also have the offset of the first match returned, e.g. agrepx("foo","sdxforzoo",max.distance=1) => 4 (there is room for quibbling over the definition of first match...). As far as I can tell, the only way to do that now is to do multiple agreps, something like min(agrep("foo",substr(rep("foobar",...))). Is there a version of R agrep that instead just returns the position of first match? The TRE package (http://www.laurikari.net/tre/) implements a much more general version of agrep (including such functionality) -- has anyone integrated it into R? Thanks, -s [[alternative HTML version deleted]] ______________________________________________ R-help< at >r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.

How to get Greenhouse-Geisser epsilons from anova?

Skotara — 2008-12-04T18:44:35

Dear all, I apologize for my basic question. I try to calculate an anova for repeated measurements with 3 factors (A,B,C) having 2, 2, and 7 levels. or with an additional fourth between subjects factor D. Everything works fine using aov(val ~ A*B*C + Error(subject/ (A*B*C) ) ) or aov(val ~ (D*A*B*C) + Error(subject/(A*B*C)) + D ) val, A, B, C, D and subject are columns in a data.frame. How can I get the estimated Greenhouse-Geisser and Huynh-Feldt epsilons? I know Peter Dalgaard described it in R-News Vol. 7/2, October 2007. However, unfortunately I am not able to apply that using my data... Furthermore, I am still confused of how SPSS calculates the epsilons since it is mentioned that perhaps there are any errors in SPSS?? I would be glad if anyone could help me! I am looking forward to hearing from you! Thank you! Nils ______________________________________________ R-help< at >r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-projec

Re: RODBC - problems connecting to oracle through linux

Simon Collins — 2008-12-04T18:43:51

Just an update to this - it is possible to connect to an 11G database via using the oracle 10G (libsqora.so.10.1) driver. Quoting "Prof Brian Ripley" stats.ox.ac.uk>: ______________________________________________ R-help< at >r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.

Re: Writing a list help

2008-12-04T18:35:49

You can get all the rownames by calling write.table once on something like # if all list elements have equal length t(data.frame(MyList)) # otherwise, let NA extend where needed t(do.call(cbind, lapply(MyList, ts))) You can take care of NAs as needed on your end. HTH -----Original Message----- From: r-help-bounces< at >r-project.org [mailto:r-help-bounces< at >r-project.org] On Behalf Of Rajasekaramya Sent: Thursday, December 04, 2008 10:52 AM To: r-help< at >r-project.org Subject: [R] Writing a list help hi there, I have a list called MyList. MyList[[1]] [1] "SRY" [2] "RPS4Y1" [3] "ZFY" [4] "ZFX /// ZFY" [5] "LOC728137" [6] "LOC100101116 /// TTTY1" [7] "AMELY" [8] "TBL1Y" [9] "PRKY" [[2]] . . . [[24]] I have list.I wanna write it to file. for(i in 1:24) { write.table(as.vector(as.matrix(MyLIST[[i]] file="geneset.txt",sep="\t",append=T) } I used a for loop to write my list.It worked but i have two concerns with it. First while writing in a txt file it writing in the names verticaly

Re: readAsciiGrid memory issues

Roger Bivand — 2008-12-04T18:25:07

runs out of memory. I have 4Gb ram This function is written for small rasters only, and is not suited to data of this size. that enable to read in large Use rgdal, preferably using a tiled approach, leaving at least 75% of your memory free initially, because operations on the tile will result in copies. See several threads on R-sig-geo on handling large data sets, including one at: https://stat.ethz.ch/pipermail/r-sig-geo/2008-November/004486.html Also visit the raster package on R-Forge, which may offer a packaged solution to this. Roger Bivand ______________________________________________ R-help< at >r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.

Re: Simulating underdispersed counts

Greg Snow — 2008-12-04T18:03:36

Yes, that is basically truncation, but the original poster said they wanted something fast and did not give detail. If they just wanted some data to put into a glm model to demonstrate, then this would work. Another idea that I had is instead of dropping all the extreme values, is take the extreme values and replace them with a new draw, this would still allow for the possibility of some extreme values, but would reduce the number of them. This fits the idea "if my results don't match what I expected, something must have gone wrong and I will just do it again" that is sometimes seen in researchers who don't fully understand the idea of variation (Mendel's gardener/assistant as a possible example). Another approach: One way to think of a regular Poisson process is to have several bins and place objects in the bins at random. If the probability of placing something into a bin is independent of how many objects are already in that bin (and the others) then the counts of objects per bin will follow a Poiss

Re: TreeMap Graphs

jim holtman — 2008-12-04T17:59:10

Checkout the portfolio package. On Thu, Dec 4, 2008 at 9:44 AM, rbs.com> wrote: ______________________________________________ R-help< at >r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.