[R] poLCA problem
Good morning everyone, I am running into errors with poLCA as follows: Error in round(mf) : non-numeric argument to mathematical function. Here is what I have. Both variables are coded as 1, 2, 3 df <- as.data.frame(data) items <- c("x1", "x2") <- there are more variables but shortened for this purpose df2 <- df[items] i <- cbind(x1, x2)~1 poLCA (i, df2, nclass=2, maxiter=100, nrep=10, verbose =TRUE) It is after the poLCA that I get the error. Any thoughts on what is causing this? Thanks, Scott [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] cooks distance for repeated measures anova
Hi all, Is there a way to get cook's distance for a repeated measures anova? Neither cooks.distance or CookD from the predictmeans package seem to allow for this. For example, if I have the model data(iris) mod<-aov(Sepal.Length ~ Petal.Length + Petal.Width + Error(Species), data=iris) both cooks.distance(mod) and library(predictmeans) CookD(mod, group=Species) give an error saying they don't support an aovlist object. I would prefer a method to get a cook's distance for each category in my repeated factor (i.e. Species), rather than each observation. Thanks! -- Walker Pedersen, Ph.D. Center for Healthy Minds University of Wisconsin -- Madison [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] Quantitative Methods Workshops in May 2020
FORWARDED – The following message has been forwarded and is not related to the University of Guelph. Apologies for the cross-posting. Good morning everyone. We sincerely hope you are all keeping safe and healthy while we all endure this pandemic. As a result of the restrictions on public gatherings, we are going to be delivering our May 2020 workshops via online live streaming now as opposed to face-to-face at McMaster University. We will be providing materials ahead of time for all participants and will still include time for individual consultation regarding your own data (which can be scheduled at the live-stream workshop). Each workshop provides a hands-on opportunity to learn using both R (with R-Studio) and Mplus. For further information and to register, please go to: https://workshops.enablytics.com/ [1] Introduction to Structural Equation Modeling This one-day hands-on workshop covers various introductory topics in structural equation modeling with continuous and categorical variables. Topics include, assumptions and data considerations, model creation, identification, and evaluation, multiple regression vs path analysis, path analysis, testing direct and indirect effects, and confirmatory factor analysis. Examples will be demonstrated in both R (using R-Studio) and Mplus. Syntax and output for both programs will be provided for all examples covered in the workshop. [2] Advanced Structural Equation Modeling This one-day hands-on workshop covers various advanced topics in structural equation modeling with continuous and categorical variables. Topics include, model creation, identification, and evaluation, testing moderation, mediation and moderated mediation, multiple group modeling, handling missing data, measurement invariance and power analysis. Examples will be demonstrated in both R (using R-Studio) and Mplus. Syntax and output for both programs will be provided for all examples covered in the workshop. [3] Growth Modeling This one-day hands-on workshop covers various topics in growth modeling (longitudinal modeling) with continuous and categorical variables. Topics include, growth modeling without covariates, growth modeling with time invariant and varying covariates, centering points, piecewise growth modeling, missing data and power analysis. Examples will be demonstrated in both R (using R-Studio) and Mplus. Syntax and output for both programs will be provided for all examples covered in the workshop. [4] Multilevel Modeling This two day hands-on workshop covers various topics in multilevel modeling with continuous and categorical variables. Topics include, when multilevel analysis is necessary, multilevel regression, random slopes and cross-level effects, multilevel confirmatory factor analysis and the MIMIC model, multilevel path analysis, multilevel mediation and moderation, multilevel latent variable modeling, longitudinal data, and power analysis. Examples will be demonstrated in both R (using R-Studio) and Mplus. Syntax and output for both programs will be provided for all examples covered in the workshop. Thank you, Scott [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] choose(n, k) as n approaches k
This struck me as incorrect: > choose(3.99, 4) [1] 0.979 > choose(3.999, 4) [1] 0 > choose(4, 4) [1] 1 > choose(4.001, 4) [1] 4 > choose(4.01, 4) [1] 1.02 Should base::choose(n, k) check whether n is within machine precision of k and return 1? Thanks, Erik *** sessionInfo() R version 3.6.0 beta (2019-04-15 r76395) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS High Sierra 10.13.6 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Creating hanging bar plot in r from dplyr
That is perfect. Thanks! -- Scott R. Colwell, PhD On 2019-04-20, 1:23 PM, "Jeff Newmiller" wrote: Not really sure I understand what you want. Here is some code to consider: library(ggplot2) library(dplyr) library(tidyr) dta <- read.table( text = "samp.N RSQMRB_uc MRB_sb MRB_bp 50 0.3 1.4237.6 37.6 50 0.4 8.6143.1 43.1 50 0.5 7.4131.6 31.6 50 0.6 5.0621.5 21.5 50 0.7 3.3814.1 14.1 50 0.8 -1.075.16 5.16 100 0.3 -6.4140.3 40.3 100 0.4 -10.621.0 21.0 100 0.5 -9.0213.2 13.2 100 0.6 -9.855.14 5.14 100 0.7 -7.942.08 2.08 100 0.8 -4.811.28 1.28 ", header = TRUE ) dta2 <- ( dta %>% mutate( samp.N = factor( samp.N ) , RSQ = factor( RSQ ) ) %>% gather( Measure, value, -c( samp.N, RSQ ) ) ) ggplot( dta2, aes( x = RSQ, y = value, fill = samp.N ) ) + geom_bar( stat = "identity", position = "dodge", colour = "black" ) + facet_wrap( ~ Measure, ncol = 1, scale = "free_y" ) + ylab( "" ) On Sat, 20 Apr 2019, Scott Colwell wrote: > I am trying to figure out how to create a hanging bar plot from dplyr. > I have used dplyr as follows: > table4 <- cr %>% > group_by(samp.N, RSQ) %>% > summarize( >MRB_uc = mean(CF.F1F2/0.40*100)-100, >MRB_sb = mean(SBC.F1F2.Alpha/0.40*100) - 100, >MRB_bp = mean(BPC.F1F2.Alpha/0.40*100) - 100 > ) > which provides me with this: > samp.N RSQ MRB_uc MRB_sb MRB_bp > > 1 50 0.3 1.42 37.6 37.6 > 2 50 0.4 8.61 43.1 43.1 > 3 50 0.5 7.41 31.6 31.6 > 4 50 0.6 5.06 21.5 21.5 > 5 50 0.7 3.38 14.1 14.1 > 6 50 0.8 -1.07 5.16 5.16 > 7100 0.3 -6.41 40.3 40.3 > 8100 0.4 -10.6 21.0 21.0 > 9100 0.5 -9.02 13.2 13.2 > 10100 0.6 -9.85 5.14 5.14 > 11100 0.7 -7.94 2.08 2.08 > 12100 0.8 -4.81 1.28 1.28 > What I want to do is create a hanging bar plot with the x-axis being samp.N value by RSQ value. The bars are then values of MRB_uc, MRB_sb, and MRB_bp. Given some values are negative, some bars will be above zero and others below (hence the hanging bar plot) > I don't have any code yet as I am completely unfamiliar with how to do this. Any suggestions would be really appreciated. > Thank you! > Scott > > > > > -- > Scott R. Colwell, PhD > > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > 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. > --- Jeff NewmillerThe . . Go Live... DCN:Basics: ##.#. ##.#. Live Go... Live: OO#.. Dead: OO#.. Playing Research Engineer (Solar/BatteriesO.O#. #.O#. with /Software/Embedded Controllers) .OO#. .OO#. rocks...1k --- __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Creating hanging bar plot in r from dplyr
I am trying to figure out how to create a hanging bar plot from dplyr. I have used dplyr as follows: table4 <- cr %>% group_by(samp.N, RSQ) %>% summarize( MRB_uc = mean(CF.F1F2/0.40*100)-100, MRB_sb = mean(SBC.F1F2.Alpha/0.40*100) - 100, MRB_bp = mean(BPC.F1F2.Alpha/0.40*100) - 100 ) which provides me with this: samp.N RSQ MRB_uc MRB_sb MRB_bp 1 50 0.3 1.42 37.6 37.6 2 50 0.4 8.61 43.1 43.1 3 50 0.5 7.41 31.6 31.6 4 50 0.6 5.06 21.5 21.5 5 50 0.7 3.38 14.1 14.1 6 50 0.8 -1.07 5.16 5.16 7100 0.3 -6.41 40.3 40.3 8100 0.4 -10.6 21.0 21.0 9100 0.5 -9.02 13.2 13.2 10100 0.6 -9.85 5.14 5.14 11100 0.7 -7.94 2.08 2.08 12100 0.8 -4.81 1.28 1.28 What I want to do is create a hanging bar plot with the x-axis being samp.N value by RSQ value. The bars are then values of MRB_uc, MRB_sb, and MRB_bp. Given some values are negative, some bars will be above zero and others below (hence the hanging bar plot) I don't have any code yet as I am completely unfamiliar with how to do this. Any suggestions would be really appreciated. Thank you! Scott -- Scott R. Colwell, PhD [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] 2019 Spring Workshops using R
Apologies for any cross-postings. Just a reminder that registration is open for four quantitative methods workshops in May of 2019. Each workshop features hands-on examples in Mplus and R, plus lots of opportunities to discuss the analysis for your own research. For more information and to register, please see https://enablytics.com/workshop_events/ [1] Introductory Structural Equation Modeling - May 5, 2019 141 Adelaide Street West Toronto, Ontario M5H 3L5 This one-day hands-on workshop covers various introductory topics in structural equation modeling with continuous and categorical variables. Topics include, assumptions and data considerations, model creation, identification, and evaluation, multiple regression vs path analysis, path analysis, testing direct and indirect effects, and confirmatory factor analysis. Syntax and output for both Mplus and R will be provided for all examples covered in the workshop. [2] Advanced Structural Equation Modeling - May 6, 2019 141 Adelaide Street West Toronto, Ontario M5H 3L5 This one-day hands-on workshop covers various advanced topics in structural equation modeling with continuous and categorical variables. Topics include, model creation, identification, and evaluation, testing moderation, mediation and moderated mediation, multiple group modeling, handling missing and messy data, measurement invariance and power analysis. Syntax and output for both Mplus and R will be provided for all examples covered in the workshop. [3] Multilevel Modeling - May 7 and 8, 2019 141 Adelaide Street West Toronto, Ontario M5H 3L5 This two day hands-on workshop covers various topics in multilevel modeling with continuous and categorical variables. Topics include, when multilevel analysis is necessary, multilevel regression, random slopes and cross-level effects, multilevel confirmatory factor analysis and the MIMIC model, multilevel path analysis, multilevel mediation and moderation, multilevel latent variable modeling, longitudinal data, and power analysis. Syntax and output for both Mplus and R will be provided for all examples covered in the workshop. [4] Growth Modeling - May 9 and 10, 2019 141 Adelaide Street West Toronto, Ontario M5H 3L5 This two-day hands-on workshop covers various topics in growth modeling (longitudinal modeling) with continuous and categorical variables. Topics include, growth modeling without covariates, growth modeling with time invariant and varying covariates, centering points, piecewise growth modeling, autoregressive latent trajectory modeling (ALT -- Scott R. Colwell, PhD, CStat, PStat [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] plot one levelplot over another
Hi, I am using levelplot() to plot a primary response surface, z1. I wish to write a custom panel function that will let me plot another surface z2 = f(x,y) over z1. The second surface z2 is either NA or 1, and at locations where z2 = 1, I will use a color with low alpha to let the the z1 surface show through. How can I do this? Regards, Scott Waichler Pacific Northwest National Laboratory Richland, WA USA __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] GLM Model Summary
R-Help: We are working with your GLM R package. The Summary(Model) now gets printed by the program as one object and we want to put the coefficient columns into Excel. We took an initial stab at this by counting the number of characters occupied by each column. But we have now learned that the number of characters in a column depends on the length of the variable names, so is not a constant number (e.g., 54 characters to a line). We therefore ask, is it possible for us to get the Summary(Model) column by column, i.e., a separate object for each column? That way we could assemble an Excel table easily rather than having to count the number of characters. Is this possible for us to do by ourselves? Or could you modify the package in some way? We appreciate your attention. Thank you! Scott Neslin Prasad Vana Dartmouth College [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] seq() problem with chron
Hi, I encountered the problem below where the last value in the chron vector created with seq() should have a time of 15:30, but instead has 15:15. What causes this and how can I make sure that the last value in the chron vector is the same as the "to" value in seq()? library(chron) dt1 <- chron("02/20/13", "00:00:00") dt2 <- chron("07/03/18", "15:30:00") dt <- seq(from=dt1, to=dt2, by=1/(24*4)) dt[length(dt)] #[1] (07/03/18 15:15:00) Thanks, Scott Waichler Pacific Northwest National Laboratory scott.waich...@pnnl.gov __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] download.file() problems with binary files containing EOF byte in Windows
Hello, I'm trying to get a package to pass win-builder and have been having a bit of trouble with Windows R and binary files (in my case a small .tar.gz used in testing). After a little debugging, I think I've narrowed it down to download.file() truncating files to the first '1a' byte (often used for EOF but I think a valid byte inside gzip files) on downloads from local "file://xxx". I'm trying to figure out if this is a known "feature" of Windows that I should just avoid or does this seem like a bug? For example: #write a file starting with byte 1a (decimal 26) writeBin(26:100,'tmp.bin',size=1) download.file('file://tmp.bin','download.bin') file.size('tmp.bin') file.size('download.bin') On Windows (session info below), I get file sizes of 75 and 0 and on Linux I get 75 and 75. As a more real world example, if I download.file() on a .gz file then a remote download seems to return different size files from a local download. For example for a gz file from a google hit about gzip (http://commandlinefanatic.com/cgi-bin/showarticle.cgi?article=art053): download.file('http://commandlinefanatic.com/gunzip.c.gz','gunzip.c.gz') download.file('file://gunzip.c.gz','dl.gz') file.size('gunzip.c.gz') file.size('dl.gz') I get a 4704 byte file for the remote download and 360 for the local download in Windows (versus 4704 and 4704 on Linux). Note that the 361st byte is 1a: readBin('gunzip.c.gz','raw',361) The various download.file options don't seem to fix this with the same 360 bytes for: download.file('file://gunzip.c.gz','dl.gz',mode='wb') file.size('dl.gz') download.file('file://gunzip.c.gz','dl.gz',mode='wb',method='internal') file.size('dl.gz') It looks like the 'auto' and 'internal' methods both resolve to the 'wininet' method on Windows and mode is automatically set to 'wb' for gz files so maybe not surprising those don't change things. Thanks, Scott ## Windows sessionInfo(): R version 3.5.1 (2018-07-02) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 8.1 x64 (build 9600) Matrix products: default locale: [1] LC_COLLATE=English_United States.1252 [2] LC_CTYPE=English_United States.1252 [3] LC_MONETARY=English_United States.1252 [4] LC_NUMERIC=C [5] LC_TIME=English_United States.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base loaded via a namespace (and not attached): [1] compiler_3.5.1 ## Linux sessionInfo(): R version 3.4.4 (2018-03-15) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 16.04.5 LTS Matrix products: default BLAS: /usr/lib/libblas/libblas.so.3.6.0 LAPACK: /usr/lib/lapack/liblapack.so.3.6.0 locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base loaded via a namespace (and not attached): [1] compiler_3.4.4 __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] [Rd] source(echo = TRUE) with a iso-8859-1 encoded file gives an error
On Fri, May 04, 2018 at 10:58:26PM +, Ista Zahn wrote: > On Fri, May 4, 2018 at 4:47 PM, Scott Kostyshak <skostys...@ufl.edu> wrote: > > I have very little knowledge about file encodings and would like to > > learn more. > > > > I've read the following pages to learn more: > > > > > > https://urldefense.proofpoint.com/v2/url?u=http-3A__stat.ethz.ch_R-2Dmanual_R-2Ddevel_library_base_html_Encoding.html=DwIFaQ=pZJPUDQ3SB9JplYbifm4nt2lEVG5pWx2KikqINpWlZM=neJ42wVqpDzuvOKMBML6-HnbH0l0aXpb0ZUFWoGb-Bo=yaDPpePO4lxR7-PBircARZlFh-GVyi5sTNtjTr_JZ7U=PSqR5opjnHspAeM6Edm1ddsaY3ok1bnV-t6W4MKtVCM= > > > > https://urldefense.proofpoint.com/v2/url?u=https-3A__stackoverflow.com_questions_4806823_how-2Dto-2Ddetect-2Dthe-2Dright-2Dencoding-2Dfor-2Dread-2Dcsv=DwIFaQ=pZJPUDQ3SB9JplYbifm4nt2lEVG5pWx2KikqINpWlZM=neJ42wVqpDzuvOKMBML6-HnbH0l0aXpb0ZUFWoGb-Bo=yaDPpePO4lxR7-PBircARZlFh-GVyi5sTNtjTr_JZ7U=1M6pNfwFR5uG5DkSAHPpXZKYETCiwV1wsJxpew6lThY= > > > > https://urldefense.proofpoint.com/v2/url?u=https-3A__developer.r-2Dproject.org_Encodings-5Fand-5FR.html=DwIFaQ=pZJPUDQ3SB9JplYbifm4nt2lEVG5pWx2KikqINpWlZM=neJ42wVqpDzuvOKMBML6-HnbH0l0aXpb0ZUFWoGb-Bo=yaDPpePO4lxR7-PBircARZlFh-GVyi5sTNtjTr_JZ7U=hAF57aL9khHQ_2Ndars7qMO-FoqxnnmOiEDIprsllko= > > > > The last one, in particular, has been very helpful. I would be > > interested in any further references that you suggest. > > > > I attach a file that reproduces the issue I would like to learn more > > about. I do not know if the file encoding will be correctly preserved > > through email, so I also provide the file (temporarily) on Dropbox here: > > > > > > https://urldefense.proofpoint.com/v2/url?u=https-3A__www.dropbox.com_s_3lbgebk7b5uaia7_encoding-5Fexport-5Fissue.R-3Fdl-3D0=DwIFaQ=pZJPUDQ3SB9JplYbifm4nt2lEVG5pWx2KikqINpWlZM=neJ42wVqpDzuvOKMBML6-HnbH0l0aXpb0ZUFWoGb-Bo=yaDPpePO4lxR7-PBircARZlFh-GVyi5sTNtjTr_JZ7U=fGtYdB-U7ktXVFeniRudE-ZmxmCP3ZUfeLOvJ0AJwqs= > > > > The file gives an error when using "source()" with the > > argument echo = TRUE: > > > > > source("encoding_export_issue.R", echo = TRUE) > > Error in nchar(dep, "c") : invalid multibyte string, element 1 > > In addition: Warning message: > > In grepl("^[[:blank:]]*$", dep[1L]) : > > input string 1 is invalid in this locale > > > > The problem comes from the "á" character in the .R file. The file > > appears to be encoded as "iso-8859-1": > > > > $ file --mime-encoding encoding_export_issue.R > > encoding_export_issue.R: iso-8859-1 > > > > Note that for me: > > > > > getOption("encoding") > > [1] "native.enc" > > > > so "native.enc" is used for the "encoding" argument of source(). > > > > The following two calls succeed: > > > > > source("encoding_export_issue.R", echo = TRUE, encoding = "unknown") > > > source("encoding_export_issue.R", echo = TRUE, encoding = "iso-8859-1") > > > > Is this file a valid "iso-8859-1" encoded file? > > The one you attached is not. The one linked to in dropbox is. > > Why does source() fail > > in the case of encoding set to "native.enc"? Is it because of the > > settings to UTF-8 in my locale (see info on my system at the bottom of > > this email). > > Yes. > > > > > I'm guessing it would be a bad idea to put > > > > options(encoding = "unknown") > > > > in my .Rprofile, because it is difficult to always correctly guess the > > encoding of files? > > My guess is that the issue is less about the difficulty of guessing > the encoding, and more about the time it takes to do so. That's not > particularly relevant for the "source" function, but the encoding > option is used by many of the file IO functions in R and so has > implications well beyond the behavior of "source". Ah I did not think about this possibility. Makes sense. > > Is there a reason why setting it to "unknown" would > > lead to more problems than leaving it set to "native.enc"? > > It depends on what you are actually doing. If you are on a UTF-8 > locale and working exclusively with UTF-8 files, setting > options(encoding = "unknown") will just slow down your file IO by > checking for the encoding every time. Good to know. Thank you for your response, Ista. Scott -- Scott Kostyshak Assistant Professor of Economics University of Florida https://people.clas.ufl.edu/skostyshak/ > > > > I've reproduced the
[R] [Rd] source(echo = TRUE) with a iso-8859-1 encoded file gives an error
I have very little knowledge about file encodings and would like to learn more. I've read the following pages to learn more: http://stat.ethz.ch/R-manual/R-devel/library/base/html/Encoding.html https://stackoverflow.com/questions/4806823/how-to-detect-the-right-encoding-for-read-csv https://developer.r-project.org/Encodings_and_R.html The last one, in particular, has been very helpful. I would be interested in any further references that you suggest. I attach a file that reproduces the issue I would like to learn more about. I do not know if the file encoding will be correctly preserved through email, so I also provide the file (temporarily) on Dropbox here: https://www.dropbox.com/s/3lbgebk7b5uaia7/encoding_export_issue.R?dl=0 The file gives an error when using "source()" with the argument echo = TRUE: > source("encoding_export_issue.R", echo = TRUE) Error in nchar(dep, "c") : invalid multibyte string, element 1 In addition: Warning message: In grepl("^[[:blank:]]*$", dep[1L]) : input string 1 is invalid in this locale The problem comes from the "á" character in the .R file. The file appears to be encoded as "iso-8859-1": $ file --mime-encoding encoding_export_issue.R encoding_export_issue.R: iso-8859-1 Note that for me: > getOption("encoding") [1] "native.enc" so "native.enc" is used for the "encoding" argument of source(). The following two calls succeed: > source("encoding_export_issue.R", echo = TRUE, encoding = "unknown") > source("encoding_export_issue.R", echo = TRUE, encoding = "iso-8859-1") Is this file a valid "iso-8859-1" encoded file? Why does source() fail in the case of encoding set to "native.enc"? Is it because of the settings to UTF-8 in my locale (see info on my system at the bottom of this email). I'm guessing it would be a bad idea to put options(encoding = "unknown") in my .Rprofile, because it is difficult to always correctly guess the encoding of files? Is there a reason why setting it to "unknown" would lead to more problems than leaving it set to "native.enc"? I've reproduced the above behavior on R-devel (r74677) and 3.4.3. Below is my session info and locale info for my system with the 3.4.3 version: > sessionInfo() R version 3.4.3 (2017-11-30) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 16.04.3 LTS Matrix products: default BLAS: /usr/lib/libblas/libblas.so.3.6.0 LAPACK: /usr/lib/lapack/liblapack.so.3.6.0 locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base loaded via a namespace (and not attached): [1] compiler_3.4.3 > Sys.getlocale() [1] "LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C" Thanks for your time, Scott P.S. Note that I had posted this question to r-devel, which was the incorrect choice. For archival purposes, I reference the thread here: https://www.mail-archive.com/search?l=mid=20180501185750.445oub53vcdnyyyx%40steph -- Scott Kostyshak Assistant Professor of Economics University of Florida https://people.clas.ufl.edu/skostyshak/ # Ch?vez quantile_type <- 4 __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] Newbie - Scrape Data From PDFs?
Hello, I’m new to R and am using it with RStudio to learn the language. I’m doing so as I have quite a lot of traffic data I would like to explore. My problem is that all the data is located on a number of PDFs. Can someone point me to info on gathering data from other sources? I’ve been to the R FAQ and didn’t see anything and would appreciate your thoughts. I am quite sure now that often, very often, in matters concerning religion and politics a man's reasoning powers are not above the monkey's. -- Mark Twain __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] convertTime package.
To whom it might concern. I am working on a project that needs the convertTime function. I am currently using version 3.4.1 and it says not available for the version. Two questions is there a work around for the function or is there another package that contains that functions. Thanks, Scott Anderwald __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] ggplot's aes_ doesn't work as expected for x=factor()
I'm learning how to use ggplot in a programming approach where I supply variable names on the fly. Using aes_ doesn't work as I expected when x = factor(x). Is this a bug or am I not understanding something? # toy dataset df <- data.frame(LogicalVar = c(FALSE, TRUE, FALSE, TRUE, FALSE, TRUE), Var1 = c(0.01, 0.01, 1, 1, 30, 30), pct1 = c(12, 88, 60, 40, 93, 7), Var2 = c(2, 2, 4, 4, 8, 8), pct2 = c(43, 57, 10, 90, 50, 50) ) varnames <- names(df) # using aes()--this works ggplot(df, aes(x=factor(Var1), y=pct1, fill=LogicalVar)) + geom_bar(stat="identity") # works # using aes_() works in this instance ggplot(df, aes_(x=as.name(varnames[2]), y=as.name(varnames[3]), fill=as.name(varnames[1]))) + geom_bar(stat="identity") # works # it doesn't work here, where only change is using x=factor() ggplot(df, aes_(x=factor(as.name(varnames[2])), y=as.name(varnames[3]), fill=as.name(varnames[1]))) + geom_bar(stat="identity") # doesn't work Error in unique.default(x, nmax = nmax) : unique() applies only to vectors # aes_ does work if I make the x variable a factor ahead of time df2 <- df df2$Var1 <- as.factor(df2$Var1) ggplot(df2, aes_(x=as.name(varnames[2]), y=as.name(varnames[3]), fill=as.name(varnames[1]))) + geom_bar(stat="identity") # works Regards, Scott Waichler Pacific Northwest National Laboratory Richland, WA, USA __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] plot3D color ramp not working as expected
Hi, I want to use a discrete color ramp with plot3D, and show NA values as white (default). I get unexpected results per the following. # as in help(slice3D) example: par(mfrow = c(2,2)) x <- y <- z <- seq(-1, 1, by = 0.1) grid <- mesh(x, y, z) colvar <- with(grid, x*exp(-x^2 - y^2 - z^2)) slice3D (x, y, z, colvar = colvar, theta = 60) # # use three discrete classes and colors instead of a continuous ramp slice3D(x, y, z, colvar = colvar, theta = 60, col = c("blue", "green", "red"), breaks = c(-0.5, -0.1, 0.1, 0.5)) # now set a vertical slice of the cube to NA colvar[10,,] <- NA # displays as expected; default NAcol = "white" slice3D (x, y, z, colvar = colvar, theta = 60) # does not display as expected--notice # the colors shifted down in value, with NA and -0.5 to -0.1 now both white. slice3D(x, y, z, colvar = colvar, theta = 60, col = c("blue", "green", "red"), breaks = c(-0.5, -0.1, 0.1, 0.5)) Please help. Thanks, Scott Scott Waichler, PhD Pacific Northwest National Laboratory scott.waich...@pnnl.gov Richland, Washington, USA __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] customizing color key with plot3D
Karline, Thank you for your help. I discovered that in addition to including clim, I needed to omit breaks. This code uses one of your other examples as a starting point and works as intended: persp3D(z = volcano, zlim = c(-60, 200), phi = 20, colkey = list(length = 0.2, width = 0.4, shift = 0.15, cex.axis = 0.8, cex.clab = 0.85), lighting = TRUE, lphi = 90, clab = c("","height","m"), bty = "f", plot = FALSE) elev.classes <- matrix(findInterval(volcano, vec = seq(50, 200, by=50)), nrow=nrow(volcano), ncol=ncol(volcano)) class.colors <- c("red", "blue", "green") # add as image with own color key, at bottom image3D(z = -60, colvar = elev.classes, add = TRUE, col = class.colors, #breaks = seq(0.5, 3.5, by=1), clim=c(1,3), colkey = list(length = 0.2, width = 0.4, shift = -0.15, cex.axis = 0.8, cex.clab = 0.85, addlines=TRUE, tick=FALSE, at = c(1.33, 2, 2.66), labels=paste("Class", 1:3)), clab = c("","Elev Classes"), plot = TRUE) Your package plot3D is a huge help to me. Previously I had to use a completely different software, VisIt, to do these kinds of composite plots. R is my preferred tool and it is great to finally be able to do these at "home". Thanks, Scott Waichler scott.waich...@pnnl.gov Pacific Northwest National Laboratory Richland, Washington, USA > -Original Message----- > From: Karline Soetaert [mailto:karline.soeta...@nioz.nl] > Sent: Wednesday, June 21, 2017 4:16 AM > To: Waichler, Scott R > Subject: RE: customizing color key with plot3D > > There is an example in the colkey help file: > > example(colkey) > > will show it ( working on the iris dataset) > > I think the basic thing is that you can use "at" to position the labels but > then > you also have to specify clim , i.e. "at = c(1.33, 2, 2.66), clim = > c(0.5,3.5), col = > jetc.col(3)" > > Here is the example: > > with(iris, scatter3D(x = Sepal.Length, y = Sepal.Width, > z = Petal.Length, colvar = as.integer(Species), >col = c("orange", "green", "lightblue"), pch = 16, cex = 2, > clim = c(1, 3), ticktype = "detailed", phi = 20, > xlab = "Sepal Length", ylab = "Sepal Width", > zlab = "Petal Length", main = "iris", > colkey = list(at = c(1.33, 2, 2.66), side = 1, >addlines = TRUE, length = 0.5, width = 0.5, >labels = c("setosa", "versicolor", "virginica") ))) > > > hope it helps, > > > Karline > > -Original Message- > From: Waichler, Scott R [mailto:scott.waich...@pnnl.gov] > Sent: woensdag 21 juni 2017 2:01 > To: R. Help <r-help@r-project.org> > Cc: Karline Soetaert <karline.soeta...@nioz.nl> > Subject: customizing color key with plot3D > > Hi, I am doing composite plots with the package plot3D. One of my > variables is qualitative and indexed to integers, and I would like the legend > for it to have labels located at the integer values (midpoints), and not at > the > breaks between classes. In the example below, the Elev Classes legend has > labels at the breaks and nothing at the midpoints. How can I show the class > labels at 1:3, and not the breaks? > > library(plot3D) > persp3D(z = volcano, zlim = c(-60, 200), phi = 20, > colkey = list(length = 0.2, width = 0.4, shift = 0.15, > cex.axis = 0.8, cex.clab = 0.85), lighting = TRUE, lphi = 90, > clab = c("","height","m"), bty = "f", plot = FALSE) # classify the > volcano > elevations with 3 classes elev.classes <- matrix(findInterval(volcano, vec = > seq(50, 200, by=50)), nrow=nrow(volcano), ncol=ncol(volcano)) class.colors > <- c("red", "blue", "green") # add as image with own color key, at bottom > image3D(z = -60, colvar = elev.classes, add = TRUE, > col = class.colors, breaks = seq(0.5, 3.5, by=1), > colkey = list(length = 0.2, width = 0.4, shift = -0.15, > cex.axis = 0.8, cex.clab = 0.85, addlines=TRUE, > tick=FALSE, > at = 1:3, labels=paste("Class", 1:3)), > clab = c("","Elev Classes"), plot = TRUE) > > Thanks, > Scott > > Scott Waichler > Pacific Northwest National Laboratory > Richland, Washington, USA > scott.waich...@pnnl.gov __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] customizing color key with plot3D
Hi, I am doing composite plots with the package plot3D. One of my variables is qualitative and indexed to integers, and I would like the legend for it to have labels located at the integer values (midpoints), and not at the breaks between classes. In the example below, the Elev Classes legend has labels at the breaks and nothing at the midpoints. How can I show the class labels at 1:3, and not the breaks? library(plot3D) persp3D(z = volcano, zlim = c(-60, 200), phi = 20, colkey = list(length = 0.2, width = 0.4, shift = 0.15, cex.axis = 0.8, cex.clab = 0.85), lighting = TRUE, lphi = 90, clab = c("","height","m"), bty = "f", plot = FALSE) # classify the volcano elevations with 3 classes elev.classes <- matrix(findInterval(volcano, vec = seq(50, 200, by=50)), nrow=nrow(volcano), ncol=ncol(volcano)) class.colors <- c("red", "blue", "green") # add as image with own color key, at bottom image3D(z = -60, colvar = elev.classes, add = TRUE, col = class.colors, breaks = seq(0.5, 3.5, by=1), colkey = list(length = 0.2, width = 0.4, shift = -0.15, cex.axis = 0.8, cex.clab = 0.85, addlines=TRUE, tick=FALSE, at = 1:3, labels=paste("Class", 1:3)), clab = c("","Elev Classes"), plot = TRUE) Thanks, Scott Scott Waichler Pacific Northwest National Laboratory Richland, Washington, USA scott.waich...@pnnl.gov __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] how apply.monthly() in package xts works
Hi, I found that apply.monthly() in xts does not work as I expected in the case of a sparse timeseries: my.dates <- as.Date(c("1992-06-01", "1992-06-24", "1992-06-30", "1993-06-22", "1994-06-07", "1995-06-08")) my.xts <- xts(1:6, my.dates) start(my.xts) # "1992-06-24" end(my.xts) # "1995-06-08" apply.monthly(my.xts, mean) # [,1] # 1995-06-08 3.5 The endpoints it chooses are based on looking at the month (June) alone. I was able to get a value for each (month, year) in the timeseries with the following use of aggregate(): my.months <- months(my.dates) my.years <- years(my.dates) df1 <- data.frame(x = coredata(my.xts), dates = my.dates, months = my.months, years = my.years) df2 <- aggregate(df1[-c(3,4)], df1[c("months", "years")], mean) xts(df2$x, df2$dates) #[,1] # 1992-06-182 # 1993-06-224 # 1994-06-075 # 1995-06-086 Two questions: 1) Is there a more elegant way to do this? 2) Shouldn't the xts documentation discuss the problem of sparse data? Regards, Scott Waichler Pacific Northwest National Laboratory Richland, WA USA __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] odd behavior of numeric()
Why does: > numeric(0.2*25) return [1] 0 0 0 0 0 but > numeric((1-0.8)*25) returns [1] 0 0 0 0 [running version 3.2.0] [Apologies if this has been asked before - it's a hard question to find specific search terms for] Thanks, Scott __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] R packages for Mac Users
Hello, Does anyone know if all the R packages that are available for Windows users are also available for Mac users? Thank you, Scott -- Scott R. Colwell, Ph.D. Associate Professor, Dept. of Mkt/Cons Studies Adjunct Professor, Dept. of Psychology University of Guelph Guelph, Ontario, Canada, N1G 2W1 __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] Output In R
?write.csv and look at with the editor of choice. On Tue, Aug 18, 2015 at 6:41 AM, Shivi82 shivibha...@ymail.com wrote: Hello All, As i am a newbie in R so most of you would have seen this question zillion times. I searched for the answer on this forum as well on other various forums however could not find the answer i am looking for. I am dplyr package and used a very basic code: select(june,city,state,mod) The data sheet i am using has more than 3 million observations but the console does not print all of them and show only few options and give a message: [ reached getOption(max.print) -- omitted 376341 rows ] What is the option that i need to add to see all values in the output. Similarly once i scroll down and then if i scroll up i am not able to see the values starting from row #1. Please suggest -- View this message in context: http://r.789695.n4.nabble.com/Output-In-R-tp4711227.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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 -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] How to 2d cross sections from a 3d finite element mesh
Hi, I have a 3D finite element mesh where each element (cell) is defined by 8 vertices. Each element is a regular polyhedron. The overall domain is a block in shape, but its horizontal principal axes are not coincident with x and y (i.e. the domain is rotated about the z-axis). I want to plot 2D cross sections of discrete and continuous values assigned to the elements. I can think of two ways to go about providing the values to plot: 1) Use the cross section plane intersection with the elements to define a 2D polygon for each intersected element, and plot each as a polygon, with the final product being a mosaic of polygons within the plane of intersection; 2) interpolate to a regular grid and then plot that. Method 1 seems preferable to plotting discrete variables, while the second would be better for contouring a continuous variable. I know how to do a 3D interpolation using Delaunay triangulation, but I wonder if there is a package out there to simplify things. I don't know at all how to go about doing it the first way. Can anyone suggest or point me to existing methods? Thanks, Scott Waichler Pacific Northwest National Laboratory Richland, WA USA __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Can't install rgl: installed package can't be loaded; 'memory not mapped'
/include -g -O2 -fvisibility=hidden -fpic -g -O2 -c select.cpp -o select.o g++ -I/files3/R/R-3.2.1_install/lib64/R/include -DNDEBUG -DHAVE_PNG_H -I/usr/include/libpng12 -DHAVE_FREETYPE -Iext/ftgl -I/usr/include/freetype2 -Iext -I/usr/local/include -g -O2 -fvisibility=hidden -fpic -g -O2 -c subscene.cpp -o subscene.o subscene.cpp: In member function 'void rgl::Subscene::setupViewport(rgl::RenderContext*)': subscene.cpp:736: warning: converting to 'int' from 'double' subscene.cpp:737: warning: converting to 'int' from 'double' subscene.cpp:738: warning: converting to 'int' from 'double' subscene.cpp:739: warning: converting to 'int' from 'double' subscene.cpp:741: warning: converting to 'int' from 'double' subscene.cpp:742: warning: converting to 'int' from 'double' subscene.cpp:743: warning: converting to 'int' from 'double' subscene.cpp:744: warning: converting to 'int' from 'double' g++ -I/files3/R/R-3.2.1_install/lib64/R/include -DNDEBUG -DHAVE_PNG_H -I/usr/include/libpng12 -DHAVE_FREETYPE -Iext/ftgl -I/usr/include/freetype2 -Iext -I/usr/local/include -g -O2 -fvisibility=hidden -fpic -g -O2 -c win32gui.cpp -o win32gui.o g++ -I/files3/R/R-3.2.1_install/lib64/R/include -DNDEBUG -DHAVE_PNG_H -I/usr/include/libpng12 -DHAVE_FREETYPE -Iext/ftgl -I/usr/include/freetype2 -Iext -I/usr/local/include -g -O2 -fvisibility=hidden -fpic -g -O2 -c win32lib.cpp -o win32lib.o g++ -I/files3/R/R-3.2.1_install/lib64/R/include -DNDEBUG -DHAVE_PNG_H -I/usr/include/libpng12 -DHAVE_FREETYPE -Iext/ftgl -I/usr/include/freetype2 -Iext -I/usr/local/include -g -O2 -fvisibility=hidden -fpic -g -O2 -c x11gui.cpp -o x11gui.o g++ -I/files3/R/R-3.2.1_install/lib64/R/include -DNDEBUG -DHAVE_PNG_H -I/usr/include/libpng12 -DHAVE_FREETYPE -Iext/ftgl -I/usr/include/freetype2 -Iext -I/usr/local/include -g -O2 -fvisibility=hidden -fpic -g -O2 -c x11lib.cpp -o x11lib.o g++ -shared -L/files3/R/R-3.2.1_install/lib64/R/lib -L/usr/local/lib64 -o rgl.so ABCLineSet.o BBoxDeco.o Background.o ClipPlane.o Color.o Disposable.o Light.o LineSet.o LineStripSet.o Material.o NULLgui.o PlaneSet.o PointSet.o PrimitiveSet.o RenderContext.o Shape.o SphereMesh.o SphereSet.o SpriteSet.o String.o Surface.o TextSet.o Texture.o Viewpoint.o api.o assert.o callbacks.o device.o devicemanager.o fps.o ftgl.o geom.o gl2ps.o glErrors.o glgui.o gui.o init.o par3d.o pixmap.o platform.o pretty.o render.o rglmath.o rglview.o scene.o select.o subscene.o win32gui.o win32lib.o x11gui.o x11lib.o -lGLU -lGL -L/usr/lib64 -lpng12 -lX11 -lfreetype -L/files3/R/R-3.2.1_install/lib64/R/lib -lR installing to /files3/R/R-3.2.1_install/rgl.Rcheck/rgl/libs ** R ** demo ** inst ** preparing package for lazy loading ** help *** installing help indices ** building package indices ** installing vignettes ** testing if installed package can be loaded sh: line 1: 11949 Segmentation fault '/files3/R/R-3.2.1_install/lib64/R/bin/R' --no-save --slave 21 '/tmp/RtmpQCpp6N/file2b115f4f8e1d' *** caught segfault *** address (nil), cause 'memory not mapped' aborting ... ERROR: loading failed * removing '/files3/R/R-3.2.1_install/rgl.Rcheck/rgl' Scott Waichler Pacific Northwest National Laboratory Richland, WA USA __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Can't install rgl: installed package can't be loaded; 'memory not mapped'
Hi, I can't install package rgl. The last lines from the install process talking about the error are: I'd guess you have an OpenGL problem. Does glxgears run? Yes, it does. I wasn't aware of the program before you mentioned it, but a display opens with 3 gears and here is some output: [root@hokulea R-3.2.1_install]# glxgears -info GL_RENDERER = Mesa GLX Indirect GL_VERSION= 1.2 (1.5 Mesa 6.5.1) GL_VENDOR = Mesa project: www.mesa3d.org GL_EXTENSIONS = GL_ARB_depth_texture GL_ARB_imaging GL_ARB_multitexture GL_ARB_point_parameters GL_ARB_point_sprite GL_ARB_shadow GL_ARB_texture_border_clamp GL_ARB_texture_cube_map GL_ARB_texture_env_add GL_ARB_texture_env_combine GL_ARB_texture_env_crossbar GL_ARB_texture_env_dot3 GL_ARB_texture_mirrored_repeat GL_ARB_texture_non_power_of_two GL_ARB_window_pos GL_EXT_abgr GL_EXT_bgra GL_EXT_blend_color GL_EXT_blend_func_separate GL_EXT_blend_minmax GL_EXT_blend_subtract GL_EXT_draw_range_elements GL_EXT_framebuffer_object GL_EXT_fog_coord GL_EXT_multi_draw_arrays GL_EXT_packed_pixels GL_EXT_rescale_normal GL_EXT_secondary_color GL_EXT_separate_specular_color GL_EXT_shadow_funcs GL_EXT_stencil_wrap GL_EXT_texture3D GL_EXT_texture_edge_clamp GL_EXT_texture_env_add GL_EXT_texture_env_combine GL_EXT_texture_env_dot3 GL_EXT_texture_lod_bias GL_EXT_texture_object GL_EXT_vertex_array GL_ATI_texture_mirror_once GL_IBM_texture_mirrored_repeat GL_NV_blend_square GL_NV_texture_rectan! gle GL_NV_texgen_reflection GL_SGIS_generate_mipmap GL_SGIS_texture_lod GL_SGIX_depth_texture GL_SGIX_shadow 24839 frames in 6.0 seconds = 4153.389 FPS 7152 frames in 6.0 seconds = 1192.594 FPS . . . Scott ** testing if installed package can be loaded sh: line 1: 11949 Segmentation fault '/files3/R/R- 3.2.1_install/lib64/R/bin/R' --no-save --slave 21 '/tmp/RtmpQCpp6N/file2b115f4f8e1d' *** caught segfault *** address (nil), cause 'memory not mapped' aborting ... ERROR: loading failed I realize this is probably not an R problem, but a Google search turns up nothing that helps, and I'm hoping someone here can help anyway. Below are my sessionInfo() output and the contents of the first file generated with R CMD check rgl_0.95.1247.tar.gz. sessionInfo() R version 3.2.1 (2015-06-18) Platform: x86_64-unknown-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base loaded via a namespace (and not attached): [1] tools_3.2.1 # R CMD check rgl_0.95.1247.tar.gz * installing *source* package 'rgl' ... ** package 'rgl' successfully unpacked and MD5 sums checked checking for gcc... gcc -std=gnu99 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 -std=gnu99 accepts -g... yes checking for gcc -std=gnu99 option to accept ISO C89... none needed checking how to run the C preprocessor... gcc -std=gnu99 -E checking for gcc... (cached) gcc -std=gnu99 checking whether we are using the GNU C compiler... (cached) yes checking whether gcc -std=gnu99 accepts -g... (cached) yes checking for gcc -std=gnu99 option to accept ISO C89... (cached) none needed checking whether __attribute__((visibility())) is supported... yes checking whether gcc -std=gnu99 accepts -fvisibility... yes checking whether accepts -fvisibility... no checking for libpng-config... yes configure: using libpng-config configure: using libpng dynamic linkage checking for X... libraries , headers checking GL/gl.h usability... yes checking GL/gl.h presence... yes checking for GL/gl.h... yes checking GL/glu.h usability... yes checking GL/glu.h presence... yes checking for GL/glu.h... yes checking for glEnd in -lGL... yes checking for gluProject in -lGLU... yes checking for freetype-config... yes configure: using Freetype and FTGL configure: creating ./config.status config.status: creating src/Makevars ** libs g++ -I/files3/R/R-3.2.1_install/lib64/R/include -DNDEBUG -DHAVE_PNG_H - I/usr/include/libpng12 -DHAVE_FREETYPE -Iext/ftgl - I/usr/include/freetype2 -Iext -I/usr/local/include -g -O2 -fvisibility=hidden -fpic -g -O2 -c ABCLineSet.cpp -o ABCLineSet.o g++ -I/files3/R/R-3.2.1_install/lib64/R/include -DNDEBUG -DHAVE_PNG_H - I/usr/include/libpng12 -DHAVE_FREETYPE -Iext/ftgl - I/usr/include/freetype2 -Iext -I/usr/local/include -g -O2 -fvisibility=hidden -fpic -g -O2 -c
[R] write.table with append=T after using cat on same file
Hi, For years I've been writing text to the beginning of files with cat(append=F) , then following that text with data written by write.table(append=T). It is now giving me an error message. I'm using R-3.1.2. What gives? df - data.frame(x = 1, y = 1:10, z = 10:1) cat(file=junk.txt, sep=, # An introductory note.\n) write.table(df, file=junk.txt, sep=,, append=T, quote=F, row.names=F, col.names=F) Error in file(file, ifelse(append, a, w)) : invalid 'open' argument Thanks, Scott Waichler Pacific Northwest National Laboratory Richland, WA USA __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] Looping and break
Hello, I apologies for bringing up next and break in loops given that there is so much on the net about it, but I've tried numerous examples found using Google and just can't seem to get this to work. This is a simple version of what I am doing with matrices but it shows the issue. I need to have the loop indexed as n to perform a calculation on the variable total. But if total is greater than 8, it goes to the next loop indexed a. For example, it does condition a = 1 for n = 1 to 50 but within n if total is greater than 8 it goes to the next condition of a which would be a = 2, and so on. for (a in 1:3){ if (a == 1) { b - c(1:5) } if (a == 2) { b - c(1:5) } if (a == 3) { b - c(1:5) } for (n in 1:50){ if (n 15) next total - 2*b if (total 8) next } } Any help would be greatly appreciated. Thanks, Scott -- View this message in context: http://r.789695.n4.nabble.com/Looping-and-break-tp4704093.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] Saving Mean Relative Difference from all.equal()
I think I have one solution. Not very pretty though. Relies on the text not changing at all. as.numeric(gsub(Mean relative difference: , , all.equal(cov2cor(ITEMCOV),cor(item.data))[2])) Is there a better way? -- View this message in context: http://r.789695.n4.nabble.com/Saving-Mean-Relative-Difference-from-all-equal-tp4703905p4703908.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] Saving Mean Relative Difference from all.equal()
Hello, Does anyone know how to save the numeric value of the mean relative difference when using the all.equal() command? For example this: all.equal(cov2cor(ITEMCOV),cor(item.data)) Gives: [1] Attributes: Length mismatch: comparison on first 1 components [2] Mean relative difference: 0.01523708 I'd like to save the value 0.01523708 in a numeric format. Thanks, -- View this message in context: http://r.789695.n4.nabble.com/Saving-Mean-Relative-Difference-from-all-equal-tp4703905.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Extracting Factor Pattern Matrix Similar to Proc Factor
Thanks everyone -- View this message in context: http://r.789695.n4.nabble.com/Extracting-Factor-Pattern-Matrix-Similar-to-Proc-Factor-tp4703704p4703904.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Extracting Factor Pattern Matrix Similar to Proc Factor
Thanks David. What do you do when the input is a covariance matrix rather than a dataset? -- View this message in context: http://r.789695.n4.nabble.com/Extracting-Factor-Pattern-Matrix-Similar-to-Proc-Factor-tp4703704p4703719.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] Extracting Factor Pattern Matrix Similar to Proc Factor
Hello, I am fairly new to R and coming from SAS IML. I am rewriting one of my MC simulations in R and am stuck on extracting a factor pattern matrix as would be done in IML using Proc Factor. I have found the princomp() command and read through the manual but can't seem to figure out how to save the factor pattern matrix. I am waiting for the R for SAS Users book to arrive. What I would use in SAS IML to get at what I am looking for is: PROC FACTOR Data=MODELCOV15(TYPE=COV) NOBS=1 N=16 CORR OUTSTAT=FAC.FACOUT15; RUN; DATA FAC.PATTERN15; SET FAC.FACOUT15; IF _TYPE_='PATTERN'; DROP _TYPE_ _NAME_; RUN; Would any SAS IML to R converts be able to help me with this? Thanks, Scott Colwell, PhD -- View this message in context: http://r.789695.n4.nabble.com/Extracting-Factor-Pattern-Matrix-Similar-to-Proc-Factor-tp4703704.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] Rewriting the Biplot Function
, pc.biplot = FALSE, yCol=NULL, ...) { if (length(choices) != 2L) stop(length of choices must be 2) if (!length(scores - x$x)) stop(gettextf(object '%s' has no scores, deparse(substitute(x))), domain = NA) if (is.complex(scores)) stop(biplots are not defined for complex PCA) lam - x$sdev[choices] n - NROW(scores) lam - lam * sqrt(n) if (scale 0 || scale 1) warning('scale' is outside [0, 1]) if (scale != 0) lam - lam^scale else lam - 1 if (pc.biplot) lam - lam/sqrt(n) colouredBiplot.internal(t(t(scores[, choices])/lam), t(t(x$rotation[, choices]) * lam), yCol, ...) invisible() } I have looked into a few alternatives but most either don't allow that type of different colouring, or else aren't compatable with my various versions of R. Help me R-Help, you're my only hope, Scott [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] Rewriting the Biplot Function
Boris, Thanks very much, this looks just like what I need! All the best, Scott From: Boris Steipe [boris.ste...@utoronto.ca] Sent: 20 January 2015 20:21 To: Scott Robinson Cc: r-help@r-project.org Subject: Re: [R] Rewriting the Biplot Function Scott - A few months ago I posted on this list a modified version of biplot that takes a colour argument (and preserves the axes information, so you can use points() ). I don't have time right now to experiment and look at your code, but perhaps this does out of the box what you need. Cheers, Boris Previous post == Since I've wanted this capability for some time, I modified the original biplot() to accept a type parameter type={t (Default) | p | n}. For t, the function behaves almost exactly as before. For p it plots points, and should accept all the usual arguments for that. For n it plots nothing except the axes. You can then add the points as desired. I also added two parameters col.arrows = red, and col.text = black to have extra control. Here is an example. (Note, you have to load the function, below, first. library(MASS) data(crabs) PRC - prcomp(crabs[, 4:8]) myBiplot(PRC) myBiplot(PRC, choices=2:3, cex = 0.7, col.text=#445599) # much as before # use filled points, color by the value found in column 4 of the data r - range(crabs[,4]) n - 15 cv - cm.colors(n) c - cv[cut(crabs[,4],n)] myBiplot(PRC, choices=2:3, type = p, pch=20, col=c, col.arrows = #FF6600) # finally: plot nothing then use points() for detailed control myBiplot(PRC, choices=2:3, type = n) # no points # blue/orange crab males/females - as given by rows in the data points(PRC$x[ 1: 50, 2:3], pch=21, bg=#0066FF) points(PRC$x[ 51:100, 2:3], pch=24, bg=#0066FF) points(PRC$x[101:150, 2:3], pch=21, bg=#FF6600) points(PRC$x[151:200, 2:3], pch=24, bg=#FF6600) == myBiplot - function (x, choices = 1L:2L, scale = 1, pc.biplot = FALSE, var.axes = TRUE, type = t, col, col.arrows = #FF, col.text = #00, cex = rep(par(cex), 2), expand = 1, xlabs = NULL, ylabs = NULL, xlim = NULL, ylim = NULL, main = NULL, sub = NULL, xlab = NULL, ylab = NULL, arrow.len = 0.1, ... ) { if (length(choices) != 2L) stop(length of choices must be 2) if (!length(scores - x$x)) stop(gettextf(object '%s' has no scores, deparse(substitute(x))), domain = NA) if (is.complex(scores)) stop(biplots are not defined for complex PCA) lam - x$sdev[choices] n - NROW(scores) lam - lam * sqrt(n) if (scale 0 || scale 1) warning('scale' is outside [0, 1]) if (scale != 0) lam - lam^scale else lam - 1 if (pc.biplot) lam - lam/sqrt(n) y - t(t(x$rotation[, choices]) * lam) x - t(t(scores[, choices])/lam) # note that from here on # x is no longer the PC object # originally pased into the function n - nrow(x) p - nrow(y) if (missing(xlabs)) { xlabs - dimnames(x)[[1L]] if (is.null(xlabs)) xlabs - 1L:n } xlabs - as.character(xlabs) dimnames(x) - list(xlabs, dimnames(x)[[2L]]) if (missing(ylabs)) { ylabs - dimnames(y)[[1L]] if (is.null(ylabs)) ylabs - paste(Var, 1L:p) } ylabs - as.character(ylabs) dimnames(y) - list(ylabs, dimnames(y)[[2L]]) if (length(cex) == 1L) cex - c(cex, cex) unsigned.range - function(x) c(-abs(min(x, na.rm = TRUE)), abs(max(x, na.rm = TRUE))) rangx1 - unsigned.range(x[, 1L]) rangx2 - unsigned.range(x[, 2L]) rangy1 - unsigned.range(y[, 1L]) rangy2 - unsigned.range(y[, 2L]) if (missing(xlim) missing(ylim)) xlim - ylim - rangx1 - rangx2 - range(rangx1, rangx2) else if (missing(xlim)) xlim - rangx1 else if (missing(ylim)) ylim - rangx2 ratio - max(rangy1/rangx1, rangy2/rangx2)/expand on.exit(par(op)) op - par(pty = s) if (!is.null(main)) op - c(op, par(mar = par(mar) + c(0, 0, 1, 0))) # first, plot scores - either normally, or as row labels if (type == p) { plot(x, type = type, xlim = xlim, ylim = ylim, col = col, xlab = xlab, ylab = ylab, sub = sub, main = main, ...) } else if (type == t) { plot(x, type = n, xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, sub = sub, main = main, ...) text(x, xlabs, cex = cex[1L], col = col.text, ...) } else if (type == n) { # plot an empty frame plot(x, type = type, xlim = xlim, ylim = ylim
Re: [R] Gender balance in R
On Mon, Nov 24, 2014 at 12:34 PM, Sarah Goslee sarah.gos...@gmail.com wrote: I took a look at apparent gender among list participants a few years ago: https://stat.ethz.ch/pipermail/r-help/2011-June/280272.html Same general thing: very few regular participants on the list were women. I don't see any sign that that has changed in the last three years. The bar to participation in the R-help list is much, much lower than that to become a developer. I plotted the gender of posters on r-help over time. The plot is here: https://twitter.com/scottkosty/status/449933971644633088 The code to reproduce that plot is here: https://github.com/scottkosty/genderAnalysis The R file there will call devtools::install_github to install a package from Github used for guessing the gender based on the first name (https://github.com/scottkosty/gender). Note also on that tweet that Gabriela de Queiroz posted it, who is the founder of R-ladies; and that David Smith showed interest in discussing the topic. So there is definitely demand for some data analysis and discussion on the topic. It would be interesting to look at the stats for CRAN packages as well. The very low percentage of regular female participants is one of the things that keeps me active on this list: to demonstrate that it's not only men who use R and participate in the community. Thank you for that! Scott -- Scott Kostyshak Economics PhD Candidate Princeton University (If you decide to do the stats for 2014, be aware that I've been out on medical leave for the past two months, so the numbers are even lower than usual.) Sarah On Mon, Nov 24, 2014 at 10:10 AM, Maarten Blaauw maarten.bla...@qub.ac.uk wrote: Hi there, I can't help to notice that the gender balance among R developers and ordinary members is extremely skewed (as it is with open source software in general). Have a look at http://www.r-project.org/foundation/memberlist.html - at most a handful of women are listed among the 'supporting members', and none at all among the 29 'ordinary members'. On the other hand I personally know many happy R users of both genders. My questions are thus: Should R developers (and users) be worried that the 'other half' is excluded? If so, how could female R users/developers be persuaded to become more visible (e.g. added as supporting or ordinary members)? Thanks, Maarten -- Sarah Goslee http://www.functionaldiversity.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. __ 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] Gender balance in R
On Tue, Nov 25, 2014 at 8:24 AM, Maarten Blaauw maarten.bla...@qub.ac.uk wrote: Nice graph, Scott, thanks! Based on your code I plotted not the absolute numbers but the ratios, which show slowly increasing relative participation of female Rhelpers over time (red = women, blue=men, black=unknown). After a c. 5% female contribution in 1998, this has grown to about 15% now. At this rate we'll reach parity around AD 2080. Interesting forecasts Maarten! Let's hope for a trend break to make them wrong. Scott -- Scott Kostyshak Economics PhD Candidate Princeton University My code: if (!require(gender)) { library(devtools) install_github(scottkosty/gender) library(gender) } rHelp - rHelpNames rHelp[is.na(rHelp$gender), gender] - unknown yr - unique(rHelp$year) helpers - list(dates, M=rep(0, length(yr)), F=rep(0, length(yr)), unkn=rep(0, length(yr))) for(i in 1:nrow(rHelp)) { j - which(yr == rHelp$year[i]) gender - rHelp$gender[i] if(gender == M) helpers$M[[j]] - helpers$M[[j]]+1 else if(gender == F) helpers$F[[j]] - helpers$F[[j]]+1 else if(gender == unknown) helpers$unkn[[j]] - helpers$unkn[[j]]+1 } plot(yr, helpers$M / (helpers$M+helpers$F+helpers$unkn), type=l, col=4, ylim=c(0,1), ylab=proportions, yaxs=i) lines(yr, helpers$F / (helpers$M+helpers$F+helpers$unkn), col=2) lines(yr, helpers$unkn / (helpers$M+helpers$F+helpers$unkn)) Cheers, Maarten On 25/11/14 12:11, Scott Kostyshak wrote: On Mon, Nov 24, 2014 at 12:34 PM, Sarah Goslee sarah.gos...@gmail.com wrote: I took a look at apparent gender among list participants a few years ago: https://stat.ethz.ch/pipermail/r-help/2011-June/280272.html Same general thing: very few regular participants on the list were women. I don't see any sign that that has changed in the last three years. The bar to participation in the R-help list is much, much lower than that to become a developer. I plotted the gender of posters on r-help over time. The plot is here: https://twitter.com/scottkosty/status/449933971644633088 The code to reproduce that plot is here: https://github.com/scottkosty/genderAnalysis The R file there will call devtools::install_github to install a package from Github used for guessing the gender based on the first name (https://github.com/scottkosty/gender). Note also on that tweet that Gabriela de Queiroz posted it, who is the founder of R-ladies; and that David Smith showed interest in discussing the topic. So there is definitely demand for some data analysis and discussion on the topic. It would be interesting to look at the stats for CRAN packages as well. The very low percentage of regular female participants is one of the things that keeps me active on this list: to demonstrate that it's not only men who use R and participate in the community. Thank you for that! Scott -- Scott Kostyshak Economics PhD Candidate Princeton University (If you decide to do the stats for 2014, be aware that I've been out on medical leave for the past two months, so the numbers are even lower than usual.) Sarah On Mon, Nov 24, 2014 at 10:10 AM, Maarten Blaauw maarten.bla...@qub.ac.uk wrote: Hi there, I can't help to notice that the gender balance among R developers and ordinary members is extremely skewed (as it is with open source software in general). Have a look at http://www.r-project.org/foundation/memberlist.html - at most a handful of women are listed among the 'supporting members', and none at all among the 29 'ordinary members'. On the other hand I personally know many happy R users of both genders. My questions are thus: Should R developers (and users) be worried that the 'other half' is excluded? If so, how could female R users/developers be persuaded to become more visible (e.g. added as supporting or ordinary members)? Thanks, Maarten -- Sarah Goslee http://www.functionaldiversity.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. -- | Dr. Maarten Blaauw | Lecturer in Chronology | | School of Geography, Archaeology Palaeoecology | Queen's University Belfast, UK | | www http://www.chrono.qub.ac.uk/blaauw | tel +44 (0)28 9097 3895 __ 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] Gender balance in R
On Tue, Nov 25, 2014 at 1:15 PM, Martin Morgan mtmor...@fredhutch.org wrote: On 11/25/2014 04:11 AM, Scott Kostyshak wrote: On Mon, Nov 24, 2014 at 12:34 PM, Sarah Goslee sarah.gos...@gmail.com wrote: I took a look at apparent gender among list participants a few years ago: https://stat.ethz.ch/pipermail/r-help/2011-June/280272.html Same general thing: very few regular participants on the list were women. I don't see any sign that that has changed in the last three years. The bar to participation in the R-help list is much, much lower than that to become a developer. I plotted the gender of posters on r-help over time. The plot is here: https://twitter.com/scottkosty/status/449933971644633088 The code to reproduce that plot is here: https://github.com/scottkosty/genderAnalysis The R file there will call devtools::install_github to install a package from Github used for guessing the gender based on the first name (https://github.com/scottkosty/gender). It would be great to include in your package the script that scraped author names from R-help archives (I guess that's what you did?). Presumably it easily applies to other mailing lists hosted at the same location (R-devel, further along the ladder from user to developer, and Bioconductor / Bioc-devel, in a different domain and perhaps confounded with a different 'feel' to the list). Also the R community is definitely international, so finding more versatile gender-assignment approaches seems important. I just put the script up on https://github.com/scottkosty/genderAnalysis I don't have much time at the moment to generalize it, but a pull request is always welcome. Alternatively, anyone is welcome (at least as far as I'm concerned) to take the script and modify it for any purpose. it might be interesting to ask about participation in mailing list forums versus other, and in particular the recent Bioconductor transition from mailing list to 'StackOverflow' style support forum (https://support.bioconductor.org) -- on the one hand the 'gamification' elements might seem to only entrench male participation, while on the other we have already seen increased (quantifiable) and broader (subjective) participation from the Bioconductor community. I'd be happy to make support site usage data available, and am interested in collaborating in an academically well-founded analysis of this data; any interested parties please feel free to contact me off-list. I would be interested in collaborating on such a project in the future also. Scott -- Scott Kostyshak Economics PhD Candidate Princeton University Martin Morgan Bioconductor Note also on that tweet that Gabriela de Queiroz posted it, who is the founder of R-ladies; and that David Smith showed interest in discussing the topic. So there is definitely demand for some data analysis and discussion on the topic. It would be interesting to look at the stats for CRAN packages as well. The very low percentage of regular female participants is one of the things that keeps me active on this list: to demonstrate that it's not only men who use R and participate in the community. Thank you for that! Scott -- Scott Kostyshak Economics PhD Candidate Princeton University (If you decide to do the stats for 2014, be aware that I've been out on medical leave for the past two months, so the numbers are even lower than usual.) Sarah On Mon, Nov 24, 2014 at 10:10 AM, Maarten Blaauw maarten.bla...@qub.ac.uk wrote: Hi there, I can't help to notice that the gender balance among R developers and ordinary members is extremely skewed (as it is with open source software in general). Have a look at http://www.r-project.org/foundation/memberlist.html - at most a handful of women are listed among the 'supporting members', and none at all among the 29 'ordinary members'. On the other hand I personally know many happy R users of both genders. My questions are thus: Should R developers (and users) be worried that the 'other half' is excluded? If so, how could female R users/developers be persuaded to become more visible (e.g. added as supporting or ordinary members)? Thanks, Maarten -- Sarah Goslee http://www.functionaldiversity.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. __ 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. -- Computational Biology / Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 Location: Arnold Building M1 B861 Phone
Re: [R] Symbolic equations to R code?
Ista, On the one hand I'd like it to be as flexible as possible so the students could really come up with whatever they like. On the other hand, restricting their choices probably would make it easier to do the backend. The goal would be to get them to realize that the apparatus of hypothesis testing (when done via simulation/randomization techniques) doesn't depend on what the statistic is. The flow of steps is the same whether the statistic is a mean, variance, or their own kooky thing. Obviously this isn't the end of the story - what the statistic is actually describing is also a crucial component to interpreting the results of a hypothesis test, but I think it teaches an important pedagogical point about where statistics come from and that if they find themselves in a situation in the future where they need to make up their own, then that is perfectly okay. So they could make up ones like: arctan( (max({x})^2)/ (min({x})^2) )-3), max({x})-min({x}), sum from i to n of (x_i - 25th%ile({x}) )^3) [that might be tricky to write in a standard equation editor] none of these is hard to write in R, but translating from an equation editor might be. Perhaps the best solution is (as Alan suggests below) to write a Shiny function builder myself so that I can control the whole process and make sure that they can't enter anything that would break the backend. Or to have them learn the rudiments of writing equations in R so that it bypasses the whole process. Thanks, Scott On 11/19/14 1:52 PM, Ista Zahn wrote: Hi Scott, Can you give a couple of examples of the equations you have in mind along with how those should be translated to R? Thanks, Ista __ 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] Symbolic equations to R code?
I'm looking for a package that would take a mathematical function written in symbolic notation and convert it into R code. What I have in mind would be something like the following: 1) Have a GUI (e.g. something like Microsoft Word equation editor, but computer readable) to produce an equation so that it would look just like it would if you were writing it on paper 2) Have some way to translate the GUI into some markup language (perhaps this would be in a non-R application) 3) Have code that translates the markup into an R function. My motivation is that I teach an intro stats class that uses randomization techniques. I would like for the students to be able to make up their own statistics and then simulate null distributions or estimate bootstrap confidence intervals using them. I'd prefer to make it as easy as possible for them to enter the formula for their statistics into the computer, and I think that will be easier if it looks from their end just like it would if they were writing it on paper. After they submit the equation, I would have an RStudio shiny applet take care of the rest. Any suggestions for R packages that could do this or pieces of it would be very much appreciated. Thanks, Scott __ 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] Symbolic equations to R code?
David, Thanks for your reply and suggestions of packages - let me clarify what I am looking for: The web editor you pointed out is the sort of equation editor I'm looking for and I had seen others like this.Let's say I use that one (step 1) and then have Latex code (step 2). Step (3) would be converting the Latex into an R function automatically. That's the key package. In searching for Latex and R there seem to be lots of ways to embed R code into Latex but not to translate Latex into R. As far as R commander, I do indeed know about that and am not looking for an R GUI for my students. They are not going be learning R in the class. As far as I know, R commander can't do steps 1-3, but please let me know if I am wrong. After students use the interface in step (1), the software in steps (2) and (3) would be for my benefit, to programmatically translate the equations of 150+ students into R code. Thanks, Scott __ 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] Testing general hypotheses on regression coefficients
Hi Chris, On Fri, Sep 5, 2014 at 7:17 PM, Chris bonsxa...@yahoo.com wrote: Hi. Say I have a model like y = a + B1*x1 + B2*x2 + B3*x3 + B4*x4 + e and I want to test H0: B2/B1 = 0 As noted by Bert, think about this. or H0: B2/B1=B4/B3 (whatever H1). How can I proceed? I now about car::linearHypothesis, but I can't figure out a way to do the tests above. Any hint? Take a look at car::deltaMethod. I suggest you study the theory of the delta method. If you happen to have taken a graduate statistics/econometrics class it should not be difficult and can provide some insights. If not, at least consider that the delta method can lead to misleading estimates (biased standard errors) in many cases for finite samples. You might want to run some simulations to get a feel for it. Best, Scott -- Scott Kostyshak Economics PhD Candidate Princeton University __ 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-es] NA no es reconocido como NA
El 13 de agosto de 2014, 17:06, neo ericconchamu...@gmail.com escribió: Estimados, cuál es la diferencia para R entre : NA NA NA NA Para saberlo, usa class, e.g. class(NA) [1] character class(NA) [1] logical class(NA_integer_) [1] integer class(NA) [1] character como demuestra el segundo y tercer ejemplos, hay varios tipos de NA. Saludos, Scott -- Scott Kostyshak Economics PhD Candidate Princeton University ___ R-help-es mailing list R-help-es@r-project.org https://stat.ethz.ch/mailman/listinfo/r-help-es
[R] Reading chunks of data from a file more efficiently
Hi, I have some very large (~1.1 GB) output files from a groundwater model called STOMP that I want to read as efficiently as possible. For each variable there are over 1 million values to read. Variables are not organized in columns; instead they are written out in sections in the file, like this: X-Direction Node Positions, m 5.93145E+05 5.93155E+05 5.93165E+05 5.93175E+05 5.93245E+05 5.93255E+05 5.93265E+05 5.93275E+05 . . . 5.94695E+05 5.94705E+05 5.94715E+05 5.94725E+05 5.94795E+05 5.94805E+05 5.94815E+05 5.94825E+05 Y-Direction Node Positions, m 1.14805E+05 1.14805E+05 1.14805E+05 1.14805E+05 1.14805E+05 1.14805E+05 1.14805E+05 1.14805E+05 . . . 1.17195E+05 1.17195E+05 1.17195E+05 1.17195E+05 1.17195E+05 1.17195E+05 1.17195E+05 1.17195E+05 Z-Direction Node Positions, m 9.55000E+01 9.55000E+01 9.55000E+01 9.55000E+01 9.55000E+01 9.55000E+01 9.55000E+01 9.55000E+01 . . . I want to read and use only a subset of the variables. I wrote the function below to find the line where each target variable begins and then scan the values, but it still seems rather slow, perhaps because I am opening and closing the file for each variable. Can anyone suggest a faster way? # Reads original STOMP plot file (plot.*) directly. Should be useful when the plot files are # very large with lots of variables, and you just want to retrieve a few of them. # Arguments: 1) plot filename, 2) number of nodes, # 3) character vector of names of target variables you want to return. # Returns a list with the selected plot output. READ.PLOT.OUTPUT6 - function(plt.file, num.nodes, var.names) { lines - readLines(plt.file) num.vars - length(var.names) tmp - list() for(i in 1:num.vars) { ind - grep(var.names[i], lines, fixed=T, useBytes=T) if(length(ind) != 1) stop(Not one line in the plot file with matching variable name.\n) tmp[[i]] - scan(plt.file, skip=ind, nmax=num.nodes, quiet=T) } return(tmp) } # end READ.PLOT.OUTPUT6() Regards, Scott Waichler Pacific Northwest National Laboratory Richland, WA, USA scott.waich...@pnnl.gov __ 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] Using axis limits with plot3D
I would like to use the functions in the plot3D package but I am having trouble getting the axis limits to work correctly. The slices plotted by the code below go beyond the bounds of the persp box and obscure the axis information. How can I show just the part of the domain within x.limits and y.limits? library(plot3D) x - z - seq(-4, 4, by=0.2) y - seq(-6, 6, by=0.2) M - mesh(x,y,z) R - with(M, sqrt(x^2 + y^2 +z^2)) p - sin(2*R)/(R+1e-3) x.limits - c(-2, 2) y.limits - c(-2, 2) slice3D(x,y,z, colvar=p, xs=0, ys=c(0, 4), zs=NULL, xlim=x.limits, ylim=y.limits, scale=F, ticktype=detailed) Thanks, Scott Waichler Pacific Northwest National Laboratory Richland, WA, USA __ 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-es] Descargar lista de paquetes zipeados
Hola, Si el objetivo final es replicar algo en el futuro con exactamente las mismas versiones de los paquetes, tal vez le sirve packrat: http://rstudio.github.io/packrat/ Saludos, Scott 2014-07-22 14:03 GMT+10:00 Julio Alejandro Di Rienzo dirienzo.ju...@gmail.com: Hola Alguien sabe como descargar una lista de librerías de R en formato zipeado. Por ejemplo quiero descargar las librerías (lme4, latticeExtras, Biobase,., etc,etc) en formato zipeado. Se que puedo hacerlo una por una desde el cran pero quisiera tener un procedimiento para hacerlo automáticamente. Prof. Julio Di Rienzo Estadística y Biometría FCA- U.N. Córdoba http://sites.google.com/site/juliodirienzo ___ R-help-es mailing list R-help-es@r-project.org https://stat.ethz.ch/mailman/listinfo/r-help-es
[R] how to subset based on other row values and multiplicity
Hi R experts, I have a dataset as sampled below. Values are only regarded as Œconfirmed¹ in an individual (Œid¹) if they occur more than once at least 30 days apart. id date value a2000-01-01 x a2000-03-01 x b2000-11-11 w c2000-11-11 y c2000-10-01 y c2000-09-10 y c2000-12-12 z c2000-10-11 z d2000-11-11 w d2000-11-10 w I wish to subset the data to retain rows where the value for the individual is confirmed more than 30 days apart. So, after deleting all rows with just one occurrence of id and value, the rest would be the earliest occurrence of each value in each case id, provided 31 or more days exist between the dates. If 1 value is present per id, each value level needs to be assessed independently. This example would then reduce to: id date value a2000-01-01 x c2000-09-10 y c2000-10-11 z I can do this via some crude loops and subsetting, but I am looking for as much efficiency as possible as the dataset has around 50 million rows to assess. Any suggestions welcomed. Thanks in advance Scott Williams MD Melbourne, Australia This email (including any attachments or links) may contain confidential and/or legally privileged information and is intended only to be read or used by the addressee. If you are not the intended addressee, any use, distribution, disclosure or copying of this email is strictly prohibited. Confidentiality and legal privilege attached to this email (including any attachments) are not waived or lost by reason of its mistaken delivery to you. If you have received this email in error, please delete it and notify us immediately by telephone or email. Peter MacCallum Cancer Centre provides no guarantee that this transmission is free of virus or that it has not been intercepted or altered and will not be liable for any delay in its receipt. __ 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 subset based on other row values and multiplicity
Thanks guys - amazingly prompt solutions from the R community as always. Yes, the c-y value reverts to just the first date event - the spirit of this is that I am trying to identify and confirm a list of diagnoses that a patient has coded in government administrative data. Once a diagnosis is made and confirmed, I am not interested in whether it is listed again and again later on. I just need that date at which it first became apparent. So in the multiple c-y case, the min date is the correct one. Some cases will have the same diagnosis listed dozens of times, hence the very bloated dataset. Time to churn through the data is not a big issue, so I will have a go with Jim¹s neat code he just sent on perhaps a few thousand rows and see how I get on. S On 17/07/2014 12:09 am, John McKown john.archie.mck...@gmail.com wrote: On Wed, Jul 16, 2014 at 8:51 AM, jim holtman jholt...@gmail.com wrote: I can reproduce what you requested, but there was the question about what happens with the multiple 'c-y' values. require(data.table) x - read.table(text = 'id date value + a2000-01-01 x + a2000-03-01 x + b2000-11-11 w + c2000-11-11 y + c2000-10-01 y + c2000-09-10 y + c2000-12-12 z + c2000-10-11 z + d2000-11-11 w + d2000-11-10 w', as.is = TRUE, header = TRUE) setDT(x) x[, date := as.Date(date)] setkey(x, id, value, date) y - x[ + , { + if (.N == 1) val - NULL # only one -- delete + else { + dif - difftime(tail(date, -1), head(date, -1), units = 'days') + # return first value if any 31 + if (any(dif = 31)) val - list(date = date[1L]) + else val - NULL + } + val + } + , keyby = 'id,value' + ] y id value date 1: a x 2000-01-01 2: c y 2000-09-10 3: c z 2000-10-11 Jim Holtman Data Munger Guru What is the problem that you are trying to solve? Tell me what you want to do, not how you want to do it. Wow, I picked up a couple of _nice_ techniques from that one post! Looks like data.table will let me do SQL like things in R. I have a warped brain. I think in result sets and matrix operations Many thanks. -- There is nothing more pleasant than traveling and meeting new people! Genghis Khan Maranatha! John McKown This email (including any attachments or links) may contain confidential and/or legally privileged information and is intended only to be read or used by the addressee. If you are not the intended addressee, any use, distribution, disclosure or copying of this email is strictly prohibited. Confidentiality and legal privilege attached to this email (including any attachments) are not waived or lost by reason of its mistaken delivery to you. If you have received this email in error, please delete it and notify us immediately by telephone or email. Peter MacCallum Cancer Centre provides no guarantee that this transmission is free of virus or that it has not been intercepted or altered and will not be liable for any delay in its receipt. __ 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 evaluating partitioning around medoids clustering method R clValid package
I have a data.frame with 300 observations of 36 numerical, categorical, and NA variables. I am trying to evaluate the partitioning around medoids clustering algorithm for a marketing segmentation study. My original dataset has over 130,000 observations, but I took a sample for easy reproducibility reasons. My machine Mac OSX 10.9.3: sessionInfo() R version 3.1.0 (2014-04-10) Platform: x86_64-apple-darwin13.1.0 (64-bit) Problem: Getting an error when doing internal and stability evaluation with the clValid CRAN package in R. Code: #Convert csv to data.frame frame -as.data.frame(Smallstore1) library(cluster) #Create dissimilarity matrix #Gower coefficient for finding distance between mixed variables daisy1 - daisy(frame, metric = gower, type = list(ordratio = c(1:36))) #k-medoid algorithm with 3 clusters kanswers - pam(daisy1, 3, diss = TRUE) #Evaluate k-mediod clustering algorithm with 2 to 6 clusters #Import clValid package library(clValid) #Internal validation internval1 - clValid(daisy1, 2:6, clMethods = pam, validation = internal) #Error in switch(class(obj), matrix = mat - obj, ExpressionSet = mat -Biobase::exprs(obj), : EXPR must be a length 1 vector #Error in summary(internval1) : #error in evaluating the argument 'object' in selecting a method for function 'summary': Error: object 'internval1' not found #External validation stabval1 - clValid(daisy1, 2:6, clMethods = pam, validation = stability) #Error in switch(class(obj), matrix = mat - obj, ExpressionSet = mat - Biobase::exprs(obj), : EXPR must be a length 1 vector Data: I put the data.frame in a dissimilarity matrix using the daisy function and used partitioning around medoids with 3 clusters. The daisy and pam functions come from the cluster CRAN package in R. Since the data.frame has mixed values, the gower distance coefficient is used. Here's the head of the first 7 variables, but I took out the names of the email for privacy reasons. head(frame) user_id emailAge Gender Household.Income Marital.Status Presence .of.children 1 12945 @bellycard.com NAMaleNANA NA 2 12947 @bellycard.com NAMaleNANA NA 3 12990 @gmail.com NANANANA NA 4 13160 @gmail.com 25-34 Male100k-125k Single No 5 13195 @gmail.com NAMale75k-100kSingle No 6 13286 @gmail.com NANANANA NA Please let me know if I can provide more information. -- Scott Davis Cell: (408)826-9561 Skype ID: Scdavis61 San Jose, CA. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Error creating daisy matrix in R cluster package - Cannot allocate vector of size 66.0 Gb
My purpose involves creating a dissimilarity matrix using the daisy package in R before applying k-mediod clustering for customer segmentation. The dataset has 133,153 observations of 35 variables in a data.frame with numerical, categorical, blank cells and missing values. Missing values refer to NA, while a blank cells means nothing present within the data.frame. Hereâs my OS: sessionInfo() R version 3.1.0 (2014-04-10) Platform x86_64-w64-mingw32/x64 (64-bit) I have 35 variables, but here is description of the first 5: head(df) user_idAgeGender Household.Income Marital.Status 1 12945 Male 2 12947 Male 3 12990 4 13160 25-34 Male 100k-125k Single 5 13195 Male 75k-100kSingle 6 13286 Since the Windows computer has 3 Gb RAM, I increased the virtual memory to 100Gb hoping that would be enough to create the matrix - it didn't work. I've looked into other R packages for solving the memory problem, but they don't work. I cannot use the `bigmemory` with the `biganalytics` package because it only accepts numeric matrices. The `clara` and `ff` packages also accept only numeric matrices. Here's the daisy script: #Load csv file Store1 - read.csv(/Users/name/Client1.csv, head = TRUE) #Convert csv to data.frame df -as.data.frame(Store1) #Increase memory allocation in R to 70 GB using the command: memory.limit(size = 7) [1] 7 #Load cluster package library(cluster) #Create daisy dissimilarity matrix #Use Gower distance coefficient for mixed variables #Set type as ratio scaled variable daisy1 - daisy(df, metric = gowerâ, type = list(ordratio = c(1:35))) #Error: cannot allocate vector of size 66.0 Gb How can I fix the error? -- Scott Davis Cell: (408)826-9561 Skype ID: Scdavis61 San Jose, CA. [[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-es] Bootstrap
2014-06-11 7:57 GMT-04:00 rubenfcasal rubenfca...@gmail.com: Hola Celia, Yo normalmente empleo el paquete boot (me parece mejor que el paquete bootstrap), pero en cualquier caso si necesitas información adicional sobre el tema, la referencia que se suele recomendar es: Davison, A.C. and Hinkley, D.V. (1997). Bootstrap Methods and Their Application. Cambridge University Press Un saludo, Rubén Fernández Casal El 11/06/2014 9:04, Celia Rubio Linares escribió: Hola! Tengo que hacer un proyecto acerca del paquete bootstrap de R, alguien podría facilitarme información completa (y en español a ser posible) acerca de este paquete? Muchas gracias de antemano, y un saludo. [[alternative HTML version deleted]] Hola Celia, Soy el encargado del paquete bootstrap. También recomiendo el paquete boot si es un trabajo serio (como otro paquete). Además, el paquete boot tiene una opción para que funcione en paralel. Todavía no hemos hecho lo mismo en el paquete boostrap, que es sobre todo para aprender. Y para este motivo creo que es mejor comenzar con el libro que inspiró el paquete, o sea: An Introduction to the Bootstrap por Bradley Efron. Este libro es una maravilla. Uno no tiene que tener much experiencia en matématicas para leerlo. Explica muy bien la intuición de por qué el bootstrap funciona (y por qué no en los casos en que no funciona). Busqué un poco pero creo que por desgracia el libro no se tardujo al español. Como puedes ver todos aquí te estamos recomiendo libros en vez de explicar el paquete bootstrap. Creo que es porque en cuanto a como usar el paquete, si entiendes el bootstrap, no hay mucha explicación necesaria. Ve los ejemplos. Cualquier duda, haznos una pregunta específica de lo que quieres hacer y lo que intentaste. Saludos, Scott ___ R-help-es mailing list R-help-es@r-project.org https://stat.ethz.ch/mailman/listinfo/r-help-es
Re: [R] Using apply with more than one matrix
I would ask you to look at this loop-free approach and ask if this is not equally valid? ans - matrix(NA, ncol=2, nrow=2) ind.not.na - which(!is.na(a1)) ans[] - condition1*a1[,,ind.not.na[1]]+ m2 # two matrices of equal dimensions, one logical. ans [,1] [,2] [1,] NA 1.66 [2,] 2.74 NA Thanks, I am learning something. I didn't know you could multiply a logical object by a numerical one. But notice the answer is not the same as mine, because I am doing an operation on the vector of values a1[i,j,] first. I tried a modification on sapply below, but it doesn't work because I haven't referenced the 3d array a1 properly. So I guess I must try to get a 2d result from a1 first, then use that in matrix arithmetic. Sapply or mapply may work, I haven't used these much and will try to learn better how to use them. Your use of sapply looks good; but I'm trying to understand if and how I can bring in the operation on a1. This doesn't work: evaluate - function(idx) { ind.not.na - which(!is.na(a1[idx,])) ])) # doesn't work; improper indexing for a1 if(length(ind.not.na) 0) { return(condition1*(a1[idx,ind.not.na[1]] + m2[idx])) # doesn't work; improper indexing for a1 } } vec - sapply(seq(length(m2)), evaluate) Scott Waichler -Original Message- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Wednesday, April 30, 2014 8:46 PM To: Waichler, Scott R Cc: Bert Gunter; r-help@r-project.org Subject: Re: [R] Using apply with more than one matrix On Apr 30, 2014, at 6:03 PM, Waichler, Scott R wrote: Here is a working example with no random parts. Thanks for your patience and if I'm still off the mark with my presentation I'll drop the matter. v - c(NA, 1.5, NA, NA, NA, 1.1, 0.5, NA, NA, 1.3, 0.4, 0.9) a1 - array(v, dim=c(2,2,3)) m1 - matrix(c(NA, 1.5, 2.1, NA), ncol=2, byrow=T) m2 - matrix(c(1.56, 1.64, 1.16, 2.92), ncol=2) condition1 - !is.na(m1) m1 m2 ans - matrix(NA, ncol=2, nrow=2) # initialize for(i in 1:2) { for(j in 1:2) { ind.not.na - which(!is.na(a1[i,j,])) if(condition1[i,j] length(ind.not.na) 0) ans[i,j] - a1[i,j,ind.not.na[1]] + m2[i,j] } } ans [,1] [,2] [1,] NA 1.66 [2,] 3.14 NA I would ask you to look at this loop-free approach and ask if this is not equally valid? ans - matrix(NA, ncol=2, nrow=2) ind.not.na - which(!is.na(a1)) ans[] - condition1*a1[,,ind.not.na[1]]+ m2 # two matrices of equal dimensions, one logical. ans [,1] [,2] [1,] NA 1.66 [2,] 2.74 NA Let me try asking again in words. If I have multiple matrices or slices of 3d arrays that are the same dimension, is there a way to pass them all to apply, and have apply take care of looping through i,j? I don't think `apply` is the correct function for this. Either `mapply` or basic matrix operation seem more likely to deliver correct results: I understand that apply has just one input object x. I want to work on more than one array object at once using a custom function that has this characteristic: in order to compute the answer at i,j I need a result from higher order array at the same i,j. If you want to iterate over matrix indices you can either use the vector version e.g. m2[3] or the matrix version, m2[2,1[. vec - sapply(seq(length(m2) , function(idx) m2[idx]*condition1[idx] ) This is what I tried to demonstrate in my example above. Thanks, Scott David Winsemius Alameda, CA, USA __ 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] Using apply with more than one matrix
Sapply or mapply may work, I haven't used these much and will try to learn better how to use them. Your use of sapply looks good; but I'm trying to understand if and how I can bring in the operation on a1. This doesn't work: evaluate - function(idx) { ind.not.na - which(!is.na(a1[idx,])) ])) # doesn't work; improper indexing for a1 if(length(ind.not.na) 0) { return(condition1*(a1[idx,ind.not.na[1]] + m2[idx])) # doesn't work; improper indexing for a1 } } vec - sapply(seq(length(m2)), evaluate) Are we to assume that the length of `which(!is.na(a1[idx,])) ]))` is guaranteed to equal the length of the two other matrices? If not then what sort of relationships should be assumed? ind.not.na is just a vector that could be any length. It is a selection on the vector a1[i,j,]. I want to get the first element of that selection. The key relationship is that the dimensions of the matrices and the first two dimensions of the 3d arrays are the same. That is, i and j have the same range for all of them. Scott __ 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] Using apply with more than one matrix
Thank you, A.K. I learned from both of your solutions. I find the one that uses alply easier to follow and understand intuitively. I guess I'll want to learn more about what plyr can do. I've been using R for years but hadn't pushed vectorization far enough in my code. Now my computing will be more efficient. Thanks also to David and the others who responded to my question. It all helped. --Scott Waichler -Original Message- From: arun [mailto:smartpink...@yahoo.com] Sent: Thursday, May 01, 2014 6:24 AM To: R. Help Cc: Waichler, Scott R Subject: Re: [R] Using apply with more than one matrix Hi, Sorry, a typo in the previous function: - if (condition1[i] !is.na(indx3)) { arr[x1][indx3] + m2[i] ## should be mat2[i] } else NA --- Also, you can try: library(plyr) evaluateNew - function(arr, mat1, mat2){ if (!all(dim(mat1) == dim(mat2))) { stop(both matrices should have equal dimensions) } indx1 - unlist(alply(arr, c(1,2), function(x) {x[!is.na(x)][1]})) condition1 - !is.na(mat1) mat1 mat2 ifelse(condition1, indx1+mat2, NA) } evaluateNew(a1, m1, m2) # [,1] [,2] #[1,] NA 1.66 #[2,] 3.14 NA A.K. On Thursday, May 1, 2014 5:49 AM, arun smartpink...@yahoo.com wrote: Hi, You may try: evaluate - function(arr, mat1, mat2) { if (!all(dim(mat1) == dim(mat2))) { stop(both matrices should have equal dimensions) } indx1 - as.matrix(do.call(expand.grid, lapply(dim(arr), sequence))) indx2 - paste0(indx1[, 1], indx1[, 2]) condition1 - !is.na(mat1) mat1 mat2 ans - sapply(seq_along(unique(indx2)), function(i) { x1 - indx1[indx2 %in% unique(indx2)[i], ] indx3 - which(!is.na(arr[x1]))[1] if (condition1[i] !is.na(indx3)) { arr[x1][indx3] + m2[i] } else NA }) dim(ans) - dim(mat1) ans } evaluate(a1,m1,m2) # [,1] [,2] #[1,] NA 1.66 #[2,] 3.14 NA A.K. On Thursday, May 1, 2014 2:53 AM, Waichler, Scott R scott.waich...@pnnl.gov wrote: I would ask you to look at this loop-free approach and ask if this is not equally valid? ans - matrix(NA, ncol=2, nrow=2) ind.not.na - which(!is.na(a1)) ans[] - condition1*a1[,,ind.not.na[1]]+ m2 # two matrices of equal dimensions, one logical. ans [,1] [,2] [1,] NA 1.66 [2,] 2.74 NA Thanks, I am learning something. I didn't know you could multiply a logical object by a numerical one. But notice the answer is not the same as mine, because I am doing an operation on the vector of values a1[i,j,] first. I tried a modification on sapply below, but it doesn't work because I haven't referenced the 3d array a1 properly. So I guess I must try to get a 2d result from a1 first, then use that in matrix arithmetic. Sapply or mapply may work, I haven't used these much and will try to learn better how to use them. Your use of sapply looks good; but I'm trying to understand if and how I can bring in the operation on a1. This doesn't work: evaluate - function(idx) { ind.not.na - which(!is.na(a1[idx,])) ])) # doesn't work; improper indexing for a1 if(length(ind.not.na) 0) { return(condition1*(a1[idx,ind.not.na[1]] + m2[idx])) # doesn't work; improper indexing for a1 } } vec - sapply(seq(length(m2)), evaluate) Scott Waichler -Original Message- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Wednesday, April 30, 2014 8:46 PM To: Waichler, Scott R Cc: Bert Gunter; r-help@r-project.org Subject: Re: [R] Using apply with more than one matrix On Apr 30, 2014, at 6:03 PM, Waichler, Scott R wrote: Here is a working example with no random parts. Thanks for your patience and if I'm still off the mark with my presentation I'll drop the matter. v - c(NA, 1.5, NA, NA, NA, 1.1, 0.5, NA, NA, 1.3, 0.4, 0.9) a1 - array(v, dim=c(2,2,3)) m1 - matrix(c(NA, 1.5, 2.1, NA), ncol=2, byrow=T) m2 - matrix(c(1.56, 1.64, 1.16, 2.92), ncol=2) condition1 - !is.na(m1) m1 m2 ans - matrix(NA, ncol=2, nrow=2) # initialize for(i in 1:2) { for(j in 1:2) { ind.not.na - which(!is.na(a1[i,j,])) if(condition1[i,j] length(ind.not.na) 0) ans[i,j] - a1[i,j,ind.not.na[1]] + m2[i,j] } } ans [,1] [,2] [1,] NA 1.66 [2,] 3.14 NA I would ask you to look at this loop-free approach and ask if this is not equally valid? ans - matrix(NA, ncol=2, nrow=2) ind.not.na - which(!is.na(a1)) ans[] - condition1*a1[,,ind.not.na[1]]+ m2 # two matrices of equal dimensions, one logical. ans [,1] [,2] [1,] NA 1.66 [2,] 2.74 NA Let me try asking again in words. If I have multiple matrices or slices of 3d arrays that are the same dimension, is there a way to pass them all to apply, and have apply take care of looping through i,j? I don't
[R] Using apply with more than one matrix
Hi, I want to apply a custom function to all the elements of one matrix. The function uses the value in the same i,j in a second matrix. How can I use apply() or similar function to avoid nested loops in i and j? Thanks, Scott Waichler Pacific Northwest National Laboratory Richland, WA, USA __ 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] Using apply with more than one matrix
Ok, here is a toy example. I want to use a custom function that depends on more than one matrix/array as input, and I can't figure out how to do that with apply. v - c(NA, 1.5, NA, NA, NA, 1.1, 0.5, NA, NA, 1.3, 0.4, 0.9) a1 - array(v, dim=c(2,2,3)) m1 - matrix(c(NA, 1.5, 2.1, NA), ncol=2, byrow=T) m2 - matrix(runif(n=4, min=1, max=3), ncol=2) condition1 - ifelse(!is.na(m1) m1 m2, T, F) ans - matrix(NA, ncol=2, nrow=2) # initialize for(i in 1:2) { for(j in 1:2) { ind.not.na - which(!is.na(a1[i,j,])) if(condition1[i,j] length(ind.not.na) 0) { ans[i,j] - a1[i,j,ind.not.na[1]] + m2[i,j] } } } Scott Waichler -Original Message- From: Bert Gunter [mailto:gunter.ber...@gene.com] Sent: Wednesday, April 30, 2014 12:18 PM To: Waichler, Scott R Cc: r-help@r-project.org Subject: Re: [R] Using apply with more than one matrix Scott: Your problem specification is rather vague: What do you mean by use? In general, matrices are merely vectors with a dim attribute, so if you can do what you want with them as vectors, then that will work for them as matrices. For example: m1- matrix(1:6, nr=3) m2 - matrix(11:16, nr=2) m2*m2 [,1] [,2] [,3] [1,] 121 169 225 [2,] 144 196 256 If this does not meet your needs, you will have to follow the posting guide and provide both code and a minimal reproducible example. Cheers, Bert Bert Gunter Genentech Nonclinical Biostatistics (650) 467-7374 Data is not information. Information is not knowledge. And knowledge is certainly not wisdom. H. Gilbert Welch On Wed, Apr 30, 2014 at 10:54 AM, Waichler, Scott R scott.waich...@pnnl.gov wrote: Hi, I want to apply a custom function to all the elements of one matrix. The function uses the value in the same i,j in a second matrix. How can I use apply() or similar function to avoid nested loops in i and j? Thanks, Scott Waichler Pacific Northwest National Laboratory Richland, WA, USA __ 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] Using apply with more than one matrix
Here is a working example with no random parts. Thanks for your patience and if I'm still off the mark with my presentation I'll drop the matter. v - c(NA, 1.5, NA, NA, NA, 1.1, 0.5, NA, NA, 1.3, 0.4, 0.9) a1 - array(v, dim=c(2,2,3)) m1 - matrix(c(NA, 1.5, 2.1, NA), ncol=2, byrow=T) m2 - matrix(c(1.56, 1.64, 1.16, 2.92), ncol=2) condition1 - !is.na(m1) m1 m2 ans - matrix(NA, ncol=2, nrow=2) # initialize for(i in 1:2) { for(j in 1:2) { ind.not.na - which(!is.na(a1[i,j,])) if(condition1[i,j] length(ind.not.na) 0) ans[i,j] - a1[i,j,ind.not.na[1]] + m2[i,j] } } ans [,1] [,2] [1,] NA 1.66 [2,] 3.14 NA Let me try asking again in words. If I have multiple matrices or slices of 3d arrays that are the same dimension, is there a way to pass them all to apply, and have apply take care of looping through i,j? I understand that apply has just one input object x. I want to work on more than one array object at once using a custom function that has this characteristic: in order to compute the answer at i,j I need a result from higher order array at the same i,j. This is what I tried to demonstrate in my example above. Thanks, Scott __ 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] problem with pip2d() from ptinpoly
Thank you for catching that, Boris. I'm surprised, given that there is no mention of sense of direction in the package documentation. Scott Waichler -Original Message- From: Boris Steipe [mailto:boris.ste...@utoronto.ca] Sent: Friday, April 18, 2014 3:49 PM To: Waichler, Scott R; r-help@r-project.org Subject: Re: [R] problem with pip2d() from ptinpoly Apparently it matters whether your polygon is defined clockwise or counterclockwise. A point outside your triangle is recognized ... q2 - matrix(c(594893.0,115435.0), ncol=2, byrow=T) pip2d(Vertices = verts, Queries = q2) [1] 1 ... and defining the triangle in counterclockwise sense gives the expected behaviour. v2 - matrix(c(594891,115309,594891,117201,59,117201), ncol=2, byrow=T) pip2d(Vertices = v2, Queries = query) [1] 1 Cheers, B. On 2014-04-18, at 6:00 PM, Waichler, Scott R wrote: Hi, pip2d() doesn't seem to work correctly for me. I have a plot of a triangle that a query point fits inside, but the point is defined as outside the polygon by pip2d. library(ptinpoly) verts - matrix(c(594891,115309,59,117201,594891,117201), ncol=2, byrow=T) query - matrix(c(594885.0,115435.0), ncol=2, byrow=T) pip2d(Vertices = verts, Queries = query) # result = -1 # contrary to -1 output of pip2d, plot shows point lies within triangle plot(c(594400, 595000), c(115000, 117500), type=n) polygon(verts, border=red) points(x=query[,1], y=query[,2], col=blue) Scott Waichler Pacific Northwest National Laboratory Richland, WA, USA __ 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] problem with pip2d() from ptinpoly
Hi, pip2d() doesn't seem to work correctly for me. I have a plot of a triangle that a query point fits inside, but the point is defined as outside the polygon by pip2d. library(ptinpoly) verts - matrix(c(594891,115309,59,117201,594891,117201), ncol=2, byrow=T) query - matrix(c(594885.0,115435.0), ncol=2, byrow=T) pip2d(Vertices = verts, Queries = query) # result = -1 # contrary to -1 output of pip2d, plot shows point lies within triangle plot(c(594400, 595000), c(115000, 117500), type=n) polygon(verts, border=red) points(x=query[,1], y=query[,2], col=blue) Scott Waichler Pacific Northwest National Laboratory Richland, WA, USA __ 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] Unexpected behaviour with seq() and %in%
Hi, During the course of putting together a function I came across some unexpected behaviour when using seq() and %in%. I am creating two numeric vectors using seq(), and then using %in% to find the values in one vector that are in the other. Sometimes all the values are found, but sometimes a value is not found. It seems that the behaviour is inconsistent and I don't understand why. Below are some reproducible examples: vec1 - seq(from = 0.1, to = 20, by = 0.1) vec2 - seq(from = 2, to = 20, by = 2) vec1[vec1 %in% vec2] [1] 2 4 6 8 10 12 14 16 20 The value 18 is not found. Starting vec1 at 0.2: vec1 - seq(from = 0.2, to = 20, by = 0.1) vec2 - seq(from = 2, to = 20, by = 2) vec1[vec1 %in% vec2] [1] 2 4 8 10 12 14 16 18 20 Now the value 6 is not found. Starting vec1 at 0.5: vec1 - seq(from = 0.5, to = 20, by = 0.1) vec2 - seq(from = 2, to = 20, by = 2) vec1[vec1 %in% vec2] [1] 2 4 6 8 10 12 14 16 18 20 Now all the values are found. Can someone please explain what is happening? Why should the starting value of vec1 make a difference? Thanks Finlay R version info: platform x86_64-pc-linux-gnu arch x86_64 os linux-gnu system x86_64, linux-gnu version.string R version 3.0.1 (2013-05-16) [[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] BH correction with p.adjust
My understanding was that the vector was ranked, the adjusted p vector was calculated and then the vector is returned to the original order - I came across a stack overflow answer saying this: http://stackoverflow.com/questions/10323817/r-unexpected-results-from-p-adjust-fdr Although the code there does not appear to be the same as when I type p.adjust into the command line. The order shouldn't matter anyway since my data was ordered by p. Yesterday I tried a short example of 5 numbers and it seemed to work out though today I tried to do another short example to demonstrate that the order in the p vector you input doesn't matter but didn't quite get a working example this time. Maybe due to a rounding to first significant figure or something? smallP - c(0.01, 0.5, 0.0001) names(smallP) - c(first, second, last) p.adjust(smallP) first second last 2e-02 5e-01 3e-04 0.01*3/2 [1] 0.015 0.5*3/3 [1] 0.5 0.0001*3/1 [1] 3e-04 In any case I reconstructed a large example which can be run without real data where the figure is way off and definitely not the result of a rounding error: exampleP - seq(from=0.001, to=0.1, by=0.0001) length(exampleP) [1] 991 examplePBH - p.adjust(exampleP, method=BH) exampleP[1] [1] 1e-07 examplePBH[1] [1] 0.1 exampleP[1]*length(exampleP)/1 [1] 0.991 Any help with this would be very much appreciated. It seems like it ought to be such a simple and commonly used method and yet I am struggling and not sure what to do about it. Thanks, Scott From: David Winsemius [dwinsem...@comcast.net] Sent: 21 July 2013 03:33 To: Scott Robinson Cc: r-help@r-project.org Subject: Re: [R] BH correction with p.adjust On Jul 20, 2013, at 10:37 AM, Scott Robinson wrote: Dear List, I have been trying to use p.adjust() to do BH multiple test correction and have gotten some unexpected results. I thought that the equation for this was: pBH = p*n/i Looking at the code for `p.adjust`, you see that the method is picked from a switch function lp - length(p) BH = { i - lp:1L o - order(p, decreasing = TRUE) ro - order(o) pmin(1, cummin(n/i * p[o]))[ro] } You may not have sorted the p-values in pList. where p is the original p value, n is the number of tests and i is the rank of the p value. However when I try and recreate the corrected p from my most significant value it does not match up to the one computed by the method p.adjust: setwd(C:/work/Methylation/IMA/GM/siteLists) hypTable - read.delim(hypernormal vs others.txt) pList - hypTable$p names(pList) - hypTable$site adjusted - p.adjust(pList, method=BH) adjusted[1] cg27433479 0.05030589 pList[1]*nrow(hypTable)/1 cg27433479 0.09269194 No data provided, so unable to pursue this further. I tried to recreate this is a small example of a vector of 5 p values but everything worked as expected there. I was wondering if there is some subtle difference about how p.adjust operates? Is there something more complicated about how to calculate 'n' or 'i' - perhaps due to identical p values being assigned the same rank or something? Does anyone have an idea what might be going on here? -- David Winsemius Alameda, CA, USA __ 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] BH correction with p.adjust
Dear List, I have been trying to use p.adjust() to do BH multiple test correction and have gotten some unexpected results. I thought that the equation for this was: pBH = p*n/i where p is the original p value, n is the number of tests and i is the rank of the p value. However when I try and recreate the corrected p from my most significant value it does not match up to the one computed by the method p.adjust: setwd(C:/work/Methylation/IMA/GM/siteLists) hypTable - read.delim(hypernormal vs others.txt) pList - hypTable$p names(pList) - hypTable$site adjusted - p.adjust(pList, method=BH) adjusted[1] cg27433479 0.05030589 pList[1]*nrow(hypTable)/1 cg27433479 0.09269194 I tried to recreate this is a small example of a vector of 5 p values but everything worked as expected there. I was wondering if there is some subtle difference about how p.adjust operates? Is there something more complicated about how to calculate 'n' or 'i' - perhaps due to identical p values being assigned the same rank or something? Does anyone have an idea what might be going on here? Many thanks, Scott __ 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] Survey imputation
I'm working with NHIS survye data. I'd like the to use muliple imputation to cover the missing data for the variables in which I'm interested. My question concerns the use of certain variables in the imputation model. For example, race would be an important predictor in the imputation model, but it has been imputed (hot deck) as have other variables that might be useful in the imputation model. What is the wisdom of using data that have been imputed in subsequent imputations models? __ 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] Survey imputation
I paln on using R to do the imputation once I figure out how to do it. - Original Message - From: Bert Gunter gunter.ber...@gene.com To: Scott Raynaud scott.rayn...@yahoo.com Cc: r-help@r-project.org r-help@r-project.org Sent: Thursday, June 13, 2013 12:45 PM Subject: Re: [R] Survey imputation Is this an R question? Seems like it belongs on a statistical or survey list, not r-help. Cheers, Bert On Thu, Jun 13, 2013 at 10:37 AM, Scott Raynaud scott.rayn...@yahoo.com wrote: I'm working with NHIS survye data. I'd like the to use muliple imputation to cover the missing data for the variables in which I'm interested. My question concerns the use of certain variables in the imputation model. For example, race would be an important predictor in the imputation model, but it has been imputed (hot deck) as have other variables that might be useful in the imputation model. What is the wisdom of using data that have been imputed in subsequent imputations models? __ 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. -- 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 __ 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] SPlus script
Ok. I tried copying the original program, changing all _ to - and running in 3.0.1. Still no error message or output. It's a mystery to me. - Original Message - From: William Dunlap wdun...@tibco.com To: Scott Raynaud scott.rayn...@yahoo.com; r-help@r-project.org r-help@r-project.org Cc: Sent: Wednesday, June 5, 2013 2:17 PM Subject: RE: [R] SPlus script Both the R and S+ versions (which seem to differ only in the use of _ for assignment in the S+ version) do nothing but define some functions. You would not expect any printed output unless you used those functions on some data. Is there another script that does that? Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Scott Raynaud Sent: Wednesday, June 05, 2013 6:21 AM To: r-help@r-project.org Subject: [R] SPlus script This originally was an SPlus script that I modifeid about a year-and-a-half ago. It worked perfectly then. Now I can't get any output despite not receiving an error message. I'm providing the SPLUS script as a reference. I'm running R15.2.2. Any help appreciated. MY MODIFICATION*** ** ## sshc.ssc: sample size calculation for historical control studies ## J. Jack Lee (jj...@mdanderson.org) and Chi-hong Tseng ## Department of Biostatistics, Univ. of Texas M.D. Anderson Cancer Center ## ## 3/1/99 ## updated 6/7/00: add loess ##-- Required Input: # # rc number of response in historical control group # nc sample size in historical control # d target improvement = Pe - Pc # method 1=method based on the randomized design # 2=Makuch Simon method (Makuch RW, Simon RM. Sample size considerations # for non-randomized comparative studies. J of Chron Dis 1980; 3:175-181. # 3=uniform power method optional Input: # # alpha size of the test # power desired power of the test # tol convergence criterion for methods 1 2 in terms of sample size # tol1 convergence criterion for method 3 at any given obs Rc in terms of difference # of expected power from target # tol2 overall convergence criterion for method 3 as the max absolute deviation # of expected power from target for all Rc # cc range of multiplicative constant applied to the initial values ne # l.span smoothing constant for loess # # Note: rc is required for methods 1 and 2 but not 3 # method 3 return the sample size need for rc=0 to (1-d)*nc # Output # for methdos 1 2: return the sample size needed for the experimental group (1 number) # for given rc, nc, d, alpha, and power # for method 3: return the profile of sample size needed for given nc, d, alpha, and power # vector $ne contains the sample size corresponding to rc=0, 1, 2, ... nc*(1-d) # vector $Ep contains the expected power corresponding to # the true pc = (0, 1, 2, ..., nc*(1-d)) / nc # #-- sshc-function(rc, nc=1092, d=.085779816, method=3, alpha=0.05, power=0.8, tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5) { ### for method 1 if (method==1) { ne1-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01) return(ne=ne1) } ### for method 2 if (method==2) { ne-nc ne1-nc+50 while(abs(ne-ne1)tol ne110){ ne-ne1 pe-d+rc/nc ne1-nef(rc,nc,pe*ne,ne,alpha,power) ## if(is.na(ne1)) print(paste('rc=',rc,',nc=',nc,',pe=',pe,',ne=',ne)) } if (ne110) return(NA) else return(ne=ne1) } ### for method 3 if (method==3) { if (tol1 tol2/10) tol1-tol2/10 ncstar-(1-d)*nc pc-(0:ncstar)/nc ne-rep(NA,ncstar + 1) for (i in (0:ncstar)) { ne[i+1]-ss.rand(i,nc,d,alpha=.05,power=.8,tol=.01) } plot(pc,ne,type='l',ylim=c(0,max(ne)*1.5)) ans-c.searchd(nc, d, ne, alpha, power, cc, tol1) ### check overall absolute deviance old.abs.dev-sum(abs(ans$Ep-power)) ##bad-0 print(round(ans$Ep,4)) print(round(ans$ne,2)) lines(pc,ans$ne,lty=1,col=8) old.ne-ans$ne ##while(max(abs(ans$Ep-power))tol2 bad==0){ unnecessary ## while(max(abs(ans$Ep-power))tol2){ ans-c.searchd(nc, d, ans$ne, alpha, power, cc, tol1) abs.dev-sum(abs(ans$Ep-power)) print(paste( old.abs.dev=,old.abs.dev)) print(paste( abs.dev=,abs.dev)) ##if (abs.dev old.abs.dev) { bad-1} old.abs.dev-abs.dev print(round(ans$Ep,4)) print(round(ans$ne,2)) lines(pc,old.ne,lty=1,col=1) lines(pc,ans$ne,lty=1,col=8) ### add convex ans$ne-convex(pc,ans$ne)$wy ### add loess ###old.ne-ans$ne loess.ne-loess(ans$ne ~ pc, span=l.span) lines(pc,loess.ne$fit,lty=1,col=4) old.ne-loess.ne$fit ###readline() } return
Re: [R] SPlus script
Ok. Now I see that the sshc function is not being called. Thanks for pointing that out. I'm not certain about the solution, however. I tried putting call(sshc) at the end of the program, but nothing happened. My memory about all of this is fuzzy. Suggestions on how to call the function appreciated. - Original Message - From: William Dunlap wdun...@tibco.com To: Scott Raynaud scott.rayn...@yahoo.com; r-help@r-project.org r-help@r-project.org Cc: Sent: Wednesday, June 5, 2013 2:17 PM Subject: RE: [R] SPlus script Both the R and S+ versions (which seem to differ only in the use of _ for assignment in the S+ version) do nothing but define some functions. You would not expect any printed output unless you used those functions on some data. Is there another script that does that? Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Scott Raynaud Sent: Wednesday, June 05, 2013 6:21 AM To: r-help@r-project.org Subject: [R] SPlus script This originally was an SPlus script that I modifeid about a year-and-a-half ago. It worked perfectly then. Now I can't get any output despite not receiving an error message. I'm providing the SPLUS script as a reference. I'm running R15.2.2. Any help appreciated. MY MODIFICATION*** ** ## sshc.ssc: sample size calculation for historical control studies ## J. Jack Lee (jj...@mdanderson.org) and Chi-hong Tseng ## Department of Biostatistics, Univ. of Texas M.D. Anderson Cancer Center ## ## 3/1/99 ## updated 6/7/00: add loess ##-- Required Input: # # rc number of response in historical control group # nc sample size in historical control # d target improvement = Pe - Pc # method 1=method based on the randomized design # 2=Makuch Simon method (Makuch RW, Simon RM. Sample size considerations # for non-randomized comparative studies. J of Chron Dis 1980; 3:175-181. # 3=uniform power method optional Input: # # alpha size of the test # power desired power of the test # tol convergence criterion for methods 1 2 in terms of sample size # tol1 convergence criterion for method 3 at any given obs Rc in terms of difference # of expected power from target # tol2 overall convergence criterion for method 3 as the max absolute deviation # of expected power from target for all Rc # cc range of multiplicative constant applied to the initial values ne # l.span smoothing constant for loess # # Note: rc is required for methods 1 and 2 but not 3 # method 3 return the sample size need for rc=0 to (1-d)*nc # Output # for methdos 1 2: return the sample size needed for the experimental group (1 number) # for given rc, nc, d, alpha, and power # for method 3: return the profile of sample size needed for given nc, d, alpha, and power # vector $ne contains the sample size corresponding to rc=0, 1, 2, ... nc*(1-d) # vector $Ep contains the expected power corresponding to # the true pc = (0, 1, 2, ..., nc*(1-d)) / nc # #-- sshc-function(rc, nc=1092, d=.085779816, method=3, alpha=0.05, power=0.8, tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5) { ### for method 1 if (method==1) { ne1-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01) return(ne=ne1) } ### for method 2 if (method==2) { ne-nc ne1-nc+50 while(abs(ne-ne1)tol ne110){ ne-ne1 pe-d+rc/nc ne1-nef(rc,nc,pe*ne,ne,alpha,power) ## if(is.na(ne1)) print(paste('rc=',rc,',nc=',nc,',pe=',pe,',ne=',ne)) } if (ne110) return(NA) else return(ne=ne1) } ### for method 3 if (method==3) { if (tol1 tol2/10) tol1-tol2/10 ncstar-(1-d)*nc pc-(0:ncstar)/nc ne-rep(NA,ncstar + 1) for (i in (0:ncstar)) { ne[i+1]-ss.rand(i,nc,d,alpha=.05,power=.8,tol=.01) } plot(pc,ne,type='l',ylim=c(0,max(ne)*1.5)) ans-c.searchd(nc, d, ne, alpha, power, cc, tol1) ### check overall absolute deviance old.abs.dev-sum(abs(ans$Ep-power)) ##bad-0 print(round(ans$Ep,4)) print(round(ans$ne,2)) lines(pc,ans$ne,lty=1,col=8) old.ne-ans$ne ##while(max(abs(ans$Ep-power))tol2 bad==0){ unnecessary ## while(max(abs(ans$Ep-power))tol2){ ans-c.searchd(nc, d, ans$ne, alpha, power, cc, tol1) abs.dev-sum(abs(ans$Ep-power)) print(paste( old.abs.dev=,old.abs.dev)) print(paste( abs.dev=,abs.dev)) ##if (abs.dev old.abs.dev) { bad-1} old.abs.dev-abs.dev print(round(ans$Ep,4)) print(round(ans$ne,2)) lines(pc,old.ne,lty=1,col=1) lines(pc,ans$ne,lty=1,col=8) ### add convex ans$ne-convex(pc
Re: [R] SPlus script
I actually had tried placing arguments in the call but it didn't work. However, I did not think about writing it to a variable and printing. That seems to have done the trick. Funny, I don't remember having to do that before, but that's not surprising. Anyway, thanks for helping to diagnose and fix the problem. - Original Message - From: Ista Zahn istaz...@gmail.com To: Scott Raynaud scott.rayn...@yahoo.com Cc: William Dunlap wdun...@tibco.com; r-help@r-project.org r-help@r-project.org Sent: Thursday, June 6, 2013 9:15 AM Subject: Re: [R] SPlus script Presumably something like r - sshc(50) print(r) But if you were getting output before than you already have a script that does something like this. It would be better to find it... Best, Ista On Thu, Jun 6, 2013 at 9:02 AM, Scott Raynaud scott.rayn...@yahoo.com wrote: Ok. Now I see that the sshc function is not being called. Thanks for pointing that out. I'm not certain about the solution, however. I tried putting call(sshc) at the end of the program, but nothing happened. My memory about all of this is fuzzy. Suggestions on how to call the function appreciated. - Original Message - From: William Dunlap wdun...@tibco.com To: Scott Raynaud scott.rayn...@yahoo.com; r-help@r-project.org r-help@r-project.org Cc: Sent: Wednesday, June 5, 2013 2:17 PM Subject: RE: [R] SPlus script Both the R and S+ versions (which seem to differ only in the use of _ for assignment in the S+ version) do nothing but define some functions. You would not expect any printed output unless you used those functions on some data. Is there another script that does that? Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Scott Raynaud Sent: Wednesday, June 05, 2013 6:21 AM To: r-help@r-project.org Subject: [R] SPlus script This originally was an SPlus script that I modifeid about a year-and-a-half ago. It worked perfectly then. Now I can't get any output despite not receiving an error message. I'm providing the SPLUS script as a reference. I'm running R15.2.2. Any help appreciated. MY MODIFICATION*** ** ## sshc.ssc: sample size calculation for historical control studies ## J. Jack Lee (jj...@mdanderson.org) and Chi-hong Tseng ## Department of Biostatistics, Univ. of Texas M.D. Anderson Cancer Center ## ## 3/1/99 ## updated 6/7/00: add loess ##-- Required Input: # # rc number of response in historical control group # nc sample size in historical control # d target improvement = Pe - Pc # method 1=method based on the randomized design # 2=Makuch Simon method (Makuch RW, Simon RM. Sample size considerations # for non-randomized comparative studies. J of Chron Dis 1980; 3:175-181. # 3=uniform power method optional Input: # # alpha size of the test # power desired power of the test # tol convergence criterion for methods 1 2 in terms of sample size # tol1 convergence criterion for method 3 at any given obs Rc in terms of difference # of expected power from target # tol2 overall convergence criterion for method 3 as the max absolute deviation # of expected power from target for all Rc # cc range of multiplicative constant applied to the initial values ne # l.span smoothing constant for loess # # Note: rc is required for methods 1 and 2 but not 3 # method 3 return the sample size need for rc=0 to (1-d)*nc # Output # for methdos 1 2: return the sample size needed for the experimental group (1 number) # for given rc, nc, d, alpha, and power # for method 3: return the profile of sample size needed for given nc, d, alpha, and power # vector $ne contains the sample size corresponding to rc=0, 1, 2, ... nc*(1-d) # vector $Ep contains the expected power corresponding to # the true pc = (0, 1, 2, ..., nc*(1-d)) / nc # #-- sshc-function(rc, nc=1092, d=.085779816, method=3, alpha=0.05, power=0.8, tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5) { ### for method 1 if (method==1) { ne1-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01) return(ne=ne1) } ### for method 2 if (method==2) { ne-nc ne1-nc+50 while(abs(ne-ne1)tol ne110){ ne-ne1 pe-d+rc/nc ne1-nef(rc,nc,pe*ne,ne,alpha,power) ## if(is.na(ne1)) print(paste('rc=',rc,',nc=',nc,',pe=',pe,',ne=',ne)) } if (ne110) return(NA) else return(ne=ne1) } ### for method 3 if (method==3) { if (tol1 tol2/10) tol1-tol2/10 ncstar-(1-d
[R] SPlus script
This originally was an SPlus script that I modifeid about a year-and-a-half ago. It worked perfectly then. Now I can't get any output despite not receiving an error message. I'm providing the SPLUS script as a reference. I'm running R15.2.2. Any help appreciated. MY MODIFICATION* ## sshc.ssc: sample size calculation for historical control studies ## J. Jack Lee (jj...@mdanderson.org) and Chi-hong Tseng ## Department of Biostatistics, Univ. of Texas M.D. Anderson Cancer Center ## ## 3/1/99 ## updated 6/7/00: add loess ##-- Required Input: # # rc number of response in historical control group # nc sample size in historical control # d target improvement = Pe - Pc # method 1=method based on the randomized design # 2=Makuch Simon method (Makuch RW, Simon RM. Sample size considerations # for non-randomized comparative studies. J of Chron Dis 1980; 3:175-181. # 3=uniform power method optional Input: # # alpha size of the test # power desired power of the test # tol convergence criterion for methods 1 2 in terms of sample size # tol1 convergence criterion for method 3 at any given obs Rc in terms of difference # of expected power from target # tol2 overall convergence criterion for method 3 as the max absolute deviation # of expected power from target for all Rc # cc range of multiplicative constant applied to the initial values ne # l.span smoothing constant for loess # # Note: rc is required for methods 1 and 2 but not 3 # method 3 return the sample size need for rc=0 to (1-d)*nc # Output # for methdos 1 2: return the sample size needed for the experimental group (1 number) # for given rc, nc, d, alpha, and power # for method 3: return the profile of sample size needed for given nc, d, alpha, and power # vector $ne contains the sample size corresponding to rc=0, 1, 2, ... nc*(1-d) # vector $Ep contains the expected power corresponding to # the true pc = (0, 1, 2, ..., nc*(1-d)) / nc # #-- sshc-function(rc, nc=1092, d=.085779816, method=3, alpha=0.05, power=0.8, tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5) { ### for method 1 if (method==1) { ne1-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01) return(ne=ne1) } ### for method 2 if (method==2) { ne-nc ne1-nc+50 while(abs(ne-ne1)tol ne110){ ne-ne1 pe-d+rc/nc ne1-nef(rc,nc,pe*ne,ne,alpha,power) ## if(is.na(ne1)) print(paste('rc=',rc,',nc=',nc,',pe=',pe,',ne=',ne)) } if (ne110) return(NA) else return(ne=ne1) } ### for method 3 if (method==3) { if (tol1 tol2/10) tol1-tol2/10 ncstar-(1-d)*nc pc-(0:ncstar)/nc ne-rep(NA,ncstar + 1) for (i in (0:ncstar)) { ne[i+1]-ss.rand(i,nc,d,alpha=.05,power=.8,tol=.01) } plot(pc,ne,type='l',ylim=c(0,max(ne)*1.5)) ans-c.searchd(nc, d, ne, alpha, power, cc, tol1) ### check overall absolute deviance old.abs.dev-sum(abs(ans$Ep-power)) ##bad-0 print(round(ans$Ep,4)) print(round(ans$ne,2)) lines(pc,ans$ne,lty=1,col=8) old.ne-ans$ne ##while(max(abs(ans$Ep-power))tol2 bad==0){ unnecessary ## while(max(abs(ans$Ep-power))tol2){ ans-c.searchd(nc, d, ans$ne, alpha, power, cc, tol1) abs.dev-sum(abs(ans$Ep-power)) print(paste( old.abs.dev=,old.abs.dev)) print(paste( abs.dev=,abs.dev)) ##if (abs.dev old.abs.dev) { bad-1} old.abs.dev-abs.dev print(round(ans$Ep,4)) print(round(ans$ne,2)) lines(pc,old.ne,lty=1,col=1) lines(pc,ans$ne,lty=1,col=8) ### add convex ans$ne-convex(pc,ans$ne)$wy ### add loess ###old.ne-ans$ne loess.ne-loess(ans$ne ~ pc, span=l.span) lines(pc,loess.ne$fit,lty=1,col=4) old.ne-loess.ne$fit ###readline() } return(list(ne=ans$ne, Ep=ans$Ep)) } } ## needed for method 1 nef2-function(rc,nc,re,ne,alpha,power){ za-qnorm(1-alpha) zb-qnorm(power) xe-asin(sqrt((re+0.375)/(ne+0.75))) xc-asin(sqrt((rc+0.375)/(nc+0.75))) ans- 1/(4*(xc-xe)^2/(za+zb)^2-1/(nc+0.5)) - 0.5 return(ans) } ## needed for method 2 nef-function(rc,nc,re,ne,alpha,power){ za-qnorm(1-alpha) zb-qnorm(power) xe-asin(sqrt((re+0.375)/(ne+0.75))) xc-asin(sqrt((rc+0.375)/(nc+0.75))) ans-(za*sqrt(1+(ne+0.5)/(nc+0.5))+zb)^2/(2*(xe-xc))^2-0.5 return(ans) } ## needed for method 3 c.searchd-function(nc, d, ne, alpha=0.05, power=0.8, cc=c(0.1,2),tol1=0.0001){ #--- # nc sample size of control group # d the differece to detect between control and experiment # ne vector of starting sample size of experiment group # corresonding to rc of 0 to nc*(1-d) # alpha size of test # power target power # cc pre-screen vector of constant c, the range should cover the # the value of cc
Re: [R] SPlus script
See my comments below. From: Pascal Oettli kri...@ymail.com To: Scott Raynaud scott.rayn...@yahoo.com Sent: Wednesday, June 5, 2013 10:02 AM Subject: Re: [R] SPlus script Hello, 1) It is always nice to say something as Hello, Tried, but could make it come out in anything other than four letters. 2) What do you want us to do with that script, without the required commented, minimal, self-contained, reproducible code? I want it to run as it did before. I cannot explain why the exact same code used to run and now doesn't. I've tried everything I can think of, so maybe brighter minds than mine can shed some light. 3) The lastest version of R is 3.0.1. Yes, I know. I was trying to keep things as constant as possible. Besides, do you really think the version wouold make a difference? The code runs without error. It just doesn't produce output. Regards, Pascal 2013/6/5 Scott Raynaud scott.rayn...@yahoo.com This originally was an SPlus script that I modifeid about a year-and-a-half ago. It worked perfectly then. Now I can't get any output despite not receiving an error message. I'm providing the SPLUS script as a reference. I'm running R15.2.2. Any help appreciated. MY MODIFICATION* ## sshc.ssc: sample size calculation for historical control studies ## J. Jack Lee (jj...@mdanderson.org) and Chi-hong Tseng ## Department of Biostatistics, Univ. of Texas M.D. Anderson Cancer Center ## ## 3/1/99 ## updated 6/7/00: add loess ##-- Required Input: # # rc number of response in historical control group # nc sample size in historical control # d target improvement = Pe - Pc # method 1=method based on the randomized design # 2=Makuch Simon method (Makuch RW, Simon RM. Sample size considerations # for non-randomized comparative studies. J of Chron Dis 1980; 3:175-181. # 3=uniform power method optional Input: # # alpha size of the test # power desired power of the test # tol convergence criterion for methods 1 2 in terms of sample size # tol1 convergence criterion for method 3 at any given obs Rc in terms of difference # of expected power from target # tol2 overall convergence criterion for method 3 as the max absolute deviation # of expected power from target for all Rc # cc range of multiplicative constant applied to the initial values ne # l.span smoothing constant for loess # # Note: rc is required for methods 1 and 2 but not 3 # method 3 return the sample size need for rc=0 to (1-d)*nc # Output # for methdos 1 2: return the sample size needed for the experimental group (1 number) # for given rc, nc, d, alpha, and power # for method 3: return the profile of sample size needed for given nc, d, alpha, and power # vector $ne contains the sample size corresponding to rc=0, 1, 2, ... nc*(1-d) # vector $Ep contains the expected power corresponding to # the true pc = (0, 1, 2, ..., nc*(1-d)) / nc # #-- sshc-function(rc, nc=1092, d=.085779816, method=3, alpha=0.05, power=0.8, tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5) { ### for method 1 if (method==1) { ne1-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01) return(ne=ne1) } ### for method 2 if (method==2) { ne-nc ne1-nc+50 while(abs(ne-ne1)tol ne110){ ne-ne1 pe-d+rc/nc ne1-nef(rc,nc,pe*ne,ne,alpha,power) ## if(is.na(ne1)) print(paste('rc=',rc,',nc=',nc,',pe=',pe,',ne=',ne)) } if (ne110) return(NA) else return(ne=ne1) } ### for method 3 if (method==3) { if (tol1 tol2/10) tol1-tol2/10 ncstar-(1-d)*nc pc-(0:ncstar)/nc ne-rep(NA,ncstar + 1) for (i in (0:ncstar)) { ne[i+1]-ss.rand(i,nc,d,alpha=.05,power=.8,tol=.01) } plot(pc,ne,type='l',ylim=c(0,max(ne)*1.5)) ans-c.searchd(nc, d, ne, alpha, power, cc, tol1) ### check overall absolute deviance old.abs.dev-sum(abs(ans$Ep-power)) ##bad-0 print(round(ans$Ep,4)) print(round(ans$ne,2)) lines(pc,ans$ne,lty=1,col=8) old.ne-ans$ne ##while(max(abs(ans$Ep-power))tol2 bad==0){ unnecessary ## while(max(abs(ans$Ep-power))tol2){ ans-c.searchd(nc, d, ans$ne, alpha, power, cc, tol1) abs.dev-sum(abs(ans$Ep-power)) print(paste( old.abs.dev=,old.abs.dev)) print(paste( abs.dev=,abs.dev)) ##if (abs.dev old.abs.dev) { bad-1} old.abs.dev-abs.dev print(round(ans$Ep,4)) print(round(ans$ne,2)) lines(pc,old.ne,lty=1,col=1) lines(pc,ans$ne,lty=1,col=8) ### add convex ans$ne-convex(pc,ans$ne)$wy ### add loess ###old.ne-ans$ne loess.ne-loess(ans$ne ~ pc, span=l.span) lines(pc,loess.ne$fit,lty=1,col=4) old.ne-loess.ne$fit ###readline() } return(list(ne=ans$ne, Ep=ans$Ep
Re: [R] SPlus script
As you can see, my original modification was to replace all _ with -. It worked after I did that. This is a simulation that generates its own data based on the given parameters. I should obtain a plot of the number of experimental subjects as a function of the repsonse rate in the historical group. Maybe I'm forgetting something here... It's been a while since I ran this. Thanks for your help. - Original Message - From: William Dunlap wdun...@tibco.com To: Scott Raynaud scott.rayn...@yahoo.com; r-help@r-project.org r-help@r-project.org Cc: Sent: Wednesday, June 5, 2013 2:17 PM Subject: RE: [R] SPlus script Both the R and S+ versions (which seem to differ only in the use of _ for assignment in the S+ version) do nothing but define some functions. You would not expect any printed output unless you used those functions on some data. Is there another script that does that? Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Scott Raynaud Sent: Wednesday, June 05, 2013 6:21 AM To: r-help@r-project.org Subject: [R] SPlus script This originally was an SPlus script that I modifeid about a year-and-a-half ago. It worked perfectly then. Now I can't get any output despite not receiving an error message. I'm providing the SPLUS script as a reference. I'm running R15.2.2. Any help appreciated. MY MODIFICATION*** ** ## sshc.ssc: sample size calculation for historical control studies ## J. Jack Lee (jj...@mdanderson.org) and Chi-hong Tseng ## Department of Biostatistics, Univ. of Texas M.D. Anderson Cancer Center ## ## 3/1/99 ## updated 6/7/00: add loess ##-- Required Input: # # rc number of response in historical control group # nc sample size in historical control # d target improvement = Pe - Pc # method 1=method based on the randomized design # 2=Makuch Simon method (Makuch RW, Simon RM. Sample size considerations # for non-randomized comparative studies. J of Chron Dis 1980; 3:175-181. # 3=uniform power method optional Input: # # alpha size of the test # power desired power of the test # tol convergence criterion for methods 1 2 in terms of sample size # tol1 convergence criterion for method 3 at any given obs Rc in terms of difference # of expected power from target # tol2 overall convergence criterion for method 3 as the max absolute deviation # of expected power from target for all Rc # cc range of multiplicative constant applied to the initial values ne # l.span smoothing constant for loess # # Note: rc is required for methods 1 and 2 but not 3 # method 3 return the sample size need for rc=0 to (1-d)*nc # Output # for methdos 1 2: return the sample size needed for the experimental group (1 number) # for given rc, nc, d, alpha, and power # for method 3: return the profile of sample size needed for given nc, d, alpha, and power # vector $ne contains the sample size corresponding to rc=0, 1, 2, ... nc*(1-d) # vector $Ep contains the expected power corresponding to # the true pc = (0, 1, 2, ..., nc*(1-d)) / nc # #-- sshc-function(rc, nc=1092, d=.085779816, method=3, alpha=0.05, power=0.8, tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5) { ### for method 1 if (method==1) { ne1-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01) return(ne=ne1) } ### for method 2 if (method==2) { ne-nc ne1-nc+50 while(abs(ne-ne1)tol ne110){ ne-ne1 pe-d+rc/nc ne1-nef(rc,nc,pe*ne,ne,alpha,power) ## if(is.na(ne1)) print(paste('rc=',rc,',nc=',nc,',pe=',pe,',ne=',ne)) } if (ne110) return(NA) else return(ne=ne1) } ### for method 3 if (method==3) { if (tol1 tol2/10) tol1-tol2/10 ncstar-(1-d)*nc pc-(0:ncstar)/nc ne-rep(NA,ncstar + 1) for (i in (0:ncstar)) { ne[i+1]-ss.rand(i,nc,d,alpha=.05,power=.8,tol=.01) } plot(pc,ne,type='l',ylim=c(0,max(ne)*1.5)) ans-c.searchd(nc, d, ne, alpha, power, cc, tol1) ### check overall absolute deviance old.abs.dev-sum(abs(ans$Ep-power)) ##bad-0 print(round(ans$Ep,4)) print(round(ans$ne,2)) lines(pc,ans$ne,lty=1,col=8) old.ne-ans$ne ##while(max(abs(ans$Ep-power))tol2 bad==0){ unnecessary ## while(max(abs(ans$Ep-power))tol2){ ans-c.searchd(nc, d, ans$ne, alpha, power, cc, tol1) abs.dev-sum(abs(ans$Ep-power)) print(paste( old.abs.dev=,old.abs.dev)) print(paste( abs.dev=,abs.dev)) ##if (abs.dev old.abs.dev) { bad-1} old.abs.dev-abs.dev print(round(ans$Ep,4)) print(round(ans$ne
[R] vector field from a 3D scalar field
I have a 3D field of a scalar variable (x, y, z, value). Is there a way to generate a vector field from this data--gradient at defined points? I found the rasterVis package for 2D data, but as yet nothing for 3D data. Thanks, Scott Waichler Pacific Northwest National Laboratory Richland, WA scott.waich...@pnnl.gov __ 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] how to get xmlToList() to retry if http fails
Hi, I am using xmlToList() in a loop with a call to a webservice, per the code below. # Loop thru target locs for(i in 1:num.target.locs) { url - paste(sep=/, http://www.earthtools.org/timezone;, lat[i], lon[i]) tmp - xmlToList(url) df$time.offset[i] - tmp$offset system(sleep 1) # wait 1 second per requirements of above web service } # end loop thru target locations Failure struck midway through my loop, with the message below. failed to load HTTP resource Error: 1: failed to load HTTP resource I presume that the webservice failed to respond in this instance. How can I trap the error and have it retry after waiting a second or two, instead of exiting? Thanks. --Scott Waichler Pacific Northwest National Laboratory Richland, WA, USA scott.waich...@pnnl.gov __ 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] Iterating through slots of an S4 object
I want to loop through slots of an S4 object and am unsure how to do so since the only way I can find to access them is individually in the form 'object@slotName'. I have guessed at a few possibilities which did not work and I have read some S4 object tutorials and things but still unsuccessful. I presume it is possible though? Any help would be much appreciated, Scott [[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] WriteXLS: 'object not found' error within function
Dear All, I am using WriteXLS to write tables with multiple sheets with the command: WriteXLS(tables, ExcelFileName = fileName, SheetNames = tableList, perl = perl, verbose = FALSE, Encoding = c(UTF-8, latin1), row.names = TRUE, col.names = TRUE, AdjWidth = TRUE, AutoFilter = FALSE, BoldHeaderRow = FALSE, FreezeRow = 0, FreezeCol = 0, envir = parent.frame()) ...where tables is the name of a list of data.frames which are the tables to go into my sheets. I am having no problem running my code unless it is within a function in which case I get the error: Error in get(as.character(x), envir = envir) : object 'tables' not found I have checked that there is not an issue with scope using these two lines within my function immediately before calling the WriteXLS function: print(class(tables)) flush.console() [1] list At least I would have thought this would have ruled out any scope-related issues. Has anyone else had this problem or have any ideas why it might be happening? Thanks, Scott PS here is the function in full: makeTables - function(design, designInfo, oldMatrix, matrix, annot, tableList) { rownames(design) - designInfo[,1] fit - lmFit(matrix, design) fit - eBayes(fit) tables-list() for (i in 1:length(tableList)) { table - makeTable(oldMatrix, matrix, annot, fit, tableList[i]) print(paste(tableList[i], =, sep=)) print(table) tables[[i]] - table } #Saving bit (set to FALSE if not using) saveIt - TRUE if(saveIt) { fileName - paste0(~,paste(colnames(design)[2:length(colnames(design))], collapse = +),.xls) print(paste(sheetnames=,length(tableList), print(class(tables)) flush.console() WriteXLS(tables, ExcelFileName = fileName, SheetNames = tableList, perl = perl, verbose = FALSE, Encoding = c(UTF-8, latin1), row.names = TRUE, col.names = TRUE, AdjWidth = TRUE, AutoFilter = FALSE, BoldHeaderRow = FALSE, FreezeRow = 0, FreezeCol = 0, envir = parent.frame()) } } __ 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] sweave and lessR
Hello I am using sweave with the Density() function from lessR package. However, unlike the older color.density() function, Density() does not work with the standard graphic output functions in R, such as pdf. Is there away to include figures fgenerated rom lessR in sweave documents? thank you __ 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] sweave and lessR
Thank you for the reply. With many calls to Density() I miss Sweave taking care of the naming of the files and the generating of the \includegraphics statements, but I guess this will work. On Dec 16, 2012, at 2:54 PM, Duncan Murdoch murdoch.dun...@gmail.com wrote: On 12-12-16 12:42 PM, Scott Cushman wrote: Hello I am using sweave with the Density() function from lessR package. However, unlike the older color.density() function, Density() does not work with the standard graphic output functions in R, such as pdf. Is there away to include figures fgenerated rom lessR in sweave documents? Tell Density to write to foo.pdf, then just use \includegraphics{foo} in your Sweave document. (This assumes you're using pdflatex; it looks as though you can't do it at all if you can't handle PDF figures.) Duncan Murdoch __ 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] analysis of bitmaps
I'm interested in using R to perform some statistical analysis of bitmaps. When I search, I see a lot of information about outputting bitmaps, but I'm not finding much about loading a bitmap into a data frame so that it can be analyzed. Would someone have a quick pointer to help me out? [[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] generating random samples of IG distribution
I don't read R-help these days so have just seen this. Both generalized inverse Gaussian and normal inverse Gaussian are in GeneralizedHyperbolic. HyperbolicDist is no longer being maintained. David Scott On 12/06/2012 5:41 a.m., David L Carlson wrote: Should have been For the normal inverse Gaussian: Package 'GeneralizedHyperbolic' For the generalized inverse Gaussian: Package 'HyperbolicDist' -- David L Carlson Associate Professor of Anthropology Texas AM University College Station, TX 77843-4352 -Original Message- From: David L Carlson [mailto:dcarl...@tamu.edu] Sent: Monday, June 11, 2012 10:26 AM To: 'shirin nezampour'; 'r-help@r-project.org' Subject: RE: [R] generating random samples of IG distribution For the normal inverse Gaussian: Package 'GeneralizedHyperbolic' For the generalized inverse Gaussian: Package 'GeneralizedHyperbolic' -- David L Carlson Associate Professor of Anthropology Texas AM University College Station, TX 77843-4352 -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- project.org] On Behalf Of shirin nezampour Sent: Sunday, June 10, 2012 11:37 AM To: r-help@r-project.org Subject: [R] generating random samples of IG distribution Dear R users, I want to generating random samples from Inverse Gaussian distribution . How can I do? and what package should I install? Thanks. Shirin [[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. -- _ David Scott Department of Statistics The University of Auckland, PB 92019 Auckland 1142,NEW ZEALAND Phone: +64 9 923 5055, or +64 9 373 7599 ext 85055 Email: d.sc...@auckland.ac.nz, Fax: +64 9 373 7018 __ 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
Here's the info Michael Weylandt requested: installed.packages()[c(lme4,nlme,Matrix),c(2,3,12)] LibPath Version Built lme4 /usr/lib/R/site-library 0.999375-40 2.13.1 nlme /usr/lib/R/library 3.1-104 2.15.0 Matrix /usr/lib/R/library 1.0-6 2.15.0 Looks like I have lme4 for version 2.13 rather than 2.15. Am I reading this right? I didn't install this and I'm not a Linux person so someone will have to tell me the command to obtain the Linux information requested. Then maybe I can figure whether to go to the Debian board or not. - Original Message - From: Duncan Murdoch murdoch.dun...@gmail.com To: Scott Raynaud scott.rayn...@yahoo.com Cc: r-help@r-project.org r-help@r-project.org Sent: Tuesday, August 28, 2012 10:40 AM Subject: Re: [R] Error message On 28/08/2012 10:31 AM, Scott Raynaud wrote: I suddenly started getting the error message below. Not sure why. If I type intalled.packages() it shows Matrix and lme4 installed. Can someone tell what's going on and what I need to do to remedy the problem? I'm running on a Linux box. Loading required package: Matrix Loading required package: lattice Error in dyn.load(file, DLLpath = DLLpath, ...) : function 'cholmod_l_start' not provided by package 'Matrix' Error: package/namespace load failed for ‘lme4’ Error: package/namespace load failed for ‘lme4’ I would guess you are using incompatible versions. Is R up to date? are both Matrix and lme4 up to date? (If you posted sessionInfo() we'd know this...) Duncan Murdoch __ 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 message
I suddenly started getting the error message below. Not sure why. If I type intalled.packages() it shows Matrix and lme4 installed. Can someone tell what's going on and what I need to do to remedy the problem? I'm running on a Linux box. Loading required package: Matrix Loading required package: lattice Error in dyn.load(file, DLLpath = DLLpath, ...) : function 'cholmod_l_start' not provided by package 'Matrix' Error: package/namespace load failed for ‘lme4’ Error: package/namespace load failed for ‘lme4’ __ 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
Here's the sessinoInfo(). I didn't do the upgrade and I don't know how to interpret the output. sessionInfo() R version 2.15.1 (2012-06-22) Platform: x86_64-pc-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=en_US.UTF-8 [9] LC_ADDRESS=en_US.UTF-8 LC_TELEPHONE=en_US.UTF-8 [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=en_US.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] Matrix_1.0-6 lattice_0.20-6 MASS_7.3-19 rkward_0.5.6 loaded via a namespace (and not attached): [1] grid_2.15.1 nlme_3.1-104 stats4_2.15.1 tools_2.15.1 - Original Message - From: Duncan Murdoch murdoch.dun...@gmail.com To: Scott Raynaud scott.rayn...@yahoo.com Cc: r-help@r-project.org r-help@r-project.org Sent: Tuesday, August 28, 2012 10:40 AM Subject: Re: [R] Error message On 28/08/2012 10:31 AM, Scott Raynaud wrote: I suddenly started getting the error message below. Not sure why. If I type intalled.packages() it shows Matrix and lme4 installed. Can someone tell what's going on and what I need to do to remedy the problem? I'm running on a Linux box. Loading required package: Matrix Loading required package: lattice Error in dyn.load(file, DLLpath = DLLpath, ...) : function 'cholmod_l_start' not provided by package 'Matrix' Error: package/namespace load failed for ‘lme4’ Error: package/namespace load failed for ‘lme4’ I would guess you are using incompatible versions. Is R up to date? are both Matrix and lme4 up to date? (If you posted sessionInfo() we'd know this...) Duncan Murdoch __ 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] RGDAL OGRwrite question
I have a quick question: It appears that in rgdal v0.7-12 (R version 2.15.1, OSX 10.6.8) writeOGR will not write a shapefile the the current directory. Is this correct? An earlier version of rgdal must have allowed this because I have a older script that used to work, but doesn't now. So, as an example, here is what I get today: shape = readOGR('.', layer='S20_G75_V00_HAASHP10_R00') OGR data source with driver: ESRI Shapefile Source: ., layer: S20_G75_V00_HAASHP10_R00 with 169 features and 23 fields Feature type: wkbPolygon with 2 dimensions writeOGR(shape, '.', layer='temp', driver='ESRI Shapefile', verbose=TRUE) Error in writeOGR(shape, ., layer = temp, driver = ESRI Shapefile) : Creation of output file failed writeOGR(shape, '/tmp', layer='temp', driver='ESRI Shapefile',verbose=TRUE) $object_type [1] SpatialPolygonsDataFrame $output_dsn [1] /tmp $output_layer [1] temp $output_diver [1] ESRI Shapefile $output_n [1] 169 $output_nfields [1] 23 $output_fields [1] ID ANID F_AREA Avg_z Manning [6] IniWL IniSal DispCoeff HydRad veg1_DW [11] veg2_IWveg3_SWveg4_DCveg5_ICveg6_SC [16] veg7_Marsh veg8_Swamp Rain_StID WetlandOpenwater [21] iniBed_m maxH_m BoxID $output_fclasses [1] 0 4 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2 2 2 4 $dataset_options NULL $layer_options NULL -- Scott M. Duke-Sylvester Assistant Professor Department of Biology Office : 300 E. St. Mary Blvd Billeaud Hall, Room 141 Lafayette, LA 70504 Mailing address : UL Lafayette Department of Biology P.O.Box 42451 Lafayette, LA 70504-2451 Phone : 337 482 5304 Fax : 337 482 5834 email : smd3...@louisiana.edu This e-mail message (including any attachments) is for t...{{dropped:11}} __ 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] RGDAL OGRwrite question
In the example I gave, the temp.shp, 'temp.bdf, etc ... do not exist in the current working directory (.) before I call: writeOGR(shape, '.', layer='temp', driver='ESRI Shapefile', verbose=TRUE) I should have specified this. Also, the working directory permissions are set to owner read/write/execute (unix: drwxr-xr-x). On Fri, Aug 17, 2012 at 2:17 PM, MacQueen, Don macque...@llnl.gov wrote: Are the output files already there? In that case, try overwrite_layer=TRUE -Don -- Don MacQueen Lawrence Livermore National Laboratory 7000 East Ave., L-627 Livermore, CA 94550 925-423-1062 On 8/17/12 9:28 AM, Scott Duke-Sylvester scottdukesylves...@gmail.com wrote: I have a quick question: It appears that in rgdal v0.7-12 (R version 2.15.1, OSX 10.6.8) writeOGR will not write a shapefile the the current directory. Is this correct? An earlier version of rgdal must have allowed this because I have a older script that used to work, but doesn't now. So, as an example, here is what I get today: shape = readOGR('.', layer='S20_G75_V00_HAASHP10_R00') OGR data source with driver: ESRI Shapefile Source: ., layer: S20_G75_V00_HAASHP10_R00 with 169 features and 23 fields Feature type: wkbPolygon with 2 dimensions writeOGR(shape, '.', layer='temp', driver='ESRI Shapefile', verbose=TRUE) Error in writeOGR(shape, ., layer = temp, driver = ESRI Shapefile) : Creation of output file failed writeOGR(shape, '/tmp', layer='temp', driver='ESRI Shapefile',verbose=TRUE) $object_type [1] SpatialPolygonsDataFrame $output_dsn [1] /tmp $output_layer [1] temp $output_diver [1] ESRI Shapefile $output_n [1] 169 $output_nfields [1] 23 $output_fields [1] ID ANID F_AREA Avg_z Manning [6] IniWL IniSal DispCoeff HydRad veg1_DW [11] veg2_IWveg3_SWveg4_DCveg5_ICveg6_SC [16] veg7_Marsh veg8_Swamp Rain_StID WetlandOpenwater [21] iniBed_m maxH_m BoxID $output_fclasses [1] 0 4 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2 2 2 4 $dataset_options NULL $layer_options NULL -- Scott M. Duke-Sylvester Assistant Professor Department of Biology Office : 300 E. St. Mary Blvd Billeaud Hall, Room 141 Lafayette, LA 70504 Mailing address : UL Lafayette Department of Biology P.O.Box 42451 Lafayette, LA 70504-2451 Phone : 337 482 5304 Fax : 337 482 5834 email : smd3...@louisiana.edu This e-mail message (including any attachments) is for t...{{dropped:11}} __ 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. -- Scott M. Duke-Sylvester Assistant Professor Department of Biology Office : 300 E. St. Mary Blvd Billeaud Hall, Room 141 Lafayette, LA 70504 Mailing address : UL Lafayette Department of Biology P.O.Box 42451 Lafayette, LA 70504-2451 Phone : 337 482 5304 Fax : 337 482 5834 email : smd3...@louisiana.edu This e-mail message (including any attachments) is for the sole use of the intended recipient(s) and may contain confidential and privileged information. If the reader of this message is not the intended recipient, you are hereby notified that any dissemination, distribution or copying of this message (including any attachments) is strictly prohibited. If you have received this message in error, please contact the sender by reply e-mail message and destroy all copies of the original message (including attachments). __ 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] Accents and special character using hwriter (on Windows)
The safest way to include special characters is to use a character code. You are at the mercy of the browser otherwise and browsers behave very differently when confronted with unusual things. You don't mention which browser you are using on Windows, which is a severe gap in the information you provided. Here is an example which has been tested on Windows using IE, Firefox, Safari and Chrome (recent versions of all, I can't be bothered checking the versions). It also renders correctly if the file is opened in Word (versions based on xml). require(hwriter) pg - openPage(specialcharacters.html) hwrite(Test Special Characters, pg, heading = 1, br = TRUE) hwrite(Ciencias Sociales y Juriacute;dicas n:74 | 33.94%, pg, br = TRUE) hwrite(Ciencias Sociales y Jur#237;dicas n:74 | 33.94%, pg, br = TRUE) closePage(pg) For other special character codes, see for example http://www.ascii.cl/htmlcodes.htm David Scott On 31/07/2012 9:21 p.m., ramonovelar wrote: Thanks Arun, Yes, I have Windows 7. I have tried 2 versions of R, 2.14.1 and 2.15.x, but it did not change anything. Right now I can't try a different version of win. Ramón On Tuesday, July 31, 2012, arun kirshna [via R] wrote: Hello, I tried your code in R 2.15 with Ubuntu 12.04. It looks okay to me. datosdv-Ciencias Sociales y JurÃdicas n:74 | 33.94% print(datosdv) #[1] Ciencias Sociales y JurÃdicas n:74 | 33.94% library(hwriter) p=openPage('test.html') hwrite(datosdv,p,br=TRUE) #test.html output Ciencias Sociales y JurÃdicas n:74 | 33.94% Probably, it must be specific with the windows. Are you using windows7? A.K. - Original Message - From: ramonovelar[hidden email]http://user/SendEmail.jtp?type=nodenode=4638501i=0 To: [hidden email]http://user/SendEmail.jtp?type=nodenode=4638501i=1 Cc: Sent: Monday, July 30, 2012 7:11 PM Subject: [R] Accents and special character using hwriter (on Windows) Hello, I have a problem with special characters such as à or ñ when using hwriter. This only happens when I use windows, it works fine on mac. If I do: print(datosdv) Ciencias Sociales y JurÃdicas n:74 | 33.94% but: hwrite(datosdv, p, br=TRUE) Ciencias Sociales y Jur�dicas n:74 | 33.94% The bad sign is in the code, is not a problem of the encoding of the html page, that is in UTF-8. Does anybody have found this? Many thanks in advance. Ramón -- View this message in context: http://r.789695.n4.nabble.com/Accents-and-special-character-using-hwriter-on-Windows-tp4638474.html Sent from the R help mailing list archive at Nabble.com. __ [hidden email]http://user/SendEmail.jtp?type=nodenode=4638501i=2mailing 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. __ [hidden email]http://user/SendEmail.jtp?type=nodenode=4638501i=3mailing 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/Accents-and-special-character-using-hwriter-on-Windows-tp4638474p4638501.html To unsubscribe from Accents and special character using hwriter (on Windows), click herehttp://r.789695.n4.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_codenode=4638474code=cmFtb24ub3ZlbGFyQGdtYWlsLmNvbXw0NjM4NDc0fC0xNzk0Mjk1MDc3 . NAMLhttp://r.789695.n4.nabble.com/template/NamlServlet.jtp?macro=macro_viewerid=instant_html%21nabble%3Aemail.namlbase=nabble.naml.namespaces.BasicNamespace-nabble.view.web.template.NabbleNamespace-nabble.view.web.template.NodeNamespacebreadcrumbs=notify_subscribers%21nabble%3Aemail.naml-instant_emails%21nabble%3Aemail.naml-send_instant_email%21nabble%3Aemail.naml __ 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. -- _ David Scott Department of Statistics The University of Auckland, PB 92019 Auckland 1142,NEW ZEALAND Phone: +64 9 923 5055, or +64 9 373 7599 ext 85055 Email: d.sc...@auckland.ac.nz, Fax: +64 9 373 7018 [[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
Re: [R] Can't get R to recognize Java for rJava installation
Milan, Merci. I did find the javah file and put it in /usr/bin, where R can now find it. However, I still get a similar error message when trying to install rJava, i.e. configure: error: One or more Java configuration variables are not set. The only field that doesn't have a value now are the cpp flags: cpp flags : '' Could this be the problem now? How can I set those, and what value should I give? Scott Waichler So I guess you need to find out what package provides this file on your distribution (which you did not mention). First check the file is currently not present. __ 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] Can't read a binary file
Hi, I've read up on readBin() and chapter 6 in the R Data Import/Export manual, but I still can't read a binary file. Here is how the creator of the file described the code that would be needed in Fortran: Every record has a return in fortran. The length of each record is nx*ny*4. To read you would use the following: nlayx = nx*ny*4 do iz=1,nz,4 read(binary file) var(1:nlayx) enddo nrest=mod(nx*ny*nz,nlayx) read(binary file) var(1:nrest) The first value in the file should be 0.05, and all of the data values are real. Here is what I get (with similar answers using double): v-readBin(plotb.251, numeric(), size=4, n=1) v [1] 1.614296e-39 v-readBin(plotb.251, numeric(), size=4, n=1, endian=swap) v [1] 1.359775e-38 Platform is Intel Linux. How can I read the file described above? Thanks, Scott Waichler, PhD Hydrology Group, Energy Environment Directorate Pacific Northwest National Laboratory scott.waich...@pnnl.gov __ 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] Can't get R to recognize Java for rJava installation
Hi, I am unable to install the package rJava. I tried doing what the output suggests, but it doesn't help. How can I get R to find/recognize my Java installation? I am running R-2.15.0. waichler@snow sudo R CMD javareconf Java interpreter : /usr/bin/java Java version : 1.6.0_22 Java home path : /usr/lib/jvm/java-1.6.0-openjdk-1.6.0.0/jre Java compiler: /usr/bin/javac Java headers gen.: Java archive tool: /usr/bin/jar Java library path: $(JAVA_HOME)/lib/i386/client:$(JAVA_HOME)/lib/i386:$(JAVA_HOME)/../lib/i386:/usr/java/packages/lib/i386:/lib:/usr/lib JNI linker flags : -L$(JAVA_HOME)/lib/i386/client -L$(JAVA_HOME)/lib/i386 -L$(JAVA_HOME)/../lib/i386 -L/usr/java/packages/lib/i386 -L/lib -L/usr/lib -ljvm JNI cpp flags: Updating Java configuration in /usr/lib/R Done. install.packages(c(rJava), dependencies = T, repos = http://cran.fhcrc.org;) . . . checking whether siglongjmp is declared... yes checking Java support in R... present: interpreter : '/usr/bin/java' archiver: '/usr/bin/jar' compiler: '/usr/bin/javac' header prep.: '' cpp flags : '' java libs : '-L/usr/lib/jvm/java-1.6.0-openjdk-1.6.0.0/jre/lib/i386/client -L/usr/lib/jvm/java-1.6.0-openjdk-1.6.0.0/jre/lib/i386 -L/usr/lib/jvm/java-1.6.0-openjdk-1.6.0.0/jre/../lib/i386 -L/usr/java/packages/lib/i386 -L/lib -L/usr/lib -ljvm' configure: error: One or more Java configuration variables are not set. Make sure R is configured with full Java support (including JDK). Run R CMD javareconf as root to add Java support to R. If you don't have root privileges, run R CMD javareconf -e to set all Java-related variables and then install rJava. ERROR: configuration failed for package ÆrJavaÇ * removing Æ/usr/lib/R/library/rJavaÇ Thanks, Scott Waichler Senior Research Scientist Hydrology Group, Energy Environment Directorate Pacific Northwest National Laboratory __ 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] Apt-get
Ok. When I type in sudo add-apt-repository ppa:marutter/rrutter I get a prompt for a password. I enter my domain password and get some print about adding the ppa and requeting that I press enter. I press enter and this is followed by print stating that the key B04C661B is being requested from hkp sevrer keyserverubuntu.com followed by connection refused HTTP fetch error 7: couldn't connect no valid Open GPG found. I'm not sure what next but I think I need to configure the keyserver. Tried this via the command line and it failed. Maybe creating a key.txt file will work but when I search http://keyserver.ubuntu.com:11371/ for E084DAB9 that fails as well. What next? From: Tyler Ritchie tyler.ritc...@gmail.com @r-project.org Sent: Thursday, March 15, 2012 1:36 PM Subject: Re: [R] Apt-get Beltrand was also on the mark, suggesting you add Michael Rutter's ppa to your repository sources. In both cases (adding the CRAN Ubuntu repositories or Michael Rutter's ppa), an additional package repository is added to your system's packages. apt then checks that repository along with the other Ubuntu repositories and exposes the relevant binary R packages that Michael Rutter is curating. So, your IS people are correct in saying that the latest version of R available through the Ubuntu packages is 2.13.1, but there are more up to date R repositories to use. If you don't have the system rights to add any additional repositories, you'll need your IS folks to add them for you. -Tyler [[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] Apt-get
My IS people insist that the latest version of R avaialble via apt-get is 2.13.1. Anything later they claim will have to be compiled. True? Will I have to compile every time I update R? Seems like a lot of work. Surely there's a way around it. - Original Message - From: Jeff Newmiller jdnew...@dcn.davis.ca.us To: Scott Raynaud scott.rayn...@yahoo.com; r-help@r-project.org r-help@r-project.org Cc: Sent: Wednesday, March 7, 2012 12:34 PM Subject: Re: [R] Apt-get Google is really useful for questions like this. http://cran.r-project.org/bin/linux/ubuntu/ --- Jeff Newmiller The . . Go Live... DCN:jdnew...@dcn.davis.ca.us 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. Scott Raynaud scott.rayn...@yahoo.com wrote: I have a box set up with Kubuntu as the OS. I didn't perform the R install but was told the version of R available via the apt-get command was 2.13.1. Is there any way to get 2.14.0 in that same manner? __ 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] Apt-get
I have a box set up with Kubuntu as the OS. I didn't perform the R install but was told the version of R available via the apt-get command was 2.13.1. Is there any way to get 2.14.0 in that same manner? __ 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] Access to OpenBLAS
My IT people have set up a Kubuntu box with an RWkard front end. I have OpenBLAS set up as a shared BLAS but I'm not sure how to get R to see it. A.3.1 of the installation docs talks about it but I'm not clear if I need a option on my startup line or if I need to find a config file. The BLAS is is in: /usr/lib/openblas-base on my machine. Recommendations? __ 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] Simulation confidence interval
The follwing is a code snippet from a power simulation program that I'm using: estbeta-fixef(fitmodel) sdebeta-sqrt(diag(vcov(fitmodel))) for(l in 1:betasize) { cibeta-estbeta[l]-sgnbeta[l]*z1score*sdebeta[l] if(beta[l]*cibeta0) powaprox[[l]]-powaprox[[l]]+1 sdepower[l,iter]-as.numeric(sdebeta[l]) } Estbeta recovers the fixed effects from a model fitted using lmer. Beta is defined elsewhere and is a user specified input that relates the data generated in the simulation to an oucome. So, it seems pretty clear that the third line from the bottom is a clever test of whether the confidence interval traps 0. My question is why use beta[l]*cibeta0 rather than estbeta[l]*cibeta0. Is that because in the long run the model parameter etimates tend toward the betas specified by the user? In other words, what really matters is the standard errors, right? __ 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] Recompile
Looked at several posts and the installtions docs and still not clear. If I compile source codes and then somewhere down the line add a new package, then I have to recompile my entire installation, correct? Seems like this is the sentiment of the emails I read. __ 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] GUI
I'm setting up a Linux box to run R. I ususally run in a Windows envrionment but after reading the docs I'm not sure what to expect in terms of the front end appearance in Linux. Does it resemble Windows or will I need Rkward or R Commander? __ 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] BLAS
I'm setting up an Ubuntu virtual machine that will use 4-Intel Xeon CPU x5650. I'd like to compile R with a BLAS but the question is whcih one. Seems like the only free ones are GotoBLAS which I'm not sure is being maintained for newer CPUs and OpenBLAS for Loongson CPUs. I saw a favorable report on OpenBLAS (http://www.rochester.edu/college/gradstudents/jolmsted/files/computing/BLAS_Comparison.pdf), but I'm not sure it's the right thing for my CPUs. The webpage for OpenBLAS says, On X86 box, compile this library for loongson3a CPU. Any opinions on whether this will work? If not, any suggestions on another free BLAS? __ 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] grplasso
So does anyone use this package? - Original Message - From: Scott Raynaud scott.rayn...@yahoo.com To: r-help@r-project.org r-help@r-project.org Cc: Sent: Tuesday, January 10, 2012 1:40 PM Subject: grplasso I want to use the grplasso package on a data set where I want to fit a linear model. My interest is in identifying significant beta coefficients. The documentation is a bit cryptic so I'd appreciate some help. I know this is a strategy for large numbers of variables but consider a simple case for pedagogical puposes. Say I have two 3 category predictors (2 dummies each), a binary predictor and a continuous predictor with a continuous outcome: y x1 x2 x3 x4 x5 x6 rows of data here .. .. Naturally, I want to select x1 and x2 as a group and x3 and x4 as another group. The documentation has a couple of examples but it's not clear how they translate to the current problem. How do I specify my groups and run the lasso regression? Looks like this is the grouping part: index-c(NA,) but I'm not sure how to specify the df for the variables past the NA for the intercept. Once that's defined the penalty can be specified: lambda - lambdamax(x, y = y, index = index, penscale = sqrt, model = LogReg()) * 0.5^(0:5) In my case I'd use LinReg for the model. Then the model: fit - grplasso(x, y = y, index = index, lambda = lambda, model = LogReg(), penscale = sqrt, control = grpl.control(update.hess = lambda, trace = 0)) again using LinReg for the model. This can be plotted against lambda, but when I do lasso regression in other software I end up with a plot of the coefficients against the tuning parameter with a cutpoint or a table and graph that tells me what to include in the model based on some selected criterion. It's not clear from the example if there's a cross-validation or some other procedure to determine what variables to include. Plot(fit) produces a graph of coefficients against lambda but nothig to indicate what to include. What is used in the package, if anything, to make that determination? __ 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] grplasso
I want to use the grplasso package on a data set where I want to fit a linear model. My interest is in identifying significant beta coefficients. The documentation is a bit cryptic so I'd appreciate some help. I know this is a strategy for large numbers of variables but consider a simple case for pedagogical puposes. Say I have two 3 category predictors (2 dummies each), a binary predictor and a continuous predictor with a continuous outcome: y x1 x2 x3 x4 x5 x6 rows of data here .. .. Naturally, I want to select x1 and x2 as a group and x3 and x4 as another group. The documentation has a couple of examples but it's not clear how they translate to the current problem. How do I specify my groups and run the lasso regression? Looks like this is the grouping part: index-c(NA,) but I'm not sure how to specify the df for the variables past the NA for the intercept. Once that's defined the penalty can be specified: lambda - lambdamax(x, y = y, index = index, penscale = sqrt, model = LogReg()) * 0.5^(0:5) In my case I'd use LinReg for the model. Then the model: fit - grplasso(x, y = y, index = index, lambda = lambda, model = LogReg(), penscale = sqrt, control = grpl.control(update.hess = lambda, trace = 0)) again using LinReg for the model. This can be plotted against lambda, but when I do lasso regression in other software I end up with a plot of the coefficients against the tuning parameter with a cutpoint or a table and graph that tells me what to include in the model based on some selected criterion. It's not clear from the example if there's a cross-validation or some other procedure to determine what variables to include. Plot(fit) produces a graph of coefficients against lambda but nothig to indicate what to include. What is used in the package, if anything, to make that determination? __ 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] R report generator (for Word)?
The html route is one I have used quite a lot, but rather than R2HTML I far prefer hwriter. I have spent some time on enhancing hwriter and you can find my hwriterPlus on R-forge. It has fairly extensive examples and a vignette in the inst directory. I am still working on some improvements to the package. David Scott From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] on behalf of Joshua Wiley [jwiley.ps...@gmail.com] Sent: Monday, January 02, 2012 9:31 AM To: Michael Cc: r-help Subject: Re: [R] R report generator (for Word)? Hi Michael, I like Sweave and LaTeX, but I can appreciate the difficulty using it with collaborators. What about something similar using HTML? Certainly integrates to any webpages nicely. There are two packages I think do this nicely, one is the R2HTML package (on CRAN). Another one that is not on CRAN yet, but I think has a lot of potential is the knitr package. You can find it on github. I am not personally familiar with any good ways to integrate R with MS Office products. Cheers, Josh On Sun, Jan 1, 2012 at 7:50 AM, Michael comtech@gmail.com wrote: Happy New Year all! I am looking for a good solution for keeping record of my experiments - could you please help me? My work is about analysing data... My current work-flow: 1. Everyday my bosses give me some small steps/tasks for analysing data - which are parts of one bigger/whole project. 2. Everyday I send tens of emails to bosses/colleagues to report my findings in each step. 3. Bosses/colleagues often respond to my findings in real-time and suggest new experiments/steps and ask what-if questions. 4. I often have to manually copy and paste the results from R console and put them into an Excel and decorate a bit and send out. 5. Every one week and 2 weeks, we need to present to more senior bosses with more nice-looking presentations which is a summary of our findings in those 1-2 weeks. It's this time that is most chaotic because my colleagues and I have to dig into all the hundreds of emails in the past 1-2 weeks and copy and paste and organize those data again and make a nice overall summary for presentation... 6. As I am a hard-working guy, I myself often run my own random/ad-hoc experiments using out-of-work time and whenever I have interesting findings, I will send to immediate bosses and colleagues to seek their comments. 7. All these experiments are in fact variations of different versions/ideas of one big/whole project. Lets say in one big project bosses/colleagues and I have come up with a few big ideas, then we have a few sub-projects: MyProjectIdea1 MyProjectIdea2 ... MyProjectIdeaN And each idea has a few variations, mostly are for answering what-if questions by varying the parameters here and there ... For example: MyProjectIdea1_Variation1_WhatIfParam1ChangedTo1.2? ... ... etc. 8. Most experiments run tens of minutes to many hours... and some of them have to run on Linux, and some others can be run on Windows. Fortunately we have universal paths accessible on both Windows and Linux, so those won't be problem... 9. Because of the time-consuming nature of these experiments, I also save the images as rData whenever I can. However, it's necessary to keep track of the context where these data were generated. Otherwise even the records of these images won't help recall the scenario we have run... --- Keeping track of these changes and all kinds of what-ifs now becomes increasingly a problem for me. Some times in order to respond to a query, although I have done it before already, but because I didn't keep record and save the result, or even though I have saved the memory image yet I am not completely sure about the cleanness of the results/data,I have to redo it and wait for another few hours. Is there a way that I can manage these whole processes better and be more productive? I have been digging and thinking about this for while and I guess Sweave is the right way to go? The problem for Sweave is that it's hard to make Latex generated pdf appealing to business managers... so if I keep records in Sweave/Latex for my own record/benefit (that's already a big benefit)... I still need to somehow manually copy/paste the data from Sweave/Latex/pdf into Word/Excel/Powerpoint in order to make a nice presentation... I know there are some Open Office and Word version of Sweave... the problem is that I couldn't find many demonstrations on these topics and my question is: are they good and can they fulfill what we needed? Your thoughts are greatly appreciated! Thanks a lot! [[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
Re: [R] rbinom
This makes sense. Guess I should have put a pencil to it. Further investigation revealed that it is indeed a possibility that the relation between x and y is nonlinear: ax+bx^2+c where a, b and c are to be determined. My question is how to code this in my simulated data. I could do something like this after appropriately defining beta. meanpred and varpred: x[,2]-rnorm(length,meanpred[2],sqrt(varpred[2])) x[,3]-rnorm(length,meanpred[3],sqrt(varpred[3])) fixpart-x%*%beta binomprob-exp(fixpart)/(1+exp(fixpart)) data$y-rbinom(n1,1,binomprob) but I'd need to square my x[,3] values before multiplying them by beta. Can I say: x[,3]-(rnorm(length,meanpred[3],sqrt(varpred[3])))^2 in lieu of x[,3]-rnorm(length,meanpred[3],sqrt(varpred[3]))? - Original Message - From: peter dalgaard pda...@gmail.com To: Scott Raynaud scott.rayn...@yahoo.com Cc: r-help@r-project.org r-help@r-project.org Sent: Tuesday, December 27, 2011 9:15 AM Subject: Re: [R] rbinom On Dec 27, 2011, at 15:47 , Scott Raynaud wrote: I have the following code (which I did not write) that generates data based on a logistic model. I'm only getting a single record with y=1. It seems implausible that in 50k cases that have a single y=1. Does that ring alarm bells for anyone else? Not really. As far as I can tell, fixpart is roughly -10.5 (= -1.5 - .25*36), so binomprob is around 2.75e-5, which - nonlinearity notwithstanding - suggests that the expected number of positives out of 50K is something like 1.4. To do this more precisely, just compute and print sum(binomprob) in the code you gave. beta-c(-1.585600,-0.246900) betasize-length(beta) meanpred-c(0,35.90) varpred-c(0,1.00) #loop code x-matrix(1,length,betasize) #length set to 50k #loop code x[,2]-rnorm(length,meanpred[2],sqrt(varpred[2])) #length set to 50k fixpart-x%*%beta binomprob-exp(fixpart)/(1+exp(fixpart)) data$y-rbinom(n1,1,binomprob) #more loop 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. -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd@cbs.dk Priv: pda...@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.
[R] rbinom
I have the following code (which I did not write) that generates data based on a logistic model. I'm only getting a single record with y=1. It seems implausible that in 50k cases that have a single y=1. Does that ring alarm bells for anyone else? beta-c(-1.585600,-0.246900) betasize-length(beta) meanpred-c(0,35.90) varpred-c(0,1.00) #loop code x-matrix(1,length,betasize) #length set to 50k #loop code x[,2]-rnorm(length,meanpred[2],sqrt(varpred[2])) #length set to 50k fixpart-x%*%beta binomprob-exp(fixpart)/(1+exp(fixpart)) data$y-rbinom(n1,1,binomprob) #more loop 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.