Re: [R] set seed for random draws
On Sun, 6 Nov 2011, R. Michael Weylandt wrote: I think it's easier than you are making it: the random seed is created in a "pretty-random" way when you first use it and then it is updated Ah: It is unless you then save the workspace. If you do, then evey subsequent session starts with the same seed until you save the workspace again. So never saving and always saving works fine, but occasional saving can lead to puzzlement. That "pretty-random" way is as unpredictable as pseado-random numbers (It uses a PRNG internally.) with each call to rDIST(). For example, set.seed(1) x1 <- .Random.seed rnorm(1) x2 <- .Random.seed rnorm(1) x3 <- .Random.seed identical(x1, x2) FALSE identical(x1, x3) FALSE identical(x2, x3) FALSE set.seed(1) identical(x1, .Random.seed) TRUE rnorm(2) identical(x3, .Random.seed) TRUE But the period for the random seed to repeat is very, very long so you don't have to think about it unless you really need to (or for reproducible simulations) Michael On Sat, Nov 5, 2011 at 7:22 PM, Md Desa, Zairul Nor Deana Binti wrote: Thank you everybody for the helpful advices. Basically, I try to figure out why I get different numbers as there are more than one seed for a loop within a loop. Well, I guest I got it now. Because every time random seed is called or specified it'll output different random numbers, as it's requested. Thanks! D From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] on behalf of Rolf Turner [rolf.tur...@xtra.co.nz] Sent: Saturday, November 05, 2011 3:22 PM To: Patrick Burns Cc: r-help@r-project.org; achim.zeil...@uibk.ac.at Subject: Re: [R] set seed for random draws On 05/11/11 22:00, Patrick Burns wrote: I'd suggest two rules of thumb when coming up against something in R that you aren't sure about: 1. If it is a mundane task, R probably takes care of it. 2. Experiment to see what happens. Of course you could read documentation, but no one does that. Fortune nomination! cheers, Rolf __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@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@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. -- Brian D. Ripley, rip...@stats.ox.ac.uk 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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] install.packages problem
This is something missing from your (unstated) Linux installation. curl-config is part of the original libcurl sources, but Linux distributors tend to separte it out. *How* they do so is non-standard: Fedora and other RPM-based distributions tend to use libcurl-devel Debian and related tend to use libcurl-dev You need to figure this out for your distribution and install the missing piece. On Sat, 5 Nov 2011, eric wrote: I'm trying to install the rdatamarket package. I did an install.packages('rdatamarket') command but got an error about half way through the install as follows: * installing *source* package ‘RCurl’ ... checking for curl-config... no Cannot find curl-config ERROR: configuration failed for package ‘RCurl’ The install continued after the error but looks like it was completed. I'm trying to figure out what the error means and how I fix it. Here's what I'm seeing ...ideas on how to address this would be appreciated : install.packages('rdatamarket') Installing package(s) into ‘/home/eric/R/i486-pc-linux-gnu-library/2.13’ (as ‘lib’ is unspecified) --- Please select a CRAN mirror for use in this session --- also installing the dependencies ‘RCurl’, ‘RJSONIO’ trying URL 'http://lib.stat.cmu.edu/R/CRAN/src/contrib/RCurl_1.7-0.tar.gz' Content type 'application/x-gzip' length 813252 bytes (794 Kb) opened URL == downloaded 794 Kb trying URL 'http://lib.stat.cmu.edu/R/CRAN/src/contrib/RJSONIO_0.96-0.tar.gz' Content type 'application/x-gzip' length 1144519 bytes (1.1 Mb) opened URL == downloaded 1.1 Mb trying URL 'http://lib.stat.cmu.edu/R/CRAN/src/contrib/rdatamarket_0.6.3.tar.gz' Content type 'application/x-gzip' length 12432 bytes (12 Kb) opened URL == downloaded 12 Kb * installing *source* package ‘RCurl’ ... checking for curl-config... no Cannot find curl-config ERROR: configuration failed for package ‘RCurl’ * removing ‘/home/eric/R/i486-pc-linux-gnu-library/2.13/RCurl’ * installing *source* package ‘RJSONIO’ ... Trying to find libjson.h header file checking for gcc... gcc checking whether the C compiler works... yes checking for C compiler default output file name... a.out checking for suffix of executables... checking whether we are cross compiling... no checking for suffix of object files... o checking whether we are using the GNU C compiler... yes checking whether gcc accepts -g... yes checking for gcc option to accept ISO C89... none needed USE_LOCAL = "" Using local libjson code. Copying files /tmp/RtmpFw9QeX/R.INSTALL4ebf657f/RJSONIO configure: creating ./config.status config.status: creating src/Makevars config.status: creating cleanup ** libs gcc -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -std=gnu99 -O3 -pipe -g -c ConvertUTF.c -o ConvertUTF.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONChildren.cpp -o JSONChildren.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONDebug.cpp -o JSONDebug.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONIterators.cpp -o JSONIterators.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONMemory.cpp -o JSONMemory.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONNode.cpp -o JSONNode.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONNode_Mutex.cpp -o JSONNode_Mutex.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONStream.cpp -o JSONStream.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONValidator.cpp -o JSONValidator.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONWorker.cpp -o JSONWorker.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONWriter.cpp -o JSONWriter.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSON_Base64.cpp -o JSON_Base64.o gcc -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -std=gnu99 -O3 -pipe -g -c JSON_parser.c -o JSON_parser.o gcc -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -std=gnu99 -O3 -pipe -g -c RJSON.c -o RJSON.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_
Re: [R] Correlation between matrices
Hi: I don't think you want to keep these objects separate; it's better to combine everything into a data frame. Here's a variation of your example - the x variable ends up being a mouse, but you may have another variable that's more appropriate to plot so take this as a starting point. One plot uses the ggplot2 package, the other uses the lattice and latticeExtra packages. library('ggplot2') regions = c('cortex', 'hippocampus', 'brain_stem', 'mid_brain', 'cerebellum') mice = paste('mouse', 1:5, sep='') elem <- c('Cu', 'Fe', 'Zn', 'Ca', 'Enzyme') # Generate a data frame from the combinations of # mice, regions and elem: d <- data.frame(expand.grid(mice = mice, regions = regions, elem = elem), y = rnorm(125)) # Create a numeric version of mice d$mouse <- as.numeric(d$mice) # A function to return regression coefficients coefun <- function(df) coef(lm(y ~ mouse), data = df) # Apply to all regions * elem combinations coefs <- ddply(d, .(regions, elem), coefun) names(coefs) <- c('regions', 'elem', 'b0', 'b1') # Generate the plot using package ggplot2: ggplot(d, aes(x = mouse, y = y)) + geom_point(size = 2.5) + geom_abline(data = coefs, aes(intercept = b0, slope = b1), size = 1) + facet_grid(elem ~ regions) # Same plot in lattice: library('lattice') library('latticeExtra') p <- xyplot(y ~ mouse | elem + regions, data = d, type = c('p', 'r'), layout = c(5, 5)) HTH, Dennis On Sat, Nov 5, 2011 at 10:49 AM, Kaiyin Zhong wrote: >> regions = c('cortex', 'hippocampus', 'brain_stem', 'mid_brain', > 'cerebellum') >> mice = paste('mouse', 1:5, sep='') >> for (n in c('Cu', 'Fe', 'Zn', 'Ca', 'Enzyme')) { > + assign(n, as.data.frame(replicate(5, rnorm(5 > + } >> names(Cu) = names(Zn) = names(Fe) = names(Ca) = names(Enzyme) = regions >> row.names(Cu) = row.names(Zn) = row.names(Fe) = row.names(Ca) = > row.names(Enzyme) = mice >> Cu > cortex hippocampus brain_stem mid_brain cerebellum > mouse1 -0.5436573 -0.31486713 0.1039148 -0.3908665 -1.0849112 > mouse2 1.4559136 1.75731752 -2.1195118 -0.9894767 0.3609033 > mouse3 -0.6735427 -0.04666507 0.9641000 0.4683339 0.7419944 > mouse4 0.6926557 -0.47820023 1.3560802 0.9967562 -1.3727874 > mouse5 0.2371585 0.20031393 -1.4978517 0.7535148 0.5632443 >> Zn > cortex hippocampus brain_stem mid_brain cerebellum > mouse1 -0.66424043 0.6664478 1.1983546 0.0319403 0.41955740 > mouse2 -1.14510448 1.5612235 0.3210821 0.4094753 1.01637466 > mouse3 -0.85954416 2.8275458 -0.6922565 -0.8182307 -0.06961242 > mouse4 0.03606034 -0.7177256 0.7067217 0.2036655 -0.25542524 > mouse5 0.67427572 0.6171704 0.1044267 -1.8636174 -0.07654666 >> Fe > cortex hippocampus brain_stem mid_brain cerebellum > mouse1 1.8337008 2.0884261 0.29730413 -1.6884804 0.8336137 > mouse2 -0.2734139 -0.5728439 0.63791556 -0.6232828 -1.1352224 > mouse3 -0.4795082 0.1627235 0.21775206 1.0751584 -0.5581422 > mouse4 1.7125147 -0.5830600 1.40597896 -0.2815305 0.3776360 > mouse5 -0.3469067 -0.4813120 -0.09606797 1.0970077 -1.1234038 >> Ca > cortex hippocampus brain_stem mid_brain cerebellum > mouse1 -0.7663354 0.8595091 1.33803798 -1.17651576 0.8299963 > mouse2 -0.7132260 -0.2626811 0.08025079 -2.40924271 0.7883005 > mouse3 -0.7988904 -0.1144639 -0.65901136 0.42462227 0.7068755 > mouse4 0.3880393 0.5570068 -0.49969135 0.06633009 -1.3497228 > mouse5 1.0077684 0.6023264 -0.57387762 0.25919461 -0.9337281 >> Enzyme > cortex hippocampus brain_stem mid_brain cerebellum > mouse1 1.3430936 0.5335819 -0.56992947 1.3565803 -0.8323391 > mouse2 1.0520850 -1.0201124 0.8965 1.4719880 1.0854768 > mouse3 -0.2802482 0.6863323 -1.37483570 -0.7790174 0.2446761 > mouse4 -0.1916415 -0.4566571 1.93365932 1.3493848 0.2130424 > mouse5 -1.0349593 -0.1940268 -0.07216321 -0.2968288 1.7406905 > > In each anatomic region, I would like to calculate the correlation between > Enzyme activity and each of the concentrations of Cu, Zn, Fe, and Ca, and > do a scatter plot with a tendency line, organizing those plots into a grid. > See the image below for the desired effect: > http://postimage.org/image/62brra6jn/ > How can I achieve this? > > Thank you in advance. > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] set seed for random draws
I think it's easier than you are making it: the random seed is created in a "pretty-random" way when you first use it and then it is updated with each call to rDIST(). For example, set.seed(1) x1 <- .Random.seed rnorm(1) x2 <- .Random.seed rnorm(1) x3 <- .Random.seed identical(x1, x2) FALSE identical(x1, x3) FALSE identical(x2, x3) FALSE set.seed(1) identical(x1, .Random.seed) TRUE rnorm(2) identical(x3, .Random.seed) TRUE But the period for the random seed to repeat is very, very long so you don't have to think about it unless you really need to (or for reproducible simulations) Michael On Sat, Nov 5, 2011 at 7:22 PM, Md Desa, Zairul Nor Deana Binti wrote: > Thank you everybody for the helpful advices. > Basically, I try to figure out why I get different numbers as there are more > than one seed for a loop within a loop. Well, I guest I got it now. Because > every time random seed is called or specified it'll output different random > numbers, as it's requested. > > Thanks! > > D > > From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] on behalf > of Rolf Turner [rolf.tur...@xtra.co.nz] > Sent: Saturday, November 05, 2011 3:22 PM > To: Patrick Burns > Cc: r-help@r-project.org; achim.zeil...@uibk.ac.at > Subject: Re: [R] set seed for random draws > > On 05/11/11 22:00, Patrick Burns wrote: > > >> I'd suggest two rules of thumb when coming >> up against something in R that you aren't >> sure about: >> >> 1. If it is a mundane task, R probably >> takes care of it. >> >> 2. Experiment to see what happens. >> >> >> Of course you could read documentation, but >> no one does that. > > > > Fortune nomination! > > cheers, > > Rolf > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > > > __ > R-help@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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Matrix element-by-element multiplication
There are a few (nasty?) side-effects to c(), one of which is stripping a matrix of its dimensionality. E.g., x <- matrix(1:4, 2) c(x) [1] 1 2 3 4 So that's probably what happened to you. R has a somewhat odd feature of not really considering a pure vector as a column or row vector but being willing to change it to either: e.g. y <- 1:2 x %*% y y %*% x y %*% y while matrix(y) %*% x throws an error, which can also trip folks up. You might also note that x * y and y*x return the same thing in this problem. Getting back to your problem: what are v and b and what are you hoping to get done? Specifically, what happened when you tried v*b (give the exact error message). It seems likely that they are non-conformable matrices, but here non-conformable for element-wise multiplication doesn't mean the same thing as it does for matrix multiplication. E.g., x <- matrix(1:4,2) y <- matrix(1:6,2) dim(x) [1] 2 2 dim(y) [1] 2 3 x * y -- here R seems to want matrices with identical dimensions, but i can't promise that. x %*% y does work. Hope this helps and yes I know it can seem crazy at first, but there really is reason behind it at the end of the tunnel, Michael On Sun, Nov 6, 2011 at 12:11 AM, Steven Yen wrote: > My earlier attempt > > dp <- v*b > > did not work. Then, > > dp <- c(v)*b > > worked. > > Confused, > > Steven > > At 09:10 PM 11/4/2011, you wrote: > > Did you even try? > > a <- 1:3 > x <- matrix(c(1,2,3,2,4,6,3,6,9),3) > a*x > > [,1] [,2] [,3] > [1,] 1 2 3 > [2,] 4 8 12 > [3,] 9 18 27 > > Michael > > On Fri, Nov 4, 2011 at 7:26 PM, Steven Yen wrote: >> is there a way to do element-by-element multiplication as in Gauss >> and MATLAB, as shown below? Thanks. >> >> --- >> a >> >> 1.000 >> 2.000 >> 3.000 >> x >> >> 1.000 2.000 3.000 >> 2.000 4.000 6.000 >> 3.000 6.000 9.000 >> a.*x >> >> 1.000 2.000 3.000 >> 4.000 8.000 12.00 >> 9.000 18.00 27.00 >> >> >> -- >> Steven T. Yen, Professor of Agricultural Economics >> The University of Tennessee >> http://web.utk.edu/~syen/ >> [[alternative HTML version deleted]] >> >> __ >> R-help@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide >> http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. >> > > -- > Steven T. Yen, Professor of Agricultural Economics > The University of Tennessee > http://web.utk.edu/~syen/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nested "for" loops
No idea how this relates to what you said originally but glad you got it all worked out. And let us all reiterate: really, don't use nested for loops...there's a better way: promise! Michael On Sat, Nov 5, 2011 at 2:20 PM, nick_pan wrote: > I found the way out - it was because the borders of the vectors was close > enough thats why I had the same result while I was adding points to the > sequence. The example I gave was irrelevant but I made in order to find out > that the problem was. > Thank you all for your answers. > > > -- > View this message in context: > http://r.789695.n4.nabble.com/nested-for-loops-tp3992089p3993917.html > Sent from the R help mailing list archive at Nabble.com. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] install.packages problem
On Sat, Nov 5, 2011 at 2:47 PM, eric wrote: > I'm trying to install the rdatamarket package. I did an > install.packages('rdatamarket') command but got an error about half way > through the install as follows: > > * installing *source* package ‘RCurl’ ... > checking for curl-config... no > Cannot find curl-config > ERROR: configuration failed for package ‘RCurl’ > > The install continued after the error but looks like it was completed. I'm > trying to figure out what the error means and how I fix it. I think you are on Linux and missing the library headers to build against the curl. On Fedora or RedHat you would do something like 'sudo yum -y install curl-devel' and on Ubuntu it may be 'sudo apt-get install curl-dev' Andrew __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Doing dist on separate objects in a text file
Perhaps split() directly or more abstractly tapply() from base or one of the d_ply() from plyr? Michael On Sat, Nov 5, 2011 at 7:20 PM, ScottDaniel wrote: > So I have a text file that looks like this: > "Label" "X" "Y" "Slice" > 1 "Field_1_R3D_D3D_PRJ_w617.tif" 348 506 1 > 2 "Field_1_R3D_D3D_PRJ_w617.tif" 359 505 1 > 3 "Field_1_R3D_D3D_PRJ_w617.tif" 356 524 1 > 4 "Field_1_R3D_D3D_PRJ_w617.tif" 2 0 1 > 5 "Field_1_R3D_D3D_PRJ_w617.tif" 412 872 1 > 6 "Field_1_R3D_D3D_PRJ_w617.tif" 422 863 1 > 7 "Field_1_R3D_D3D_PRJ_w617.tif" 429 858 1 > 8 "Field_1_R3D_D3D_PRJ_w617.tif" 429 880 1 > 9 "Field_1_R3D_D3D_PRJ_w617.tif" 437 865 1 > 10 "Field_1_R3D_D3D_PRJ_w617.tif" 447 855 1 > 11 "Field_1_R3D_D3D_PRJ_w617.tif" 450 868 1 > 12 "Field_1_R3D_D3D_PRJ_w617.tif" 447 875 1 > 13 "Field_1_R3D_D3D_PRJ_w617.tif" 439 885 1 > 14 "Field_1_R3D_D3D_PRJ_w617.tif" 2 8 1 > > What it represents are the locations of centromeres per nucleus in a > microscope image. What I need to do is do a dist() on each grouping (the > grouping being separated by the low values of x and y's) and then compute an > average. The part that I'm having trouble with is writing code that will > allow R to separate these objects. Do I have to find some way of creating > separate data frames for each object? Or is there a way to parse the file > and generate a single data frame of all the pairwise distances? Any > suggestions or example code would be much appreciated. Thanks! > > -- > View this message in context: > http://r.789695.n4.nabble.com/Doing-dist-on-separate-objects-in-a-text-file-tp3994515p3994515.html > Sent from the R help mailing list archive at Nabble.com. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nested "for" loops
Bert, this is not helpful. Since for loops and apply functions are not vectorized, why are you admonishing Carl that vectorizing doesn't always speed up algorithms? He didn't reference apply functions as being vectorized. But you seem to be doing so. I would assert that vectorizing DOES always speed up algorithms, but things people sometimes think are vectorizing are not really. --- Jeff Newmiller The . . Go Live... DCN: Basics: ##.#. ##.#. Live Go... Live: OO#.. Dead: OO#.. Playing Research Engineer (Solar/Batteries O.O#. #.O#. with /Software/Embedded Controllers) .OO#. .OO#. rocks...1k --- Sent from my phone. Please excuse my brevity. Bert Gunter wrote: Carl: "Almost anything you can do in a for() loop can be done easier and faster with vectorization.-- " That is false: while this is certainly true for a great many basic vectorized operations, it is certainly false for most other things -- simulations are a typical example. Note that __ply type operations in base R, plyr , and other packages are (generally) _not_vectorized; they are "disguised" loops. Explicit for() loops are often just as fast or even a bit faster, though many of us prefer what we consider the more transparent functional style code of the __ply's. Cheers, Bert Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm [[alternative HTML version deleted]] _ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] linear against nonlinear alternatives - quantile regression
Dear David, Indeed rq() accepts a vector fo tau. I used the example given by Frank to run fitspl4 <- summary(rq(b1 ~ rcs(x,4), tau=c(a1,a2,a3,a4))) and it works. I even can use anova() to test equality of slopes jointly across quantiles. however, it would be interesting to test among different specifications, e.g. rcs(x,4) against rcs(x,3). but it does not work. Thanks for all suggestions! Julia > From: dwinsem...@comcast.net > Date: Sat, 5 Nov 2011 13:42:34 -0400 > To: f.harr...@vanderbilt.edu > CC: r-help@r-project.org > Subject: Re: [R] linear against nonlinear alternatives - quantile regression > > I suppose this constitutes thread drift, but your simple example, Frank, made > wonder if Rq() accepts a vector argument for tau. I seem to remember that > Koencker's rq() does.. Normally I would consult the help page, but the power > is still out here in Central Connecticut and I am corresponding with a less > capable device. I am guessing that if Rq() does accept such a vector that the > form of the nonlinearity would be imposed at all levels of tau. > > -- > David > > On Nov 5, 2011, at 10:43 AM, Frank Harrell wrote: > > > Just to address a piece of this - in the case in which you are currently > > focusing on only one quantile, the rms package can help by fitting > > restricted cubic splines for covariate effects, and then run anova to test > > for nonlinearity (sometimes a dubious practice because if you then remove > > nonlinear terms you are mildly cheating). > > > > require(rms) > > f <- Rq(y ~ x1 + rcs(x2,4), tau=.25) > > anova(f) # tests associations and nonlinearity of x2 > > > > Frank > > > > Julia Lira wrote: > >> > >> Dear all, > >> > >> I would like to know whether any specification test for linear against > >> nonlinear model hypothesis has been implemented in R using the quantreg > >> package. > >> > >> I could read papers concerning this issue, but they haven't been > >> implemented at R. As far as I know, we only have two specification tests > >> in this line: anova.rq and Khmaladze.test. The first one test equality and > >> significance of the slopes across quantiles and the latter one test if the > >> linear specification is model of location or location and scale shift. > >> > >> Do you have any suggestion? > >> > >> Thanks a lot! > >> > >> Best regards, > >> > >> Julia > >> > >>[[alternative HTML version deleted]] > >> > >> __ > >> R-help@ 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. > >> > > > > > > - > > Frank Harrell > > Department of Biostatistics, Vanderbilt University > > -- > > View this message in context: > > http://r.789695.n4.nabble.com/linear-against-nonlinear-alternatives-quantile-regression-tp3993327p3993416.html > > Sent from the R help mailing list archive at Nabble.com. > > > > __ > > R-help@r-project.org mailing list > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > > and provide commented, minimal, self-contained, reproducible code. > > __ > R-help@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. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Correlation between matrices
> regions = c('cortex', 'hippocampus', 'brain_stem', 'mid_brain', 'cerebellum') > mice = paste('mouse', 1:5, sep='') > for (n in c('Cu', 'Fe', 'Zn', 'Ca', 'Enzyme')) { + assign(n, as.data.frame(replicate(5, rnorm(5 + } > names(Cu) = names(Zn) = names(Fe) = names(Ca) = names(Enzyme) = regions > row.names(Cu) = row.names(Zn) = row.names(Fe) = row.names(Ca) = row.names(Enzyme) = mice > Cu cortex hippocampus brain_stem mid_brain cerebellum mouse1 -0.5436573 -0.31486713 0.1039148 -0.3908665 -1.0849112 mouse2 1.4559136 1.75731752 -2.1195118 -0.9894767 0.3609033 mouse3 -0.6735427 -0.04666507 0.9641000 0.4683339 0.7419944 mouse4 0.6926557 -0.47820023 1.3560802 0.9967562 -1.3727874 mouse5 0.2371585 0.20031393 -1.4978517 0.7535148 0.5632443 > Zn cortex hippocampus brain_stem mid_brain cerebellum mouse1 -0.66424043 0.6664478 1.1983546 0.0319403 0.41955740 mouse2 -1.14510448 1.5612235 0.3210821 0.4094753 1.01637466 mouse3 -0.85954416 2.8275458 -0.6922565 -0.8182307 -0.06961242 mouse4 0.03606034 -0.7177256 0.7067217 0.2036655 -0.25542524 mouse5 0.67427572 0.6171704 0.1044267 -1.8636174 -0.07654666 > Fe cortex hippocampus brain_stem mid_brain cerebellum mouse1 1.8337008 2.0884261 0.29730413 -1.6884804 0.8336137 mouse2 -0.2734139 -0.5728439 0.63791556 -0.6232828 -1.1352224 mouse3 -0.4795082 0.1627235 0.21775206 1.0751584 -0.5581422 mouse4 1.7125147 -0.5830600 1.40597896 -0.2815305 0.3776360 mouse5 -0.3469067 -0.4813120 -0.09606797 1.0970077 -1.1234038 > Ca cortex hippocampus brain_stem mid_brain cerebellum mouse1 -0.7663354 0.8595091 1.33803798 -1.17651576 0.8299963 mouse2 -0.7132260 -0.2626811 0.08025079 -2.40924271 0.7883005 mouse3 -0.7988904 -0.1144639 -0.65901136 0.42462227 0.7068755 mouse4 0.3880393 0.5570068 -0.49969135 0.06633009 -1.3497228 mouse5 1.0077684 0.6023264 -0.57387762 0.25919461 -0.9337281 > Enzyme cortex hippocampus brain_stem mid_brain cerebellum mouse1 1.3430936 0.5335819 -0.56992947 1.3565803 -0.8323391 mouse2 1.0520850 -1.0201124 0.8965 1.4719880 1.0854768 mouse3 -0.2802482 0.6863323 -1.37483570 -0.7790174 0.2446761 mouse4 -0.1916415 -0.4566571 1.93365932 1.3493848 0.2130424 mouse5 -1.0349593 -0.1940268 -0.07216321 -0.2968288 1.7406905 In each anatomic region, I would like to calculate the correlation between Enzyme activity and each of the concentrations of Cu, Zn, Fe, and Ca, and do a scatter plot with a tendency line, organizing those plots into a grid. See the image below for the desired effect: http://postimage.org/image/62brra6jn/ How can I achieve this? Thank you in advance. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] ANESRAKE package: Inappropriate error message, given the data
Hi, Attempting to use the "anesrake" raking (also referred to as rim weighting) package to weight survey data (Note: DATA is listed at the bottom along with SESSION INFO AND COMMANDS ARGUMENTS USED to allow you to reproduce the error). Presents following error message: *"Error in selecthighestpcts(discrep1, inputter, pctlim) : * * No variables are off by more than 5 percent using the method you have chosen, either weighting is unnecessary or a smaller pre-raking limit should be chosen."* However the sample proportions for the demographics are in fact far more than 5% off the target proportions (see data below) no matter what choosemethod argument is used. Have tried different "pctlim" settings as well (even small ones such as 0.0001) and still get same message. Also tried different "choosemethod" - both "max" and "total". Tried to get help from the following but see nothing on this error in these references: *Package description at: cran.r-project.org/web/packages/anesrake/anesrake .pdf* *Paper by Josh Pasek: http://www.stanford.edu/~jpasek/Josh_Pasek/Software_files/RakingDescription.pdf * Here is the* session info*: > sessionInfo() R version 2.14.0 (2011-10-31) Platform: x86_64-pc-mingw32/x64 (64-bit) locale: [1] LC_COLLATE=English_South Africa.1252 LC_CTYPE=English_South Africa.1252 [3] LC_MONETARY=English_South Africa.1252 LC_NUMERIC=C [5] LC_TIME=English_South Africa.1252 attached base packages: [1] splines stats graphics grDevices utils datasets methods [8] base other attached packages: [1] anesrake_0.70Hmisc_3.8-3 survival_2.36-10 loaded via a namespace (and not attached): [1] cluster_1.14.1 grid_2.14.0lattice_0.20-0 HERE ARE THE ARGUMENTS > medufo <- read.csv("C:/Users/lenovo/Desktop/Medufo_weights/medufo.csv", header = TRUE, row.names="ID") > age <- c(25, 25, 28, 22) > gender <- c(49, 51) > ethnic <- c(58, 14, 5, 23) > income <- c(1, 4, 7, 15, 18, 16, 19, 17, 1, 1,1) > targets <- list(age, gender, ethnic, income) > names(targets) <- c("age", "gender", "ethnic", "income") > outsave <- anesrake(targets, medufo, caseid=1, verbose = TRUE) > write.csv(print(outsave), "C:/Users/lenovo/Desktop/Medufo_weights/outsave.csv") HERE IS THE DATA LISTING > list(medufo) [[1]] age gender ethnic income 202 1 1 7 233 1 4 6 282 1 1 10 312 1 1 8 324 1 4 8 351 1 3 8 383 1 4 9 403 1 4 9 432 1 1 8 463 1 3 6 482 1 1 4 583 1 4 8 632 1 3 8 643 1 1 7 722 1 1 8 813 1 4 8 884 1 4 6 893 1 1 8 904 1 4 8 931 1 1 8 943 1 4 8 982 1 1 8 104 2 1 1 10 109 3 1 2 7 110 1 1 2 7 114 1 1 1 5 119 1 1 4 4 124 1 1 4 7 132 3 1 4 8 141 3 1 4 10 142 4 1 4 7 3 2 2 1 10 5 1 2 2 4 6 4 2 1 8 8 2 2 3 6 9 2 2 1 8 103 2 4 8 112 2 4 8 133 2 1 7 142 2 4 8 152 2 1 8 162 2 4 8 192 2 4 8 213 2 1 7 223 2 1 6 244 2 4 8 254 2 4 5 263 2 4 8 274 2 4 8 303 2 4 10 334 2 2 5 342 2 2 5 363 2 4 7 372 2 1 4 393 2 4 7 412 2 4 8 422 2 3 1 453 2 2 8 473 2 2 5 493 2 4 7 513 2 4 8 521 2 2 4 533 2 4 8 544 2 4 4 552 2 4 11 562 2 4 4 573 2 3 7 592 2 3 7 603 2 4 7 613 2 2 8 652 2 1 8 663 2 2 5 673 2 4 9 681 2 4 7 691 2 1 4 703 2 4 8 732 2 1 5 742 2 2 8 752 2 1 5 764 2 1 7 772 2 4 7 782 2 4 7 791 2 2 6 802 2 1 6 824 2 4 8 833 2 4 8 852 2 1 6 862 2 4 8 872 2 1 3 913 2 4 8 922 2 4 8 954 2 4 7 962 2 1 9 973 2 4
[R] Doing dist on separate objects in a text file
So I have a text file that looks like this: "Label" "X" "Y" "Slice" 1 "Field_1_R3D_D3D_PRJ_w617.tif" 348 506 1 2 "Field_1_R3D_D3D_PRJ_w617.tif" 359 505 1 3 "Field_1_R3D_D3D_PRJ_w617.tif" 356 524 1 4 "Field_1_R3D_D3D_PRJ_w617.tif" 2 0 1 5 "Field_1_R3D_D3D_PRJ_w617.tif" 412 872 1 6 "Field_1_R3D_D3D_PRJ_w617.tif" 422 863 1 7 "Field_1_R3D_D3D_PRJ_w617.tif" 429 858 1 8 "Field_1_R3D_D3D_PRJ_w617.tif" 429 880 1 9 "Field_1_R3D_D3D_PRJ_w617.tif" 437 865 1 10 "Field_1_R3D_D3D_PRJ_w617.tif" 447 855 1 11 "Field_1_R3D_D3D_PRJ_w617.tif" 450 868 1 12 "Field_1_R3D_D3D_PRJ_w617.tif" 447 875 1 13 "Field_1_R3D_D3D_PRJ_w617.tif" 439 885 1 14 "Field_1_R3D_D3D_PRJ_w617.tif" 2 8 1 What it represents are the locations of centromeres per nucleus in a microscope image. What I need to do is do a dist() on each grouping (the grouping being separated by the low values of x and y's) and then compute an average. The part that I'm having trouble with is writing code that will allow R to separate these objects. Do I have to find some way of creating separate data frames for each object? Or is there a way to parse the file and generate a single data frame of all the pairwise distances? Any suggestions or example code would be much appreciated. Thanks! -- View this message in context: http://r.789695.n4.nabble.com/Doing-dist-on-separate-objects-in-a-text-file-tp3994515p3994515.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nested "for" loops
I found the way out - it was because the borders of the vectors was close enough thats why I had the same result while I was adding points to the sequence. The example I gave was irrelevant but I made in order to find out that the problem was. Thank you all for your answers. -- View this message in context: http://r.789695.n4.nabble.com/nested-for-loops-tp3992089p3993917.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] install.packages problem
I'm trying to install the rdatamarket package. I did an install.packages('rdatamarket') command but got an error about half way through the install as follows: * installing *source* package ‘RCurl’ ... checking for curl-config... no Cannot find curl-config ERROR: configuration failed for package ‘RCurl’ The install continued after the error but looks like it was completed. I'm trying to figure out what the error means and how I fix it. Here's what I'm seeing ...ideas on how to address this would be appreciated : install.packages('rdatamarket') Installing package(s) into ‘/home/eric/R/i486-pc-linux-gnu-library/2.13’ (as ‘lib’ is unspecified) --- Please select a CRAN mirror for use in this session --- also installing the dependencies ‘RCurl’, ‘RJSONIO’ trying URL 'http://lib.stat.cmu.edu/R/CRAN/src/contrib/RCurl_1.7-0.tar.gz' Content type 'application/x-gzip' length 813252 bytes (794 Kb) opened URL == downloaded 794 Kb trying URL 'http://lib.stat.cmu.edu/R/CRAN/src/contrib/RJSONIO_0.96-0.tar.gz' Content type 'application/x-gzip' length 1144519 bytes (1.1 Mb) opened URL == downloaded 1.1 Mb trying URL 'http://lib.stat.cmu.edu/R/CRAN/src/contrib/rdatamarket_0.6.3.tar.gz' Content type 'application/x-gzip' length 12432 bytes (12 Kb) opened URL == downloaded 12 Kb * installing *source* package ‘RCurl’ ... checking for curl-config... no Cannot find curl-config ERROR: configuration failed for package ‘RCurl’ * removing ‘/home/eric/R/i486-pc-linux-gnu-library/2.13/RCurl’ * installing *source* package ‘RJSONIO’ ... Trying to find libjson.h header file checking for gcc... gcc checking whether the C compiler works... yes checking for C compiler default output file name... a.out checking for suffix of executables... checking whether we are cross compiling... no checking for suffix of object files... o checking whether we are using the GNU C compiler... yes checking whether gcc accepts -g... yes checking for gcc option to accept ISO C89... none needed USE_LOCAL = "" Using local libjson code. Copying files /tmp/RtmpFw9QeX/R.INSTALL4ebf657f/RJSONIO configure: creating ./config.status config.status: creating src/Makevars config.status: creating cleanup ** libs gcc -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -std=gnu99 -O3 -pipe -g -c ConvertUTF.c -o ConvertUTF.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONChildren.cpp -o JSONChildren.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONDebug.cpp -o JSONDebug.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONIterators.cpp -o JSONIterators.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONMemory.cpp -o JSONMemory.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONNode.cpp -o JSONNode.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONNode_Mutex.cpp -o JSONNode_Mutex.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONStream.cpp -o JSONStream.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONValidator.cpp -o JSONValidator.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONWorker.cpp -o JSONWorker.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSONWriter.cpp -o JSONWriter.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c JSON_Base64.cpp -o JSON_Base64.o gcc -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -std=gnu99 -O3 -pipe -g -c JSON_parser.c -o JSON_parser.o gcc -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -std=gnu99 -O3 -pipe -g -c RJSON.c -o RJSON.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c internalJSONNode.cpp -o internalJSONNode.o g++ -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -O3 -pipe -g -c libjson.cpp -o libjson.o gcc -I/usr/share/R/include -I. -Ilibjson -Ilibjson/Source -DNDEBUG=1 -DJSON_NO_EXCEPTIONS=1 -fpic -std=gnu99 -O3 -pipe -g -c rlibjson.c -o rlibjson.o g++ -shared -o RJSONIO.so ConvertUTF.o JSONChil
Re: [R] nested "for" loops
Carl: "Almost anything you can do in a for() loop can be done easier and faster with vectorization.-- " That is false: while this is certainly true for a great many basic vectorized operations, it is certainly false for most other things -- simulations are a typical example. Note that __ply type operations in base R, plyr , and other packages are (generally) _not_vectorized; they are "disguised" loops. Explicit for() loops are often just as fast or even a bit faster, though many of us prefer what we consider the more transparent functional style code of the __ply's. Cheers, Bert Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] set seed for random draws
Thank you everybody for the helpful advices. Basically, I try to figure out why I get different numbers as there are more than one seed for a loop within a loop. Well, I guest I got it now. Because every time random seed is called or specified it'll output different random numbers, as it's requested. Thanks! D From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] on behalf of Rolf Turner [rolf.tur...@xtra.co.nz] Sent: Saturday, November 05, 2011 3:22 PM To: Patrick Burns Cc: r-help@r-project.org; achim.zeil...@uibk.ac.at Subject: Re: [R] set seed for random draws On 05/11/11 22:00, Patrick Burns wrote: > I'd suggest two rules of thumb when coming > up against something in R that you aren't > sure about: > > 1. If it is a mundane task, R probably > takes care of it. > > 2. Experiment to see what happens. > > > Of course you could read documentation, but > no one does that. Fortune nomination! cheers, Rolf __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@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] testing significance of axis loadings from multivariate dudi.mix
Hi all I´m trying to tests the significance of loadings from a ordination of 46 variables (caategorical, ordinal and nominal). I used dudi.mix from ade4 for the ordination. A years ago Jari Oksanen wrote this script implementing Peres-Neto et al. 2003 (Ecology) bootstraping method: netoboot <- function (x, permutations=1000, ...) { pcnull <- princomp(x, cor = TRUE, ...) res <- pcnull$loadings out <- matrix(0, nrow=nrow(res), ncol=ncol(res)) N <- nrow(x) for (i in 1:permutations) { pc <- princomp(x[sample(N, replace=TRUE), ], cor = TRUE ...) pred <- predict(pc, newdata = x) r <- cor(pcnull$scores, pred) k <- apply(abs(r), 2, which.max) reve <- sign(diag(r[k,])) sol <- pc$loadings[ ,k] sol <- sweep(sol, 2, reve, "*") out <- out + ifelse(res > 0, sol <= 0, sol >= 0) } out/permutations } I tried to aply it to the case of dudi.mix instead of princomp this way: netoboot1<-function (x, permutations=1000,...) { dudinull <- dudi.mix(x, scannf = FALSE, nf = 3) res <- dudinull$c1 out <- matrix(0, nrow=nrow(res), ncol=ncol(res)) N <- nrow(x) for (i in 1:permutations) { dudi <- dudi.mix(x[sample(N, replace=TRUE), ], scannf = FALSE, nf = 3) pred <- predict(dudi, newdata = x) r <- cor(dudinull$li, pred) k <- apply(abs(r), 2, which.max) reve <- sign(diag(r[k,])) sol <- dudi$c1[ ,k] sol <- sweep(sol, 2, reve, "*") out <- out + ifelse(res > 0, sol <= 0, sol >= 0) } out/permutations } But a problem arised with the predict function: it doesn´t seem to work with an object from dudi.mix and I dont understand why. Can somebody tell me why? Any suggestions to modify the script or to use other method? Thanks in advance. Francisco Francisco Mora Ardila Laboratorio de Biodiversidad y Funcionamiento del Ecosistema Centro de Investigaciones en Ecosistemas UNAM-Campus Morelia Tel 3222777 ext. 42621 Morelia , MIchoacán, México. -- Open WebMail Project (http://openwebmail.org) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] set seed for random draws
On 05/11/11 22:00, Patrick Burns wrote: I'd suggest two rules of thumb when coming up against something in R that you aren't sure about: 1. If it is a mundane task, R probably takes care of it. 2. Experiment to see what happens. Of course you could read documentation, but no one does that. Fortune nomination! cheers, Rolf __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nested "for" loops
You need to define "l" as a dimensioned object , either a vector or an array, and then assign the value of your calculation to the correctly indexed "location" in that object. Otherwise you are just overwriting the value each time through the loop. Use these help pages (and review "Introduction to R" ?array ?vector ?"[" On Nov 4, 2011, at 7:49 PM, nick_pan wrote: > Hi all , I have written a code with nested "for" loops . > The aim is to estimate the maximum likelihood by creating 3 vectors with the > same length( sequence ) > and then to utilize 3 "for" loops to make combinations among the 3 vectors , > which are (length)^3 in number , and find the one that maximize the > likelihood ( maximum likelihood estimator). > > > The code I created, runs but I think something goes wrong...because when I > change the length of the vectors but not the bounds the result is the > same!!! > > I will give a simple example(irrelevant but proportional to the above) to > make it more clear... > > Lets say we want to find the combination that maximize the multiplication of > the entries of some vectors. > > V1<-c(1,2,3) > V2<-c(5, 2 , 4) > V3<-c( 4, 3, 6) > > The combination we look for is ( 3 , 5 , 6) that give us 3*5*6 = 90 > > If I apply the following in R , I won't take this result > > V1<-c(1,2,3) > V2<-c(5, 2 , 4) > V3<-c( 4, 3, 6) > > for( i in V1){ > for( j in V2) { > for( k in V3){ > > l<- i*j*k > > } > } > } > l > > Then " l<- i*j*k " is number and not vector(of all multiplications of all > the combinations) , and is 3*4*6 = 72. > > How can I fix the code? > > > > > -- > View this message in context: > http://r.789695.n4.nabble.com/nested-for-loops-tp3992089p3992089.html > Sent from the R help mailing list archive at Nabble.com. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] error message
Humm... I was using adegenet / ade4 packages and both R and R studio prompted the questions. Sorry, I did so many searches on R help and Adegenet help that I forgot to mention the packages I was using... Juliana -- View this message in context: http://r.789695.n4.nabble.com/error-message-tp3223412p3994158.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] error message
What prompted the questions? If it was a specifically biological/genetic package the bioconductor list can probably provide better help than this list. Michael On Nov 5, 2011, at 3:12 PM, JulianaMF wrote: > Hello Michael, > Sorry, I am just starting to lear all this. > Here is one of my input files (from a .str file) in which the first column > are the individuals, the second is the pop info (in this case I am stating > that I have one pop because I am still trying to find out which are the > clusters in my sample) and the others are the 10 loci. The genotypes for > each individual are in 2 rows and missing data is -9. I tried both R and R > studio and I am working from a Mac. I followed all the routine > SimStru <- read.structure(file = "SimStru.str") > and then I answered all the prompted questions (I am writing from my memory > because the files are back at the lab): 82 individuals (820 genotypes), 10 > loci, no column with marker name, column with pop info = 2, other column = 1 > (names for each sample), row with marker names = 1, missing data = -9. I > didn't change any of the other default settings. And once I answered all the > questions, I got that error message. > Thank you so much to be willing to help and sorry about my ignorance, I am > new to this! > Ind Pop 134 220 28 18 414 24 42 58423 12 > S3 1 349 163 267 316 287 412 275 234 164 351 > S3 1 369 165 267 336 287 424 275 238 188 351 > S5 1 345 163 271 316 287 360 187 234 152 343 > S5 1 365 163 283 336 287 388 187 246 152 615 > S9 1 353 163 275 300 287 400 231 234 164 347 > S9 1 361 163 279 336 287 416 275 234 170 351 > S101 325 -9 259 316 287 384 299 234 140 331 > S101 357 -9 279 328 287 400 299 234 152 351 > S151 377 163 267 316 287 416 259 234 134 339 > S151 385 163 283 344 287 416 267 254 164 363 > S171 333 163 263 316 285 380 179 234 164 355 > S171 381 163 287 328 287 400 179 238 164 399 > S211 333 163 271 356 285 388 219 250 158 335 > S211 377 165 271 360 285 416 219 270 200 355 > S221 373 163 251 316 285 404 211 234 158 335 > S221 377 163 279 352 287 404 211 234 158 355 > S261 377 163 259 324 285 424 -9 254 170 327 > S261 405 163 283 324 287 424 -9 254 188 351 > S281 333 163 267 324 287 380 187 246 164 351 > S281 333 163 267 336 287 416 295 246 164 355 > S321 321 163 259 300 285 396 291 250 152 343 > S321 365 165 259 348 285 396 291 250 164 351 > S331 325 163 263 -9 287 408 231 238 140 371 > S331 357 163 263 -9 287 432 251 246 146 371 > S371 361 163 267 320 285 416 195 254 164 343 > S371 361 165 275 320 287 420 195 254 170 355 > S381 377 163 267 348 287 416 223 234 164 339 > S381 385 165 275 388 287 416 223 254 164 363 > S401 373 163 255 336 287 384 191 234 158 347 > S401 381 163 267 348 287 416 239 250 170 495 > S441 333 163 279 300 287 408 255 238 146 351 > S441 333 163 283 316 287 416 255 246 146 387 > S451 345 163 -9 -9 287 412 187 234 158 -9 > S451 389 163 -9 -9 287 416 215 234 176 -9 > S491 349 163 271 312 285 388 191 238 164 351 > S491 357 163 283 344 287 400 191 238 182 495 > S521 353 163 259 300 287 380 195 258 -9 335 > S521 385 163 279 320 287 424 195 258 -9 351 > S561 325 163 267 300 287 400 263 226 146 355 > S561 389 163 283 316 287 400 263 238 188 407 > S571 357 163 263 300 287 412 199 234 146 343 > S571 389 163 271 316 287 412 199 242 176 391 > S611 369 163 287 316 287 380 175 246 164 339 > S611 393 163 287 316 287 432 191 250 182 347 > S621 377 163 275 316 287 392 -9 254 164 363 > S621 385 165 275 348 287 392 -9 254 164 363 > S631 333 163 267 332 285 388 175 234 158 347 > S631 401 163 271 340 287 408 199 234 164 355 > S671 377 163 251 320 287 412 187 246 158 351 > S671 377 163 275 352 287 416 187 246 158 379 > S681 377 163 287 332 287 404 291 238 164 339 > S681 405 163 287 348 287 420 291 246 176 591 > S691 325 163 255 308 285 364 191 226 -9 343 > S691 353 163 279 316 285 416 307 242 -9 359 > S711 385 163 255 328 285 384 251 230 164 343 > S711 385 165 267 352 287 424 251 230 168 363 > S731 341 159 259 308 285 416 191 230 164 335 > S731 369 163 259 324 285 416 191 230 176 467 > S741 -9 -9 271 320 287 384 183 226 152 347 > S741 -9 -9 271 348 287 408 227 22
Re: [R] error message
Hello Michael, Sorry, I am just starting to lear all this. Here is one of my input files (from a .str file) in which the first column are the individuals, the second is the pop info (in this case I am stating that I have one pop because I am still trying to find out which are the clusters in my sample) and the others are the 10 loci. The genotypes for each individual are in 2 rows and missing data is -9. I tried both R and R studio and I am working from a Mac. I followed all the routine SimStru <- read.structure(file = "SimStru.str") and then I answered all the prompted questions (I am writing from my memory because the files are back at the lab): 82 individuals (820 genotypes), 10 loci, no column with marker name, column with pop info = 2, other column = 1 (names for each sample), row with marker names = 1, missing data = -9. I didn't change any of the other default settings. And once I answered all the questions, I got that error message. Thank you so much to be willing to help and sorry about my ignorance, I am new to this! Ind Pop 134 220 28 18 414 24 42 58423 12 S3 1 349 163 267 316 287 412 275 234 164 351 S3 1 369 165 267 336 287 424 275 238 188 351 S5 1 345 163 271 316 287 360 187 234 152 343 S5 1 365 163 283 336 287 388 187 246 152 615 S9 1 353 163 275 300 287 400 231 234 164 347 S9 1 361 163 279 336 287 416 275 234 170 351 S101 325 -9 259 316 287 384 299 234 140 331 S101 357 -9 279 328 287 400 299 234 152 351 S151 377 163 267 316 287 416 259 234 134 339 S151 385 163 283 344 287 416 267 254 164 363 S171 333 163 263 316 285 380 179 234 164 355 S171 381 163 287 328 287 400 179 238 164 399 S211 333 163 271 356 285 388 219 250 158 335 S211 377 165 271 360 285 416 219 270 200 355 S221 373 163 251 316 285 404 211 234 158 335 S221 377 163 279 352 287 404 211 234 158 355 S261 377 163 259 324 285 424 -9 254 170 327 S261 405 163 283 324 287 424 -9 254 188 351 S281 333 163 267 324 287 380 187 246 164 351 S281 333 163 267 336 287 416 295 246 164 355 S321 321 163 259 300 285 396 291 250 152 343 S321 365 165 259 348 285 396 291 250 164 351 S331 325 163 263 -9 287 408 231 238 140 371 S331 357 163 263 -9 287 432 251 246 146 371 S371 361 163 267 320 285 416 195 254 164 343 S371 361 165 275 320 287 420 195 254 170 355 S381 377 163 267 348 287 416 223 234 164 339 S381 385 165 275 388 287 416 223 254 164 363 S401 373 163 255 336 287 384 191 234 158 347 S401 381 163 267 348 287 416 239 250 170 495 S441 333 163 279 300 287 408 255 238 146 351 S441 333 163 283 316 287 416 255 246 146 387 S451 345 163 -9 -9 287 412 187 234 158 -9 S451 389 163 -9 -9 287 416 215 234 176 -9 S491 349 163 271 312 285 388 191 238 164 351 S491 357 163 283 344 287 400 191 238 182 495 S521 353 163 259 300 287 380 195 258 -9 335 S521 385 163 279 320 287 424 195 258 -9 351 S561 325 163 267 300 287 400 263 226 146 355 S561 389 163 283 316 287 400 263 238 188 407 S571 357 163 263 300 287 412 199 234 146 343 S571 389 163 271 316 287 412 199 242 176 391 S611 369 163 287 316 287 380 175 246 164 339 S611 393 163 287 316 287 432 191 250 182 347 S621 377 163 275 316 287 392 -9 254 164 363 S621 385 165 275 348 287 392 -9 254 164 363 S631 333 163 267 332 285 388 175 234 158 347 S631 401 163 271 340 287 408 199 234 164 355 S671 377 163 251 320 287 412 187 246 158 351 S671 377 163 275 352 287 416 187 246 158 379 S681 377 163 287 332 287 404 291 238 164 339 S681 405 163 287 348 287 420 291 246 176 591 S691 325 163 255 308 285 364 191 226 -9 343 S691 353 163 279 316 285 416 307 242 -9 359 S711 385 163 255 328 285 384 251 230 164 343 S711 385 165 267 352 287 424 251 230 168 363 S731 341 159 259 308 285 416 191 230 164 335 S731 369 163 259 324 285 416 191 230 176 467 S741 -9 -9 271 320 287 384 183 226 152 347 S741 -9 -9 271 348 287 408 227 226 170 455 S751 373 163 251 308 285 380 195 234 158 351 S751 385 163 279 348 287 400 195 234 164 591 S791 333 163 275 304 287 400 211 218 170 339 S791 369 163 279 348 287 400 267 234 176 339 S801 393 163 259 316 287 388 255 230 170 407 S801 401 163 267 332 287 388 255 266 176 587 S
Re: [R] Export to .txt
Look at the txtStart function in the TeachingDemos package. It works like sink but also includes commands as well as output. Though I have never tried it with browser() (and it does not always include the results of errors). Another option in to use some type of editor that links with R such as emacs/ESS or tinn-R (or other) and then save the entire transcript. -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of stat.kk Sent: Tuesday, November 01, 2011 4:15 PM To: r-help@r-project.org Subject: [R] Export to .txt Hi, I would like to export all my workspace (even with the evaluation of commands) to the text file. I know about the sink() function but it doesnt work as I would like. My R-function looks like this: there are instructions for user displayed by cat() command and browser() commands for fulfilling them. While using the sink() command the instructions dont display :( Can anyone help me with a equivalent command to File - Save to file... option? Thank you very much. -- View this message in context: http://r.789695.n4.nabble.com/Export-to-txt-tp3965699p3965699.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] List of user installed packages
Running rownames(installed.packages()) will tell you the names of all packages of the version of R in which you are running the command. http://cran.r-project.org/doc/FAQ/R-FAQ.html#R-Add_002dOn-Packages tells you the names of the packages which were installed with R itself. On Nov 5, 2011, at 2:37 PM, Cem Girit wrote: > Hello, > > > >I am going to install the new version of R 2.14.1. After the > installation, I want to copy my installed packages to the new library. But > since over time I forgot which ones I installed I want to get a list of all > the packages I installed among the packages installed initially by the > R-installer. Is this possible? > > > > Cem > > > > Cem Girit, PhD > > > > Biopticon Corporation > > 182 Nassau Street, Suite 204 > > Princeton, NJ 08542 > > Tel: (609)-853-0231 > > Email:gi...@biopticon.com > > > > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. signature.asc Description: Message signed with OpenPGP using GPGMail __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nested "for" loops
If in fact this is homework, you will do yourself, your classmates, and possibly your teacher if you let them know that, at least in R, almost anything you can do in a for() loop can be done easier and faster with vectorization. If you teacher can't comprehend this, get him fired. a<-c(4,6,3) b<- c( 9,4,1) d <- c(4,7,2,8) winning.value <- max(outer(a,outer(b,d,"*"),"*")) From: R. Michael Weylandt Date: Sat, 05 Nov 2011 10:21:05 -0400 Why do you need to do it with nested for loops? It is of course possible - and I hinted how to do it in my first email - but there's no reason as far as I can see to do so, particularly as a means of MLE. Sounds suspiciously like homework... Michael On Nov 4, 2011, at 10:14 PM, nick_pan wrote: > Thank you , this works but I have to do it with nested for loops... > > Could you suggest me a way ? > -- Sent from my Cray XK6 "Pendeo-navem mei anguillae plena est." __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] linear against nonlinear alternatives - quantile regression
I suppose this constitutes thread drift, but your simple example, Frank, made wonder if Rq() accepts a vector argument for tau. I seem to remember that Koencker's rq() does.. Normally I would consult the help page, but the power is still out here in Central Connecticut and I am corresponding with a less capable device. I am guessing that if Rq() does accept such a vector that the form of the nonlinearity would be imposed at all levels of tau. -- David On Nov 5, 2011, at 10:43 AM, Frank Harrell wrote: > Just to address a piece of this - in the case in which you are currently > focusing on only one quantile, the rms package can help by fitting > restricted cubic splines for covariate effects, and then run anova to test > for nonlinearity (sometimes a dubious practice because if you then remove > nonlinear terms you are mildly cheating). > > require(rms) > f <- Rq(y ~ x1 + rcs(x2,4), tau=.25) > anova(f) # tests associations and nonlinearity of x2 > > Frank > > Julia Lira wrote: >> >> Dear all, >> >> I would like to know whether any specification test for linear against >> nonlinear model hypothesis has been implemented in R using the quantreg >> package. >> >> I could read papers concerning this issue, but they haven't been >> implemented at R. As far as I know, we only have two specification tests >> in this line: anova.rq and Khmaladze.test. The first one test equality and >> significance of the slopes across quantiles and the latter one test if the >> linear specification is model of location or location and scale shift. >> >> Do you have any suggestion? >> >> Thanks a lot! >> >> Best regards, >> >> Julia >> >>[[alternative HTML version deleted]] >> >> __ >> R-help@ 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. >> > > > - > Frank Harrell > Department of Biostatistics, Vanderbilt University > -- > View this message in context: > http://r.789695.n4.nabble.com/linear-against-nonlinear-alternatives-quantile-regression-tp3993327p3993416.html > Sent from the R help mailing list archive at Nabble.com. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] unsuscribe
Although it is possible to communicate with the list server via email, most people have better luck with the web interface. That applies with even greater force when the person is unable to spell the server commands. https://stat.ethz.ch/mailman/options/r-help -- David. On Nov 5, 2011, at 10:40 AM, Jimmy Barrera wrote: > > > ^_^ > > >> From: dwinsem...@comcast.net >> Date: Fri, 4 Nov 2011 22:02:13 -0400 >> To: peter.langfel...@gmail.com >> CC: r-help@r-project.org; christian.langk...@gmxpro.de >> Subject: Re: [R] 12th Root of a Square (Transition) Matrix >> >> This is just one of many 12-th roots. (Peter knows this i'm sure.) The >> negative of this would also be an nth root, and I read that there are quite >> few others that arise from solutions based on permuting negatives of eigen >> values of a triangularized form. But as I said , I'm not a matrix mechanic, >> so no code for that. >> >> -- >> David. >> >> On Nov 4, 2011, at 6:10 PM, Peter Langfelder >> wrote: >> >>> On Fri, Nov 4, 2011 at 2:37 PM, David Winsemius >>> wrote: >>> The 12th (matrix) root of M: e^( 1/n * log(M) ) > require(Matrix) > M1.12 <- expm( (1/12)*logm(M) ) >>> >>> I like this - haven't thought of the matrix algebra functions in Matrix. >>> >>> Thanks, >>> >>> Peter >> >> [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] 3-D ellipsoid equations update2. Error message when I run R code.
+ Hello, I want to delete prior questions online but am getting an error message? Please see R code in enclosed file. I don't understand the error message. The parametric equations of an ellipsoid can be written in terms of spherical coordinates. The three spherical coordinates are converted to Cartesian coordinates by X=a cos (α) sin(θ) Y=b sin(α) sin(θ) Z=c cos(θ) The parameter α varies from 0 to 2 π and θ varies from 0 to π . Here ( X o , Y o ,Z o ) is the center of the ellipsoid, and θ is the angle of rotation. I need to come up with an expression for the ellipsoid expressed parametrically as the path of a point in 3- space. I think that it is something like the following: x (alpha)<- x0 + a * cos(theta) * cos(alpha) - b * sin(theta) * sin(alpha) y(alpha) <- y0 + a * cos(theta) * sin(alpha) + b * sin(theta) * cos(alpha) z (alpha)<- z0 + a * cos(theta) * sin(alpha) + c * sin(theta) * cos(alpha) Do I have these equations correct? Most of the books I have read use eigenvectors. The eigenvectors of course consist of the direction cosines. My difficulty is going from that approach to the approach that Alberto Monteiro took in his message on the 9 October 2006. I understand the R code and am using it for a two-dimensional ellipse problem. There does not seem to be allowance for the new coordinates of the center of the ellipsoid under the transformation when using direction cosines. By that I mean adding the centroid coordinates would not be necessary. I need to come up with an example where I do it both ways(as above and using direction cosines). My confusion lies in the fact that rather than one rotational angle theta there are 9 direction cosines. Can you assist with this. Sincerely, Mary A. Marion mat <- matrix(c(1,1,1,1,3,2,1,2,2), 3, 3) eigens <- eigen(mat) evs <- eigens$values evecs <- eigens$vectors a <- evs[1] b <- evs[2] c <- evs[3] x0 <- 5 y0 <- 10 z0 <- 3 alpha <- atan(evecs[ , 1][2] / evecs[ , 1][1]) theta <- seq(0, 2 * pi, length=(1000)) # See Timm, Multivariate Analysis, page 62, Figure 1.6.2 # See Lay, Linear Algebra, third edition, page 84 x <- x0 + a * cos(theta) * cos(alpha) - b * sin(theta) * sin(alpha) y <- y0 + a * cos(theta) * sin(alpha) + b * sin(theta) * cos(alpha) z <- z0 + a * cos(theta) * sin(alpha) + c * sin(theta) * cos(alpha) plot(x, y, z, type = "l", main = expression("x y z"), asp = 1, xlim=c(-5,20),ylim=c(-5,20), zlim=-5,20,xaxs = 'i', yaxs = 'i', zaxs='i') error: Error in strsplit(log, NULL) : non-character argument __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] List of user installed packages
I think the installed.packages() function can give you what you need, specifically look at the priority argument. Also check this out http://stackoverflow.com/questions/1401904/painless-way-to-install-a-new-version-of-r Michael On Sat, Nov 5, 2011 at 9:37 AM, Cem Girit wrote: > Hello, > > > > I am going to install the new version of R 2.14.1. After the > installation, I want to copy my installed packages to the new library. But > since over time I forgot which ones I installed I want to get a list of all > the packages I installed among the packages installed initially by the > R-installer. Is this possible? > > > > Cem > > > > Cem Girit, PhD > > > > Biopticon Corporation > > 182 Nassau Street, Suite 204 > > Princeton, NJ 08542 > > Tel: (609)-853-0231 > > Email:gi...@biopticon.com > > > > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] similar package in R like "SKEW CALCULATOR"?
David Winsemius comcast.net> writes: > > From that hand waving description it would be difficult to tell. Sounds like a reinvention of the Pareto > Index, for which you can find many packages that provide facilities: > > http://finzi.psych.upenn.edu/cgi-bin/namazu.cgi?query=Pareto+index&max=100&; > result=normal&sort=score&idxname=functions&idxname=vignettes&idxname=views > I think the author is looking for specific measures of "reproductive skew", a term from behavioral ecology/evolutionary biology. (PS for first-time questioners: you should not assume that general readers of the R-help list know much about your particular subject area. Short definitions and web references are helpful.) Based on a quick RSiteSearch("{reproductive skew}") library(sos) findFn('{reproductive skew}') googling "reproductive skew CRAN" searching for "reproductive skew" at http://rseek.org I don't think so ... There isn't an R list for evolutionary or behavioral biology, as far as I know, but you might try asking this question on the r-sig-ecology mailing list. I suspect it would not be terribly hard to implement these methods in R, but I can't find any evidence that anyone has done it and made it publicly available. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] List of user installed packages
Hello, I am going to install the new version of R 2.14.1. After the installation, I want to copy my installed packages to the new library. But since over time I forgot which ones I installed I want to get a list of all the packages I installed among the packages installed initially by the R-installer. Is this possible? Cem Cem Girit, PhD Biopticon Corporation 182 Nassau Street, Suite 204 Princeton, NJ 08542 Tel: (609)-853-0231 Email:gi...@biopticon.com [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] linear against nonlinear alternatives - quantile regression
Just to address a piece of this - in the case in which you are currently focusing on only one quantile, the rms package can help by fitting restricted cubic splines for covariate effects, and then run anova to test for nonlinearity (sometimes a dubious practice because if you then remove nonlinear terms you are mildly cheating). require(rms) f <- Rq(y ~ x1 + rcs(x2,4), tau=.25) anova(f) # tests associations and nonlinearity of x2 Frank Julia Lira wrote: > > Dear all, > > I would like to know whether any specification test for linear against > nonlinear model hypothesis has been implemented in R using the quantreg > package. > > I could read papers concerning this issue, but they haven't been > implemented at R. As far as I know, we only have two specification tests > in this line: anova.rq and Khmaladze.test. The first one test equality and > significance of the slopes across quantiles and the latter one test if the > linear specification is model of location or location and scale shift. > > Do you have any suggestion? > > Thanks a lot! > > Best regards, > > Julia > > [[alternative HTML version deleted]] > > __ > R-help@ 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. > - Frank Harrell Department of Biostatistics, Vanderbilt University -- View this message in context: http://r.789695.n4.nabble.com/linear-against-nonlinear-alternatives-quantile-regression-tp3993327p3993416.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] unsuscribe
^_^ > From: dwinsem...@comcast.net > Date: Fri, 4 Nov 2011 22:02:13 -0400 > To: peter.langfel...@gmail.com > CC: r-help@r-project.org; christian.langk...@gmxpro.de > Subject: Re: [R] 12th Root of a Square (Transition) Matrix > > This is just one of many 12-th roots. (Peter knows this i'm sure.) The > negative of this would also be an nth root, and I read that there are quite > few others that arise from solutions based on permuting negatives of eigen > values of a triangularized form. But as I said , I'm not a matrix mechanic, > so no code for that. > > -- > David. > > On Nov 4, 2011, at 6:10 PM, Peter Langfelder > wrote: > > > On Fri, Nov 4, 2011 at 2:37 PM, David Winsemius > > wrote: > > > >> > >> The 12th (matrix) root of M: e^( 1/n * log(M) ) > >> > >>> require(Matrix) > >>> M1.12 <- expm( (1/12)*logm(M) ) > > > > I like this - haven't thought of the matrix algebra functions in Matrix. > > > > Thanks, > > > > Peter > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] zoo performance regression noticed (1.6-5 is faster...)
On Fri, Nov 4, 2011 at 1:02 PM, Gabor Grothendieck wrote: > On Fri, Nov 4, 2011 at 12:56 PM, Gabor Grothendieck > wrote: >> On Fri, Nov 4, 2011 at 12:34 PM, James Marca >> wrote: >>> Good morning, >>> >>> I have discovered what I believe to be a performance regression >>> between Zoo 1.6x and Zoo 1.7-6 in the application of rollapply. >>> On zoo 1.6x, rollapply of my function over my data takes about 20 >>> minutes. Using 1.7-6, the same code takes about 6 hours. >>> >>> R --version >>> R version 2.13.1 (2011-07-08) >>> Copyright (C) 2011 The R Foundation for Statistical Computing >>> ISBN 3-900051-07-0 >>> Platform: x86_64-pc-linux-gnu (64-bit) >>> >>> Two versions of zoo 1.6 run *fast* On one machine I am running >>> >>> less /usr/lib64/R/library/zoo/DESCRIPTION >>> Package: zoo >>> Version: 1.6-3 >>> Date: 2010-04-23 >>> Title: Z's ordered observations >>> ... >>> Packaged: 2010-04-23 07:28:47 UTC; zeileis >>> Repository: CRAN >>> Date/Publication: 2010-04-23 07:43:54 >>> Built: R 2.10.1; ; 2010-04-25 06:41:34 UTC; unix >>> >>> (Thankfully I forgot to upgrade.packages() on this machine!) >>> >>> On the other >>> >>> Package: zoo >>> Version: 1.6-5 >>> Date: 2011-04-08 >>> ... >>> Packaged: 2011-04-08 17:13:47 UTC; zeileis >>> Repository: CRAN >>> Date/Publication: 2011-04-08 17:27:47 >>> Built: R 2.13.1; ; 2011-11-04 15:49:54 UTC; unix >>> >>> I have stripped out zoo 1.7-6 from all my machines. >>> >>> I tried to ensure all libraries were identical on the two machines >>> (using lsof), and after finally downgrading zoo I got the second >>> machine to be as fast as the first, so I am quite certain the >>> difference in speed is down to the Zoo version used. >>> >>> My code runs a fairly simple function over a time series using the >>> following call to process a year of 30s data (9 columns, about a >>> million rows): >>> >>> vals <- rollapply(data=ts.data[,c(n.3.cols, o.3.cols,volocc.cols)] >>> ,width=40 >>> >>> ,FUN=rolling.function.fn(n.cols=n.3.cols,o.cols=o.3.cols,vo.cols=volocc.cols) >>> ,by.column=FALSE >>> ,align='right') >>> >>> >>> (The rolling.function.fn call returns a function that is initialized >>> with the initial call above (a trick I learned from Javascript)) >>> >>> If this is a known situation with the new 1.7 generation Zoo, my >>> apologies and I'll go away. If my code could be turned into a useful >>> test, I'd be happy to help out as much as I'm able. Given the extreme >>> runtime difference though, I thought I should offer my help in this >>> case, since zoo is such a useful package in my work. >> >> This was a known problem and was fixed but if its still there then >> there must be some other condition under which it can occur as well. >> If you can provide a small self contained reproducible example it >> would help in tracking it down. >> >> -- >> Statistics & Software Consulting >> GKX Group, GKX Associates Inc. >> tel: 1-877-GKX-GROUP >> email: ggrothendieck at gmail.com >> > > Also, as a workaround you can try this to use an old rollapply in a > new version of zoo: > > library(zoo) > source("http://r-forge.r-project.org/scm/viewvc.php/*checkout*/pkg/zoo/R/rollapply.R?revision=817&root=zoo";) > rollapply(...whatever...) > Have looked at it and there is now a performance improvement in the development version of rollapply that gives an order of magnitude performance boost in the following case: > library(zoo) > n <- 1 > z <- zoo(cbind(a = 1:n, b = 1:n)) > system.time(rollapplyr(z, 10, sum, by.column = FALSE)) user system elapsed 8.800.028.97 > > # download rollapply rev 913 from svn repo and rerun > source("http://r-forge.r-project.org/scm/viewvc.php/*checkout*/pkg/zoo/R/rollapply.R?revision=913&root=zoo";) > system.time(rollapplyr(z, 10, sum, by.column = FALSE)) user system elapsed 0.520.020.53 -- Statistics & Software Consulting GKX Group, GKX Associates Inc. tel: 1-877-GKX-GROUP email: ggrothendieck at gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nested "for" loops
Why do you need to do it with nested for loops? It is of course possible - and I hinted how to do it in my first email - but there's no reason as far as I can see to do so, particularly as a means of MLE. Sounds suspiciously like homework... Michael On Nov 4, 2011, at 10:14 PM, nick_pan wrote: > Thank you , this works but I have to do it with nested for loops... > > Could you suggest me a way ? > > > > > > -- > View this message in context: > http://r.789695.n4.nabble.com/nested-for-loops-tp3992089p3992324.html > Sent from the R help mailing list archive at Nabble.com. > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Error in eigen(a$hessian) : infinite or missing values in 'x'
Dear R-users, I'm estimating a two- dimensional state-space model using the FKF package. The resulting log likelihood function is maximized using auglag from the Alabama package. The procedure works well for a subset of my data, but if I try to use the entire data set I get the following error message. Error in eigen(a$hessian) : infinite or missing values in 'x' What's even more confusing is that if I estimate the model for a sample say data[1:200,] then there's convergence. If I estimate it for data[1:300, ] then I get the error message. But if I estimate the model for data[201,300] it once again converges. Can anyone please enlighten me; where does this error stem from and what can I do about it? Thank you in advance. Kristian [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] linear against nonlinear alternatives - quantile regression
Dear all, I would like to know whether any specification test for linear against nonlinear model hypothesis has been implemented in R using the quantreg package. I could read papers concerning this issue, but they haven't been implemented at R. As far as I know, we only have two specification tests in this line: anova.rq and Khmaladze.test. The first one test equality and significance of the slopes across quantiles and the latter one test if the linear specification is model of location or location and scale shift. Do you have any suggestion? Thanks a lot! Best regards, Julia [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] similar package in R like "SKEW CALCULATOR"?
>From that hand waving description it would be difficult to tell. Sounds like a >reinvention of the Pareto Index, for which you can find many packages that >provide facilities: http://finzi.psych.upenn.edu/cgi-bin/namazu.cgi?query=Pareto+index&max=100&result=normal&sort=score&idxname=functions&idxname=vignettes&idxname=views -- David. On Nov 5, 2011, at 9:07 AM, Knut Krueger wrote: > > Hi to all > is there a similar package like the SKEW CALCULATOR from > Peter Nonacs (University of California - Department of Ecology and > Evolutionary Biology) > > http://www.eeb.ucla.edu/Faculty/Nonacs/shareware.htm > > > Kind Regards Knut > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to write a shapefile with projection
Hi, Sorry i have put such a detailed question to the list about writing a shapefile with projection. I realized that if i use writeOGR from rgdal and not the other write shapefile functions i can get a shapefile with projection recognized by ArcGIS. The command is (in case anybody wonders): ?writeOGR(crest.sp, "I:\\LA_levee\\Shape", "llev_crest_pts6", driver = "ESRI Shapefile") where crest.sp is a spatial point data frame with projection. Thanks, Monica Indeed. writePointsShape() does not write the projection file, but using the function showWKT from rgdal, you can also create one like that: writePointsShape(crest.sp,"crest") cat(showWKT(proj4string(crest.sp)),file="crest.prj") Patrick __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] similar package in R like "SKEW CALCULATOR"?
Hi to all is there a similar package like the SKEW CALCULATOR from Peter Nonacs (University of California - Department of Ecology and Evolutionary Biology) http://www.eeb.ucla.edu/Faculty/Nonacs/shareware.htm Kind Regards Knut __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] 3-D ellipsoid equations
+ Hello, The parametric equations of an ellipsoid can be written in terms of spherical coordinates. The three spherical coordinates are converted to Cartesian coordinates by X=a cos (α) sin(θ) Y=b sin(α) sin(θ) Z=c cos(θ) for α and θ The parameter α varies from 0 to 2 Ï and θ varies from 0 to Ï . Here ( X o , Y o ,Z o ) is the center of the ellipsoid, and θ is the angle of rotation. I need to come up with an expression for the ellipsoid expressed parametrically as the path of a point in 3- space. My first try is that it is something like the following: X(alpha)=Xo+a cos(α) cos( θ )-b sin(α) cos( θ ) + c cos( θ ) Y(alpha)=Yo+ cos (α) sin(θ)+b sin(α) cos (θ) Z(alpha)=Zo+a cos (α) sin(θ) +b sin(α) cos( θ ) Most of the books I have read use eigenvectors. The eigenvectors of course consist of the direction cosines. My difficulty is going from that approach to the approach that Alberto Monteiro took in his message on the 9 October 2006. I understand the R code and am using it for a two-dimensional ellipse problem. There does not seem to be allowance for the new coordinates of the center of the ellipsoid under the transformation when using direction cosines. By that I mean adding the centroid coordinates would not be necessary as is done in my "first try". Can you help me extend this to 3 dimensions? Sincerely, Mary A. Marion + [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] set seed for random draws
I'm suspecting this is confusion about default behavior. R automatically updates the random seed when random numbers are generated (or other random operations are performed). The original poster may have experienced systems where it is up to the user to change the seed. I'd suggest two rules of thumb when coming up against something in R that you aren't sure about: 1. If it is a mundane task, R probably takes care of it. 2. Experiment to see what happens. Of course you could read documentation, but no one does that. On 05/11/2011 00:59, R. Michael Weylandt wrote: This might be more fundamental, but why do you feel the need to reset the seed each loop? There's nothing that suggests you need to... Michael On Fri, Nov 4, 2011 at 8:38 PM, Md Desa, Zairul Nor Deana Binti wrote: Hello, all! I need help on these two problems: 1) If I want to randomly draw numbers from standard normal (or other distributions) in loops e.g.: ty=0; ks=0 for (i in 1:5) { set.seed(14537+i) k<-rnorm(1) ks[i]<-.3*k+.9 if (ty==0) { while ((ks<.2)||(ks>3)) { #set.seed(13237+i*100) k<-rnorm(1) ks[i]-.3*k+.9 } } } } Question: Here I draw initial a, then if the drawn initial a satisfied 2 conditions I redraw a. I set.seed(13237) in the first draw of a, should I set.seed() in the redraw part? 2) I also have more loops after this i loop that also draw from normal(0,1). I want to randomly draws from normal(0,1) for loop j (inside loop j I draw another random numbers from N(0,1)) My question: Should I or shouldn't I set seed again and again for each loop? Why or why not. I guess this problem concerned about setting seed as I want to have different number for each i. Thanks! Deana __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@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. -- Patrick Burns pbu...@pburns.seanet.com twitter: @portfolioprobe http://www.portfolioprobe.com/blog http://www.burns-stat.com (home of 'Some hints for the R beginner' and 'The R Inferno') __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nested "for" loops
Thank you , this works but I have to do it with nested for loops... Could you suggest me a way ? -- View this message in context: http://r.789695.n4.nabble.com/nested-for-loops-tp3992089p3992324.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] HoltWinters in R 2.14.0
You are 100% correct by my estimation however I suppose I am looking for the conditions in the data that might cause the optim() or optimize() functions to fail. I took a brief tour of the HoltWinters source but the code available (readable) online seemed outdated (by way of conflicting descriptions in versioning.). I'll have another poke around the source - that is unless there is someone out there that can clearly state why optimize() fails within the context of the HoltWinters class v. 2.14.0. On Nov 4, 2011, at 8:38 PM, "Prof Brian Ripley [via R]" wrote: > On Fri, 4 Nov 2011, R. Michael Weylandt wrote: > > > I believe there were some changes to Holt-Winters, specifically in re > > optimization that probably lead to your problem, but you'll have to > > provide more details. See the NEWS file for citations about the > > change. If you put example code/data others may be able to help you -- > > I haven't updated yet so I can't be of much help. > > > > Michael > > > > > > On Fri, Nov 4, 2011 at 2:55 PM, TimothyDalbey <[hidden email]> wrote: > >> Hey All, > >> > >> First time on these forums. Thanks in advance. > >> > >> S... I have a process that was functioning well before the 2.14 > >> update. > >> Now the HoltWinters function is throwing an error whereby I get the > >> following: > >> > >> Error in HoltWinters(sales.ts) : optimization failure > Most likely it was incorrect before. You cannot assume that it was > actually 'functioning well': all the cases where we have seen this > message it was giving incorrect answers before and not detecting them. > And in all those cases the model was a bad fit and using starting > values for the optimization helped. > > >> I've been looking around to determine why this happens (see if I can test > >> the data beforehand) but I haven't come across anything. > >> > >> Any help appreciated! > > -- > Brian D. Ripley, [hidden email] > 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 > __ > [hidden email] 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. > > > If you reply to this email, your message will be added to the discussion > below: > http://r.789695.n4.nabble.com/HoltWinters-in-R-2-14-0-tp3991247p3992395.html > To unsubscribe from HoltWinters in R 2.14.0, click here. -- View this message in context: http://r.789695.n4.nabble.com/HoltWinters-in-R-2-14-0-tp3991247p3992497.html Sent from the R help mailing list archive at Nabble.com. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.