Re: [R] Substituting Greek symbols in some tick labels
On Jul 4, 2013, at 8:14 PM, Eric Archer - NOAA Federal wrote: I have a character vector that I'm using to label ticks in a dotchart. Some of the elements in the vector have an asterisk (*) where a Greek Delta needs to be placed when the plot is generated. Here's a simple example: x - 1:4 x.lab - c(a*a, bbb, c*c, ddd) dotchart(x, labels = x.lab) The first and third labels should be 'aDeltaa' and 'cDeltac'. I've tried things like, x.lab - strsplit(x.lab, [*]) x.lab - lapply(x.lab, function(y) expression(paste(y, sep = Delta))) The plotmath function paste has no sep argument. Do you want to do this by hand? (Since you have not offered values of 'y'.) x.lab - expression( a*Delta*a, bbb, c*Delta*c, ddd) # Note use of * and no quotes in an expression vector. x - 1:4 dotchart(x, labels = x.lab) -- David. but because 'y' is unevaluated, the resulting list elements won't work as tick labels. I've tried to modify it by using bquote and substitute, but couldn't get anything closer. Any suggestions? Thanks! Cheers, eric -- Eric Archer, Ph.D. Southwest Fisheries Science Center NMFS, NOAA 8901 La Jolla Shores Drive La Jolla, CA 92037 USA 858-546-7121 (work) 858-546-7003 (FAX) Marine Mammal Genetics Group: swfsc.noaa.gov/prd-mmgenetics ETP Cetacean Assessment Program: swfsc.noaa.gov/prd-etp The universe doesn't care what you believe. The wonderful thing about science is that it doesn't ask for your faith, it just asks for your eyes. - Randall Munroe Lighthouses are more helpful than churches. - Benjamin Franklin ...but I'll take a GPS over either one. - John C. Craig George [[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. David Winsemius, MD 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] Substituting Greek symbols in some tick labels
David, thanks for the reply. No, I don't want to do it by hand because the actual data that I need to label is much larger and variable. The deltas occur in the label character vector for some labels and not others and are in different positions in the character strings. In the example I gave, the 'y' is the argument in the function to the lapply where I was trying to iterate over the results from 'strsplit' to insert the 'Delta' - the results showed me it wasn't the way to go, but I was giving an example of what I'd tried. Cheers, eric On Thu, Jul 4, 2013 at 11:31 PM, David Winsemius dwinsem...@comcast.netwrote: On Jul 4, 2013, at 8:14 PM, Eric Archer - NOAA Federal wrote: I have a character vector that I'm using to label ticks in a dotchart. Some of the elements in the vector have an asterisk (*) where a Greek Delta needs to be placed when the plot is generated. Here's a simple example: x - 1:4 x.lab - c(a*a, bbb, c*c, ddd) dotchart(x, labels = x.lab) The first and third labels should be 'aDeltaa' and 'cDeltac'. I've tried things like, x.lab - strsplit(x.lab, [*]) x.lab - lapply(x.lab, function(y) expression(paste(y, sep = Delta))) The plotmath function paste has no sep argument. Do you want to do this by hand? (Since you have not offered values of 'y'.) x.lab - expression( a*Delta*a, bbb, c*Delta*c, ddd) # Note use of * and no quotes in an expression vector. x - 1:4 dotchart(x, labels = x.lab) -- David. but because 'y' is unevaluated, the resulting list elements won't work as tick labels. I've tried to modify it by using bquote and substitute, but couldn't get anything closer. Any suggestions? Thanks! Cheers, eric -- Eric Archer, Ph.D. Southwest Fisheries Science Center NMFS, NOAA 8901 La Jolla Shores Drive La Jolla, CA 92037 USA 858-546-7121 (work) 858-546-7003 (FAX) Marine Mammal Genetics Group: swfsc.noaa.gov/prd-mmgenetics ETP Cetacean Assessment Program: swfsc.noaa.gov/prd-etp The universe doesn't care what you believe. The wonderful thing about science is that it doesn't ask for your faith, it just asks for your eyes. - Randall Munroe Lighthouses are more helpful than churches. - Benjamin Franklin ...but I'll take a GPS over either one. - John C. Craig George [[alternative HTML version deleted]] __** R-help@r-project.org mailing list https://stat.ethz.ch/mailman/**listinfo/r-helphttps://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/** posting-guide.html http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius, MD Alameda, CA, USA -- Eric Archer, Ph.D. Southwest Fisheries Science Center NMFS, NOAA 8901 La Jolla Shores Drive La Jolla, CA 92037 USA 858-546-7121 (work) 858-546-7003 (FAX) Marine Mammal Genetics Group: swfsc.noaa.gov/prd-mmgenetics ETP Cetacean Assessment Program: swfsc.noaa.gov/prd-etp The universe doesn't care what you believe. The wonderful thing about science is that it doesn't ask for your faith, it just asks for your eyes. - Randall Munroe Lighthouses are more helpful than churches. - Benjamin Franklin ...but I'll take a GPS over either one. - John C. Craig George [[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] g2 test...
Fernando Marmolejo Ramos fernando.marmolejoramos at adelaide.edu.au writes: [snip] is it appropriate to use a Log likelihood ratio (G-test) test of independence when dealing with repeated categorical responses (e.g. 2 by 2 table) instead of the McNemar test? [snip] Not an R question. Try http://stats.stackexchange.com instead? __ 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] g2 test...
On Jul 5, 2013, at 08:49 , Ben Bolker wrote: Fernando Marmolejo Ramos fernando.marmolejoramos at adelaide.edu.au writes: [snip] is it appropriate to use a Log likelihood ratio (G-test) test of independence when dealing with repeated categorical responses (e.g. 2 by 2 table) instead of the McNemar test? [snip] Not an R question. Try http://stats.stackexchange.com instead? It probably has a no-homework rule too, though -- 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] Error when building a custom package
Hi I'm have trouble building a custom package in R. Building the package on my colleagues computer (who made the R code content) has been (and is still) working fine for a long time, but I simply can't duplicate the setup right somehow (even with his help) I have installed the same versions of software as my colleague: Windows 7, R 2.13.1 (with the required additional packages), Rtools 2.13, Miktex 2.9. All of which are included in the Windows PATH variable as instructed by him (and in various online tutorials) I only get an error message, when trying to build a binary file (or runnning /R CMD check package name/), not the standard tar.gz. Here is what the output looks like in the command prompt: / * installing *source* package 'slo' ... ** R ** inst ** preparing package for lazy loading Loading required package: car Warning: package 'car' was built under R version 2.13.2 Loading required package: MASS Loading required package: nnet Loading required package: survival Loading required package: splines Creating a new generic function for plot in slo Error in file(file, r, encoding = encoding) : cannot open the connection ERROR: lazy loading failed for package 'slo' * removing 'C:/Users/CVO/workspace/slo.Rcheck/slo'/ These are the lines which follow, when my colleague run the build successfully: /Creating a new generic function for plot in slo in .GlobalEnv ** help ... / I have tried following a tutorial to see if the trouble was with building packages in general, but that does not seem to be the case, as I could get the same outputs as described in this tutorial: http://stevemosher.wordpress.com/ten-steps-to-building-an-r-package-under-windows/ http://stevemosher.wordpress.com/ten-steps-to-building-an-r-package-under-windows/ (which builds a simple package using the skeleton package). I have gotten the impression that generally the error with the Error in file... is a problem with a wrong path or problems with permissions, but am unable to find a solution. I'm running the commands in the Command prompt with Administrator rights and have Admin right on the computer. Any suggestions on how I can solve this issue or simply get more detailed information about the error will be greatly appreciated. Thanks, Christer V., DK -- View this message in context: http://r.789695.n4.nabble.com/Error-when-building-a-custom-package-tp4670866.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] problem on reading many files
Hi, I have to run almost 120 stations files of temperatura (mx and min), separated, and rain. I am following the instructions on RHtests tutorial as on http://www.cmc.org.ve/mediawiki/index.php?title=Preparando_los_datos link. Bit I have no success on running on multiple files. I have the .ls file which include the names of the file station, which I include on this msg for your consideration. The instruction that I wrote on the console is listatmn - readLines (tmaxcopia.ls, warn=false) for(ifile in listatmn) FindU(paste(./,ifile,sep=),paste(salidas/,ifile,sep=),-999.9) Of course the tmaxcopia.ls is the list of stations for tmax data. But once running on the console, this is the error msg I got listatmn - readLines (tmaxcopia.ls, warn=false) Error: inesperado entrada in listatmn - readLines ( for(ifile in listatmn) FindU(paste(./,ifile,sep=),paste(salidas/,ifile,sep=),-999.9) Error: objeto 'ifnames' no encontrado What I'd like to do is to make easier and fast running the RHtest Thank's a lot in advance for your help. Dr Victor M Rodríguez M Doctor en Ciencias, en Ciencias de la Tierra / Geociencias Ambientales INIFAP. CEPAB. (465) 9580167 y 86 Ext 108 y 220 / Este msg y su(s) atado(s) tiene carácter de CONFIDENCIAL y solo es para los fines descritos en su contenido. Queda expresamente prohibida bajo la ley de protección de datos del Gobierno de los Estados Unidos Mexicanos y para su difusion o extensión a terceros, se requiere del consentimiento expresso del titular de la cuenta de correo electrónico y/o del administrador del sitio en extenso de sus responsabilidades . / [[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] coding variables which are independent or grouped
Hi everyone, I have observations of the habitats of birds that I consider independent. I also have observations of habitats from several forest fragments with several observations within each fragment that are non-independent. I'm curious to know the best strategy to code these observations for analysis (JAGS or WinBUGS through R) Using R 3.0.1 in Windows. Here's what I have so far: code indicates the individual fragment (four digits) and the bird species (one digit and nine species), type indicates fragment size (in hectares) and species (now two digits). species just lets me know what species are which (taken out later) and leafno is the average number of dead leaves in a plot. Let me know if there's anything I can clarify. Thanks so much, Jeff structure(list(code = c(1112L, 1112L, 2107L, 2107L, 2107L, 2107L, 2108L, 2108L, 2108L, 2108L, 3114L, 3114L, 3114L, 3114L, 3114L, 3114L, 4101L, 4101L, 4101L, 4101L, 1212L, 1212L, 1212L, 1212L, 1212L, 1212L, 1212L, 1212L, 2206L, 2206L, 2206L, 2206L, 2206L, 2206L, 2206L, 2206L, 3209L, 3209L, 3209L, 3209L, 3209L, 3209L, 3209L, 3209L, 4201L, 4201L, 4201L, 4201L, 4201L, 4201L, 4201L, 4201L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 2303L, 2303L, 2303L, 2303L, 2303L, 2303L, 2303L, 2303L, 2303L, 2303L, 2303L, 2303L, 2303L, 2303L, 2303L, 3304L, 3304L, 3304L, 3304L, 3304L, 3304L, 3304L, 3304L, 3304L, 3304L, 3304L, 3304L, 3304L, 3304L, 3304L, 3304L), type = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 91L, 91L, 91L, 91L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 93L, 93L, 93L, 93L, 93L, 93L, 93L, 93L, 93L, 93L, 93L, 93L, 93L, 93L, 93L, 93L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 95L, 95L, 95L, 95L, 96L, 96L, 96L, 97L, 97L, 97L, 97L, 97L, 97L, 97L, 97L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 98L, 99L, 99L, 99L, 99L, 99L, 99L, 99L, 99L, 99L, 99L, 99L, 99L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L), species = structure(c(11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c(CA, CT, FA, FC, GV, HFR, HM, MC, MF, MT, OFR, TFR), class = factor), leafno = c(1.8, 2.65, 1.25, 2.6, 1.55, 2.35, 2.3, 1.7, 2.2, 3.2, 3.1, 1.25, 2.35, 1.9, 1.75, 2.8, 1.95, 2.45, 1.55, 1.65, 0.85, 1, 1.85, 1.45, 2.55, 4.05, 3.4, 2.1, 2.05, 1.95, 2.75, 2.6, 1.95, 2.15, 2.2, 2.1, 1.7, 2.2, 1.65, 2.2, 2.95, 2, 2.05, 1.45, 1.55, 0.7, 1.45, 2.05, 1.6, 1.4, 1.25, 2.25, 2.35, 2.4, 2.15, 1.65, 1.8, 2.25, 1.45, 2.15, 2.05, 1.5, 2, 1.7, 1.15, 1.25, 1, 1.15, 1.95, 1.05, 2.4, 1.1, 1.75, 1.7, 1.7, 1.75, 1.1, 1.05, 1.8, 1.45, 2.05, 1.9, 1.15, 2, 1.15, 2.9, 1.6, 2.65, 2.1, 2.35, 1.25, 1.8, 0.85, 1.6, 2.35, 1.75, 1.5, 1.75, 1.3, 1.8, 1.05, 1.45, 0.85, 0.9, 1.7, 2.05, 1.55, 1.55, 1.4, 1.8, 1.5, 1.2, 1.3, 2, 1.4, 2.6, 2.05, 1.05, 2.55, 2.4, 0.4, 2.1, 1.65, 1.65, 1.5, 1.5, 1.55, 2.4, 1.8, 2.2, 1.7, 1.05, 1.45, 1.5, 1.1, 0.95, 0.6, 1.3, 1.95, 1.3, 1.9, 1.55, 1.75, 1.65, 2.35, 3.2, 1.55, 1.65, 1.05, 2.6, 1.75, 1.6, 2.35, 1.2, 1.35, 2.5, 0.95, 2.15, 2.15, 2.55, 1.6, 2.1, 1.9, 0.95, 2.1, 2.4, 2.4, 1.7, 2.45, 1.85, 1.15, 1.25, 0.65, 1.35, 1.8, 1.4, 1.15, 2.75, 2.25, 2, 1.75, 1.4, 2.25)), .Names = c(code, type, species, leafno), class = data.frame, row.names = c(NA, -183L)) -- Jeffrey A. Stratford, PhD Department of Biology and Health Sciences Wilkes University, PA 18766 USA 570-332-2942 http://web.wilkes.edu/jeffrey.stratford/
Re: [R] Subsetting multiple rows of a data frame at once
Hi Anika,
?merge() is a better solution.
To get the row.names intact, you could do:
carbon.fit- within(carbon.fit,{x-round(x,10);y- round(y,10)}) #Using Bill's
solution
dat1- data.frame(x=round(xt,10),y=round(yt,10))
carbon.fit1-
data.frame(carbon.fit,rNames=row.names(carbon.fit),stringsAsFactors=FALSE)
#changed here
res1-merge(dat1,carbon.fit1,by=c(x,y))
row.names(res1)- res1[,3]
res1- res1[,-3]
A.K.
- Original Message -
From: William Dunlap wdun...@tibco.com
To: arun smartpink...@yahoo.com; Shaun ♥ Anika pro_pa...@hotmail.com
Cc: R help r-help@r-project.org
Sent: Thursday, July 4, 2013 8:02 PM
Subject: RE: [R] Subsetting multiple rows of a data frame at once
xt- c(1.05, 2.85, 3.40, 4.25, 0.25, 3.05, 3.70, 0.20, 0.30, 0.70, 1.05,
1.20, 1.40, 1.90,
2.70, 3.25, 3.55, 4.60, 2.05, 2.15, 3.70, 4.85, 4.90, 1.60, 2.45, 3.20, 3.90,
4.45)
yt- c(0.25, 0.10, 0.90, 0.25, 1.05, 1.70, 2.05, 2.90, 2.35, 2.60, 2.55,
2.15, 2.75, 2.05,
2.70, 2.25, 2.55, 2.05, 3.65, 3.05, 3.00, 3.50, 3.75, 4.85, 4.50, 4.50, 3.35,
4.90)
carbon.fit = expand.grid(list(x=seq(0, 5, 0.01), y=seq(0, 5, 0.01)))
trees-do.call(rbind,lapply(seq_along(xt),function(i)
subset(carbon.fit,x==xt[i]y==yt[i])))
## xt is 28 integers long and when i run the above code it only returns the
values of 18
out of the 28 (xt,yt) pairs that i want.
You are running into the problem that two different computational methods that
give
the same result when applied to real numbers often give different results when
applied
to 64-bit floating point numbers. (In your case you expect seq(0,5,.01) to
contain, e.g.,
the floating point number generate by parsing the string 3.05.) Hence x==y
is not true
when you expect it to be. Here is where your 18 came from:
R table(xt %in% carbon.fit$x, yt %in% carbon.fit$y)
FALSE TRUE
FALSE 1 6
TRUE 3 18
Round your number to the nearest 10^-10 and you get
table(round(xt,10) %in% round(carbon.fit$x,10), round(yt,10) %in%
round(carbon.fit$y,10))
TRUE
TRUE 28
By the way, you may prefer using the merge() function rather than the
do.call(rbind,lapply(...)))
business. I think the following call to merge will do about what you want (the
row names differ -
if they are important it is possible to get them with some minor trickery):
merge(data.frame(x=xt,y=yt), carbon.fit)
(You still want to round your numbers as before.)
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 arun
Sent: Wednesday, July 03, 2013 10:15 PM
To: Shaun ♥ Anika
Cc: R help
Subject: Re: [R] Subsetting multiple rows of a data frame at once
Hi,
carbon.fit = expand.grid(list(x=seq(0, 5, 0.01), y=seq(0, 5, 0.01)))
dim(carbon.fit)
#[1] 251001 2
xtNew-sprintf(%.2f,xt)
ytNew- sprintf(%.2f,yt)
carbon.fit[]- lapply(carbon.fit,function(x) sprintf(%.2f,x))
res-do.call(rbind,lapply(seq_along(xtNew),function(i)
subset(carbon.fit,x==xtNew[i]y==ytNew[i])))
nrow(res)
#[1] 28
res
# x y
#12631 1.05 0.25
#5296 2.85 0.10
#45431 3.40 0.90
#12951 4.25 0.25
#52631 0.25 1.05
#85476 3.05 1.70
#103076 3.70 2.05
#145311 0.20 2.90
#117766 0.30 2.35
#130331 0.70 2.60
#127861 1.05 2.55
#107836 1.20 2.15
#137916 1.40 2.75
#102896 1.90 2.05
#135541 2.70 2.70
#113051 3.25 2.25
#128111 3.55 2.55
#103166 4.60 2.05
#183071 2.05 3.65
#153021 2.15 3.05
#150671 3.70 3.00
#175836 4.85 3.50
#188366 4.90 3.75
#243146 1.60 4.85
#225696 2.45 4.50
#225771 3.20 4.50
#168226 3.90 3.35
#245936 4.45 4.90
A.K.
From: Shaun ♥ Anika pro_pa...@hotmail.com
To: smartpink...@yahoo.com smartpink...@yahoo.com
Sent: Thursday, July 4, 2013 12:08 AM
Subject: RE: Subsetting multiple rows of a data frame at once
Hi There,
i can give you the data needed to perform this task...
library(akima)
library(fields)
xt- c(1.05, 2.85, 3.40, 4.25, 0.25, 3.05, 3.70, 0.20, 0.30, 0.70, 1.05,
1.20, 1.40, 1.90,
2.70, 3.25, 3.55, 4.60, 2.05, 2.15, 3.70, 4.85, 4.90, 1.60, 2.45, 3.20, 3.90,
4.45)
yt- c(0.25, 0.10, 0.90, 0.25, 1.05, 1.70, 2.05, 2.90, 2.35, 2.60, 2.55,
2.15, 2.75, 2.05,
2.70, 2.25, 2.55, 2.05, 3.65, 3.05, 3.00, 3.50, 3.75, 4.85, 4.50, 4.50, 3.35,
4.90)
xs- c(0.45, 1.05, 2.75, 3.30, 4.95, 0.40, 1.05, 2.30, 3.45, 4.60, 0.05,
1.95, 2.95, 3.70,
4.55, 0.75, 1.60, 2.10, 3.60, 4.90, 0.05, 1.35, 2.60, 3.40, 4.25)
ys- c(0.45, 0.95, 0.75, 0.95, 0.10, 1.90, 1.45, 1.25, 1.45, 1.05, 2.85,
2.60, 2.05, 2.60,
2.55, 3.75, 3.30, 3.95, 3.45, 3.70, 4.95, 4.35, 4.55, 4.40, 4.95)
carbon- c(1.43, 1.82, 1.40, 1.43, 1.96, 1.61, 1.91, 1.53, 1.17, 1.83, 2.43,
2.02, 1.66,
2.45, 2.46, 1.39, 1.10, 1.38, 1.91, 2.13, 1.88, 1.26, 2.15, 1.89, 1.69)
carbon.df=data.frame(x=xs,y=ys,z=carbon)
carbon.loess= loess(z~x*y, data= carbon.df, degree= 2)
[R] Subscript out of bound error
Hi All I have a huge matrix m (10276 X 10276) dimension with same column name and row names. (its a gene correlation matrix). I have another text file which has 2700 names, basically locus ID of genes, which are also rownames/colnames in m. Now I want to select all those columns from m whose names match with the names in the text file and return ONLY those columns and ALL the rows in a matrix format. That is, the output should a matrix of dimension 10276 X 2700. I have played around with subset and also tried all different methods like reading the text file, converting to matrix and character and DF and tried all possible ways I found on the web. The problem is that some of them work on a smaller simulated dataset but says Subscript out of bounds when I try on the bigger matrix(m). Please help.. :( Thanks! -- *Chirag Gupta* Department of Crop, Soil, and Environmental Sciences, 115 Plant Sciences Building, Fayetteville, Arkansas 72701 [[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] problem on reading many files
Hello, I'm not sure whether it is dur to the HTML version or not, but there is a problem with quotation marks in your first line. Regards, Pascal 2013/7/5 RODRIGUEZ MORENO VICTOR MANUEL rodriguez.vic...@inifap.gob.mx Hi, I have to run almost 120 stations files of temperatura (mx and min), separated, and rain. I am following the instructions on RHtests tutorial as on http://www.cmc.org.ve/mediawiki/index.php?title=Preparando_los_datoslink. Bit I have no success on running on multiple files. I have the .ls file which include the names of the file station, which I include on this msg for your consideration. The instruction that I wrote on the console is listatmn - readLines (tmaxcopia.ls, warn=false) for(ifile in listatmn) FindU(paste(./,ifile,sep=),paste(salidas/,ifile,sep=),-999.9) Of course the tmaxcopia.ls is the list of stations for tmax data. But once running on the console, this is the error msg I got listatmn - readLines (tmaxcopia.ls, warn=false) Error: inesperado entrada in listatmn - readLines ( for(ifile in listatmn) FindU(paste(./,ifile,sep=),paste(salidas/,ifile,sep=),-999.9) Error: objeto 'ifnames' no encontrado What I'd like to do is to make easier and fast running the RHtest Thank's a lot in advance for your help. Dr Victor M Rodríguez M Doctor en Ciencias, en Ciencias de la Tierra / Geociencias Ambientales INIFAP. CEPAB. (465) 9580167 y 86 Ext 108 y 220 / Este msg y su(s) atado(s) tiene carácter de CONFIDENCIAL y solo es para los fines descritos en su contenido. Queda expresamente prohibida bajo la ley de protección de datos del Gobierno de los Estados Unidos Mexicanos y para su difusion o extensión a terceros, se requiere del consentimiento expresso del titular de la cuenta de correo electrónico y/o del administrador del sitio en extenso de sus responsabilidades . / [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R and MatLab implementations of the same model differs
On 05-07-2013, at 09:53, Jannetta Steyn janne...@henning.org wrote: I don't quite know how to explain the doesn't work in more detail without any visual aid. You said that R got into an indefinite loop, whatever that maybe. When I change the solver to ode45 the script never stops running. I have even left it over night at one stage but the next day it was still busy. Is there a way to see what it is doing and to determine why it seems to be in an infinite loop? The script just ran using ode45!! For the first time ever. Please keep the conversation on R-help. Don't reply to me personally. Which script? With the corrections suggested by me and Peter? (gK_axon_AB=52.5-3 == gK_axon_AB=52.5e-3) For what it's worth: lsdoda seems quicker. Variable v_acon_AB now converges to -60 (the equilibrium state?) Berend __ 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] Unique in discerning missing values NA
Hi, I am trying to remove duplicate Patient numbers in a clinical record, I used unique menPatients[1:40,1] [1] abr1160(C)/001 ABR1363(A)/001 ABR1363(A)/001 ABR1363(A)/001 abr1772(B)/001 [6] AFR0003/001AFR0003/001afr0290(C)/001 afr1861(B)/001 Aga0007/001 [11] AGA1548(A)/001 AGA1548(A)/001 AGA1548(A)/001 AGU1680(A)/001 AGU1680(A)/001 [16] AIS0492/001AIS0492/001AKO4268(C)/001 AKO4268(C)/001 AKT0042(B)/001 [21] AKT0042(B)/001 AKT0042(B)/001 AKT0042(B)/001 AKT0042(B)/001 AKT0042(B)/001 [26] AKT0042(B)/001 alb4423(C)/001 ALF1651(A)/001 alf1722(B)/001 ALF1735(A)/001 [31] ALF1735(A)/001 ALP4321(C)/001 NA NA ALU4262(B)/001 [36] ALV4286(C)/001 ALW2579(C)/001 NA ALW4330(B)/001 AMA0011/001 3886 Levels: 0750/002 0751/001 0984/002 ABE2560(C)/001 ... zul1737(B)/001 testData-menPatients[1:40,1] I then used unique, please note the NA at position 32 in testData testUnique-unique(testData) testUnique [1] abr1160(C)/001 ABR1363(A)/001 abr1772(B)/001 AFR0003/001afr0290(C)/001 [6] afr1861(B)/001 Aga0007/001AGA1548(A)/001 AGU1680(A)/001 AIS0492/001 [11] AKO4268(C)/001 AKT0042(B)/001 alb4423(C)/001 ALF1651(A)/001 alf1722(B)/001 [16] ALF1735(A)/001 ALP4321(C)/001 NA ALU4262(B)/001 ALV4286(C)/001 [21] ALW2579(C)/001 ALW4330(B)/001 AMA0011/001 The missing value NA originally at position 32 in testdata is still there, it is in position 18. Why is this? How can I prevent this? I tried using incomprables=c(NA), but this did not work. Thanks Pancho Mulongeni Research Assistant PharmAccess Foundation 1 Fouché Street Windhoek West Windhoek Namibia Tel: +264 61 419 000 Fax: +264 61 419 001/2 Mob: +264 81 4456 286 __ 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] Subscript out of bound error
use 'match' to convert the names to column indices and then use that for indexing indx - match(subCols, names(yourMatrix) mySubset - yourMatrix[, indx] Sent from my iPad On Jul 5, 2013, at 2:22, Chirag Gupta cxg...@email.uark.edu wrote: Hi All I have a huge matrix m (10276 X 10276) dimension with same column name and row names. (its a gene correlation matrix). I have another text file which has 2700 names, basically locus ID of genes, which are also rownames/colnames in m. Now I want to select all those columns from m whose names match with the names in the text file and return ONLY those columns and ALL the rows in a matrix format. That is, the output should a matrix of dimension 10276 X 2700. I have played around with subset and also tried all different methods like reading the text file, converting to matrix and character and DF and tried all possible ways I found on the web. The problem is that some of them work on a smaller simulated dataset but says Subscript out of bounds when I try on the bigger matrix(m). Please help.. :( Thanks! -- *Chirag Gupta* Department of Crop, Soil, and Environmental Sciences, 115 Plant Sciences Building, Fayetteville, Arkansas 72701 [[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.
[R] Revolutions blog: June roundup
Revolution Analytics staff write about R every weekday at the Revolutions blog: http://blog.revolutionanalytics.com and every month I post a summary of articles from the previous month of particular interest to readers of r-help. In case you missed them, here are some articles related to R from the month of June: You can create a Word document from a template and an R script with the R2DOCX package: http://bit.ly/17XigAJ Joe Rickert reviews books and other resources for learning about time series analysis in R: http://bit.ly/17Xiisn Timely Portfolio covers 15 years of history of time series plotting with R: http://bit.ly/17Xiisp An online beer recommendation application serves as an in-depth example of building a recommendation system with R: http://bit.ly/17XigAK Software company SAP says that skills around the open source R programming language and advanced analytics are rapidly shifting from 'niche' to 'standard' requirements: http://bit.ly/17Xiiso A primer on maximum likelihood estimation in R with the bbmle package: http://bit.ly/17XigAI American Century Investments describe how they created their own R package to optimize investments and model supplier relationships in a video presentation: http://bit.ly/17Xiisq How to draw attractive decision trees using the rpart.plot package: http://bit.ly/17Xiisr What is a data scientist, what skills should they have, and how is the practice evolving? My take: http://bit.ly/17XigAL Computerworld's 6-part beginner's guide to R: http://bit.ly/17XigAM RStudio has made CRAN download statistics available, enabling a ranking of the top R packages by downloads: http://bit.ly/17XigAN Big Data and statistical modeling is changing video games: helping to identify bottlenecks for new players, detect fraud, and maximize in-app purchases. Details in this webinar replay: http://bit.ly/17Xiiss A mini-tutorial on using the Quandl package to import public financial data sets into R: http://bit.ly/17Xiist Dirk Eddelbuettel's book, Seamless R and C++ Integration with Rcpp, is now available: http://bit.ly/17XigAO An interactive application based on R and Shiny maps dialect differences across the USA: http://bit.ly/17XigAP Generating parallel random number streams with the RevoScaleR package: http://bit.ly/17Xiisu R again shows strong growth in the annual KDNuggets software poll: http://bit.ly/17XigAQ Some very useful guidelines for setting up a reproducible R project: http://bit.ly/17Xiisv Some non-R stories in the past month included: A new Hadoop appliance from Teradata (http://bit.ly/17XigAT), a new Big Data Innovation Center in Singapore (http://bit.ly/17XigAU), a drum/keyboard Billie Jean cover (http://bit.ly/17XigAW), women in movies (http://bit.ly/17XigAV), a Daft Punk / Soul Train mashup (http://bit.ly/17Xiisw), and It Gets Better at NASA (http://bit.ly/17Xiisx). Meeting times for local R user groups (http://blog.revolutionanalytics.com/local-r-groups.html) can be found on the updated R Community Calendar at: http://blog.revolutionanalytics.com/calendar.html If you're looking for more articles about R, you can find summaries from previous months at http://blog.revolutionanalytics.com/roundups/. Join the Revolution mailing list at http://revolutionanalytics.com/newsletter to be alerted to new articles on a monthly basis. As always, thanks for the comments and please keep sending suggestions to me at da...@revolutionanalytics.com . Don't forget you can also follow the blog using an RSS reader, or by following me on Twitter (I'm @revodavid). Cheers, # David -- David M Smith da...@revolutionanalytics.com VP of Marketing, Revolution Analytics http://blog.revolutionanalytics.com Tel: +1 (650) 646-9523 (Seattle WA, USA) Twitter: @revodavid We're hiring! www.revolutionanalytics.com/careers __ 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] Meta-analysis on a repeated measures design with multiple trials per subject using metafor
Dear Wolfgang and other readers of the r-help list, Thank you very much for your suggestion. Unfortunately, the data that I have can not be described with a table such as the one you have made, because there's no identical trial under both treatment 1 and treatment 2. To explain, let me explain a bit more about the experiments: * All subjects were presented with the same number of trials * Half of these trials were preceded by a prime from category 1 (treatment 1) and half of these trials with a prime from category 2 (treatment 2) * Subjects were asked to respond to these trials (a unique stimulus for each trial) by pressing one of two keys on the keyboard. Because everything was randomized, I can only calculate the total number of times a certain response was used under each type of trial. There is no pairing of trials under two treatments, so I am forced to use the marginal totals from your table. I have uploaded a simplified version of the data for one experiment to illustrate this (the actual experiments have five treatments and some have moderators): https://www.dropbox.com/s/rhgo12cm1asl6x8/exampledata.csv This is the script that I used to generate the data: https://www.dropbox.com/s/7uyeaexhnqiiy55/exampledata.R The problem thus appears to lie mainly in estimating the variance of the proportion difference from only the marginal totals, is that correct? Is there a way to calculate it from only the marginal totals? One alternative that I have tried over the last few days, is to use the b parameter of interest and it's corresponding standard error from the lme4 regression output that I use to analyse the individual experiments. Then, I use rma(yi, sei) to do a meta-analysis on these parameters. I am not sure this is correct though, since it takes into account between-subjects variance (through a random effect for subject), and it is sensitive to the covariates/moderators I include in the models that I get the b parameters from. Thanks again for your help, and for any suggestions for solving this problem! Regards, Marc On 07/04/2013 11:21 PM, Viechtbauer Wolfgang (STAT) wrote: Dear Marc, Let me see if I understand the type of data you have. You say that you have 5 experiments. And within each experiment, you have n subjects and for each subject, you have data in the form described in your post. Now for each subject, you want to calculate some kind of measure that quantifies how much more likely it was that subjects gave/chose response 2 under treatment 2 versus treatment 1. So, you would have n such values. And then you want to pool those values over the n subjects within a particular experiment and then ultimately over the 5 experiments. Is that correct so far? Assuming I got this right, let me ask you about those data that you have for each subject. In particular, are these paired data? In other words, is there are 1:1 relationship between the 30 trials under treatment 1 versus treatment 2? Or phrased yet another way, can you construct a table like this for every subject: trt 2 resp1 resp2 trt 1 resp1 a b 10 resp2 c d 20 2010 30 Note that I added the marginal counts based on your example data, but this is not sufficient to reconstruct how often response 1 was chosen for the same trial under both treatment 1 and treatment 2 (cell a). And so on for the other 3 cells. If all of this applies, then essentially you are dealing with dependent proportions and you can calculate the difference y = (20/30)-(10/30) as you have done. The corresponding sampling variance can be estimated with v = var(y) = (a+b)*(c+d)/t^3 + (a+c)*(b+d)/t^3 - 2*(a*d/t^3 - b*c/t^3) (where t is the number of trials, i.e., 30 in the example above). See, for example, section 10.1.1. in Agresti (2002) (Categorical data analysis, 2nd ed.). So, ultimately, you will have n values of y and v for a particular experiment and then the same thing for all 5 experiments. You can then pool those values with rma(yi, vi) in metafor (yi and vi being the vectors of the y and v values). You probably want to add a factor to the model that indicates which experiment those values came from. So, something like: rma(yi, vi, mods = ~ factor(experiment)). Well, I hope that I understood your data correctly. Best, Wolfgang -- Wolfgang Viechtbauer, Ph.D., Statistician Department of Psychiatry and Psychology School for Mental Health and Neuroscience Faculty of Health, Medicine, and Life Sciences Maastricht University, P.O. Box 616 (VIJV1) 6200 MD Maastricht, The Netherlands +31 (43) 388-4170 | http://www.wvbauer.com From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] On Behalf Of Marc Heerdink [m.w.heerd...@uva.nl] Sent: Wednesday, July 03, 2013 2:15 PM To: r-help@r-project.org Subject: [R] Meta-analysis on a repeated measures design
[R] To approach time series using a data set
Dear R community, I have to approach a time serie as a linear combination of many time series in a dataset. Of course, I would like that this linear combination is optimal (I would like to use the minimal number of times series) although I can loose some information. Moreover, each variable of the dataset has associated a cost which range from 0 to 9. I also would like to find this optimal combination and to make minimal de sum of this cost. Could you give to me some suggestions? Do your know some interesting reference for starting? Thanks in advance __ 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] Unique in discerning missing values NA
Hi, testUnique - unique(testData[!is.na(testData)]) or testUnique - unique(na.omit(testData)) And probably some other solutions. Regards, Pascal 2013/7/5 Pancho Mulongeni p.mulong...@namibia.pharmaccess.org Hi, I am trying to remove duplicate Patient numbers in a clinical record, I used unique menPatients[1:40,1] [1] abr1160(C)/001 ABR1363(A)/001 ABR1363(A)/001 ABR1363(A)/001 abr1772(B)/001 [6] AFR0003/001AFR0003/001afr0290(C)/001 afr1861(B)/001 Aga0007/001 [11] AGA1548(A)/001 AGA1548(A)/001 AGA1548(A)/001 AGU1680(A)/001 AGU1680(A)/001 [16] AIS0492/001AIS0492/001AKO4268(C)/001 AKO4268(C)/001 AKT0042(B)/001 [21] AKT0042(B)/001 AKT0042(B)/001 AKT0042(B)/001 AKT0042(B)/001 AKT0042(B)/001 [26] AKT0042(B)/001 alb4423(C)/001 ALF1651(A)/001 alf1722(B)/001 ALF1735(A)/001 [31] ALF1735(A)/001 ALP4321(C)/001 NA NA ALU4262(B)/001 [36] ALV4286(C)/001 ALW2579(C)/001 NA ALW4330(B)/001 AMA0011/001 3886 Levels: 0750/002 0751/001 0984/002 ABE2560(C)/001 ... zul1737(B)/001 testData-menPatients[1:40,1] I then used unique, please note the NA at position 32 in testData testUnique-unique(testData) testUnique [1] abr1160(C)/001 ABR1363(A)/001 abr1772(B)/001 AFR0003/001 afr0290(C)/001 [6] afr1861(B)/001 Aga0007/001AGA1548(A)/001 AGU1680(A)/001 AIS0492/001 [11] AKO4268(C)/001 AKT0042(B)/001 alb4423(C)/001 ALF1651(A)/001 alf1722(B)/001 [16] ALF1735(A)/001 ALP4321(C)/001 NA ALU4262(B)/001 ALV4286(C)/001 [21] ALW2579(C)/001 ALW4330(B)/001 AMA0011/001 The missing value NA originally at position 32 in testdata is still there, it is in position 18. Why is this? How can I prevent this? I tried using incomprables=c(NA), but this did not work. Thanks Pancho Mulongeni Research Assistant PharmAccess Foundation 1 Fouché Street Windhoek West Windhoek Namibia Tel: +264 61 419 000 Fax: +264 61 419 001/2 Mob: +264 81 4456 286 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] kruskal.test followed by kruskalmc
Hi all, After running kruskal.test I have got results (p0,005) pointing to reject the hypothesis that the samples were draw from the same population. Howerver when I run the kruskalmc there are no significant differences in any of the multiple comparisons. Is that possible? Some clarification? Thanks, Humber https://sites.google.com/site/humberandrade [[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] Unique in discerning missing values NA
Hello, Your data example is difficult to read into an R session. Next time, post the output of ?dput. Like this: dput(menPatients[1:40, 1]) # post the output of this The help page for unique says that Missing values are regarded as equal so you should expect one NA to still be present in the final result. If you want to remove NAs, use ?is.na. With fake data, x1 - c(1:3, NA, 4, NA, 2:9) x2 - unique(x1) x3 - x2[!is.na(x2)] x3 Hope this helps, Rui Barradas Em 05-07-2013 10:28, Pancho Mulongeni escreveu: Hi, I am trying to remove duplicate Patient numbers in a clinical record, I used unique menPatients[1:40,1] [1] abr1160(C)/001 ABR1363(A)/001 ABR1363(A)/001 ABR1363(A)/001 abr1772(B)/001 [6] AFR0003/001AFR0003/001afr0290(C)/001 afr1861(B)/001 Aga0007/001 [11] AGA1548(A)/001 AGA1548(A)/001 AGA1548(A)/001 AGU1680(A)/001 AGU1680(A)/001 [16] AIS0492/001AIS0492/001AKO4268(C)/001 AKO4268(C)/001 AKT0042(B)/001 [21] AKT0042(B)/001 AKT0042(B)/001 AKT0042(B)/001 AKT0042(B)/001 AKT0042(B)/001 [26] AKT0042(B)/001 alb4423(C)/001 ALF1651(A)/001 alf1722(B)/001 ALF1735(A)/001 [31] ALF1735(A)/001 ALP4321(C)/001 NA NA ALU4262(B)/001 [36] ALV4286(C)/001 ALW2579(C)/001 NA ALW4330(B)/001 AMA0011/001 3886 Levels: 0750/002 0751/001 0984/002 ABE2560(C)/001 ... zul1737(B)/001 testData-menPatients[1:40,1] I then used unique, please note the NA at position 32 in testData testUnique-unique(testData) testUnique [1] abr1160(C)/001 ABR1363(A)/001 abr1772(B)/001 AFR0003/001afr0290(C)/001 [6] afr1861(B)/001 Aga0007/001AGA1548(A)/001 AGU1680(A)/001 AIS0492/001 [11] AKO4268(C)/001 AKT0042(B)/001 alb4423(C)/001 ALF1651(A)/001 alf1722(B)/001 [16] ALF1735(A)/001 ALP4321(C)/001 NA ALU4262(B)/001 ALV4286(C)/001 [21] ALW2579(C)/001 ALW4330(B)/001 AMA0011/001 The missing value NA originally at position 32 in testdata is still there, it is in position 18. Why is this? How can I prevent this? I tried using incomprables=c(NA), but this did not work. Thanks Pancho Mulongeni Research Assistant PharmAccess Foundation 1 Fouché Street Windhoek West Windhoek Namibia Tel: +264 61 419 000 Fax: +264 61 419 001/2 Mob: +264 81 4456 286 __ 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] Unique in discerning missing values NA
Yes thanks, this is what I ended up doing, but I though there would be a 'internal' way to disregard NAs in unique. Thanks for the tip on dput -Original Message- From: Rui Barradas [mailto:ruipbarra...@sapo.pt] Sent: 05 July 2013 11:39 To: Pancho Mulongeni Cc: r-help@r-project.org Subject: Re: [R] Unique in discerning missing values NA Hello, Your data example is difficult to read into an R session. Next time, post the output of ?dput. Like this: dput(menPatients[1:40, 1]) # post the output of this The help page for unique says that Missing values are regarded as equal so you should expect one NA to still be present in the final result. If you want to remove NAs, use ?is.na. With fake data, x1 - c(1:3, NA, 4, NA, 2:9) x2 - unique(x1) x3 - x2[!is.na(x2)] x3 Hope this helps, Rui Barradas Em 05-07-2013 10:28, Pancho Mulongeni escreveu: Hi, I am trying to remove duplicate Patient numbers in a clinical record, I used unique menPatients[1:40,1] [1] abr1160(C)/001 ABR1363(A)/001 ABR1363(A)/001 ABR1363(A)/001 abr1772(B)/001 [6] AFR0003/001AFR0003/001afr0290(C)/001 afr1861(B)/001 Aga0007/001 [11] AGA1548(A)/001 AGA1548(A)/001 AGA1548(A)/001 AGU1680(A)/001 AGU1680(A)/001 [16] AIS0492/001AIS0492/001AKO4268(C)/001 AKO4268(C)/001 AKT0042(B)/001 [21] AKT0042(B)/001 AKT0042(B)/001 AKT0042(B)/001 AKT0042(B)/001 AKT0042(B)/001 [26] AKT0042(B)/001 alb4423(C)/001 ALF1651(A)/001 alf1722(B)/001 ALF1735(A)/001 [31] ALF1735(A)/001 ALP4321(C)/001 NA NA ALU4262(B)/001 [36] ALV4286(C)/001 ALW2579(C)/001 NA ALW4330(B)/001 AMA0011/001 3886 Levels: 0750/002 0751/001 0984/002 ABE2560(C)/001 ... zul1737(B)/001 testData-menPatients[1:40,1] I then used unique, please note the NA at position 32 in testData testUnique-unique(testData) testUnique [1] abr1160(C)/001 ABR1363(A)/001 abr1772(B)/001 AFR0003/001 afr0290(C)/001 [6] afr1861(B)/001 Aga0007/001AGA1548(A)/001 AGU1680(A)/001 AIS0492/001 [11] AKO4268(C)/001 AKT0042(B)/001 alb4423(C)/001 ALF1651(A)/001 alf1722(B)/001 [16] ALF1735(A)/001 ALP4321(C)/001 NA ALU4262(B)/001 ALV4286(C)/001 [21] ALW2579(C)/001 ALW4330(B)/001 AMA0011/001 The missing value NA originally at position 32 in testdata is still there, it is in position 18. Why is this? How can I prevent this? I tried using incomprables=c(NA), but this did not work. Thanks Pancho Mulongeni Research Assistant PharmAccess Foundation 1 Fouché Street Windhoek West Windhoek Namibia Tel: +264 61 419 000 Fax: +264 61 419 001/2 Mob: +264 81 4456 286 __ 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] kruskal.test followed by kruskalmc
Humber, Have a look at this: http://r.789695.n4.nabble.com/Multiple-Comparisons-Kruskal-Wallis-Test-kruskal-agricolae-and-kruskalmc-pgirmess-don-t-yield-the-sa-td4639004.html Hope it helps. Kind regards, José Prof. José Iparraguirre Chief Economist Age UK -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Humber Andrade Sent: 05 July 2013 11:38 To: r-help@r-project.org Subject: [R] kruskal.test followed by kruskalmc Hi all, After running kruskal.test I have got results (p0,005) pointing to reject the hypothesis that the samples were draw from the same population. Howerver when I run the kruskalmc there are no significant differences in any of the multiple comparisons. Is that possible? Some clarification? Thanks, Humber https://sites.google.com/site/humberandrade [[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. The Wireless from Age UK | Radio for grown-ups. www.ageuk.org.uk/thewireless If you’re looking for a radio station that offers real variety, tune in to The Wireless from Age UK. Whether you choose to listen through the website at www.ageuk.org.uk/thewireless, on digital radio (currently available in London and Yorkshire) or through our TuneIn Radio app, you can look forward to an inspiring mix of music, conversation and useful information 24 hours a day. --- Age UK is a registered charity and company limited by guarantee, (registered charity number 1128267, registered company number 6825798). Registered office: Tavis House, 1-6 Tavistock Square, London WC1H 9NA. For the purposes of promoting Age UK Insurance, Age UK is an Appointed Representative of Age UK Enterprises Limited, Age UK is an Introducer Appointed Representative of JLT Benefit Solutions Limited and Simplyhealth Access for the purposes of introducing potential annuity and health cash plans customers respectively. Age UK Enterprises Limited, JLT Benefit Solutions Limited and Simplyhealth Access are all authorised and regulated by the Financial Services Authority. -- This email and any files transmitted with it are confidential and intended solely for the use of the individual or entity to whom they are addressed. If you receive a message in error, please advise the sender and delete immediately. Except where this email is sent in the usual course of our business, any opinions expressed in this email are those of the author and do not necessarily reflect the opinions of Age UK or its subsidiaries and associated companies. Age UK monitors all e-mail transmissions passing through its network and may block or modify mails which are deemed to be unsuitable. Age Concern England (charity number 261794) and Help the Aged (charity number 272786) and their trading and other associated companies merged on 1st April 2009. Together they have formed the Age UK Group, dedicated to improving the lives of people in later life. The three national Age Concerns in Scotland, Northern Ireland and Wales have also merged with Help the Aged in these nations to form three registered charities: Age Scotland, Age NI, Age Cymru. __ 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 and MatLab implementations of the same model differs
On 5 July 2013 09:44, Berend Hasselman b...@xs4all.nl wrote: On 05-07-2013, at 09:53, Jannetta Steyn janne...@henning.org wrote: I don't quite know how to explain the doesn't work in more detail without any visual aid. You said that R got into an indefinite loop, whatever that maybe. When I change the solver to ode45 the script never stops running. I have even left it over night at one stage but the next day it was still busy. Is there a way to see what it is doing and to determine why it seems to be in an infinite loop? The script just ran using ode45!! For the first time ever. Please keep the conversation on R-help. Don't reply to me personally. Apologies for that. It was meant to go to the list. Which script? With the corrections suggested by me and Peter? (gK_axon_AB=52.5-3 == gK_axon_AB=52.5e-3) The R script. With the correction of AB=52.5e-3 For what it's worth: lsdoda seems quicker. Yes, lsoda is quicker but I want to use ode45 because that is what they used in the paper I am referencing and I am trying to get the same results as they did. Variable v_acon_AB now converges to -60 (the equilibrium state?) I'm trying to determine now whether this is correct. Regards Jannetta -- === Web site: http://www.jannetta.com Email: janne...@henning.org === [[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] iterative methods
Dear R users,
Please help me with some documentation for newbie about R
programming, algorithms, create iterative C++ function (like
for, while, if , etc).
Start at
http://www.R-project.org/posting-guide.html
and especiallly the resources cited towards the end of that page, just after
the bit titled Do your homework before posting. It lists useful links.
Introduction to R and the R language definition (sec 3.2 especially) seem
relevant.
For books, browse Google for 'R programming books'; there are lots. Some are
aimed at particular types of problem, but you know your problem best.
S Ellison
***
This email and any attachments are confidential. Any use...{{dropped:8}}
__
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 when building a custom package
I only get an error message, when trying to build a binary
file (or runnning /R CMD check package name/), not the
standard tar.gz. Here is what the output looks like in the
command prompt:
...
I have gotten the impression that generally the error with
the Error in file... is a problem with a wrong path or
problems with permissions,
Not sure whether it will solve your problem, but specifying a writeable
location for R CMD INSTALL --build (as described in 1.3.3 Building binary
packages in Writing R extensions) should dodge a permission problem if that
is the issue. Personally I prefer to do that anyway to avoid overwriting an
existing (working) package installation.
As an aside - probably unrelated and apparently not a warning anyway - I don't
see why you would need to be be creating a new plot generic. Isn't the existing
plot generic sufficient?
S Ellison
***
This email and any attachments are confidential. Any use...{{dropped:8}}
__
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] kruskal.test followed by kruskalmc
Thank you Prof. José Iparraguirre. Maybe I am wrong but I think the issues are not the same. His data doesn't showed significant differences after kruskal.test(), that was not my case. Anyway follow below my the results I've got and the database. Thank you, # kruskal.test(data$resp,data$group) Kruskal-Wallis rank sum test data: data$resp and data$group Kruskal-Wallis chi-squared = 32.3546, df = 14, p-value = 0.003566 kruskalmc(data$resp,data$group) Multiple comparison test after Kruskal-Wallis p.value: 0.05 Comparisons obs.dif critical.dif difference A-B 4.8303571 62.37688 FALSE A-C 3.8928571 62.37688 FALSE A-D 0.4821429 62.37688 FALSE . . . M-P 14.250 60.26179 FALSE N-O 1.375 60.26179 FALSE N-P 6.125 60.26179 FALSE O-P 4.750 60.26179 FALSE # database group resp A 0.1 A 581.8 A 90.5 A 70.1 A 820.1 A 1159.2 A 2478.1 A 2475.3 B 351.8 B 370.1 B 326.1 B 751.9 B 931 B 588.2 B 70.1 B 1754.9 C 289.8 C 254.1 C 370.3 C 459.8 C 412.5 C 591.5 C 986.9 C 890 D 425.6 D 397.4 D 464 D 370.9 D 417.3 D 455 D 568.2 D 599.4 E 405.1 E 626.2 E 299 E 493.8 E 362.6 E 309.8 E 522.7 E 433.3 F 698.6 F 42.5 F 7.4 F 10.6 F 95.8 F 497.5 F 987.9 F 925.1 G 492.9 G 376 G 413 G 278.3 G 344.2 G 292.2 G 429.4 G 368 H 241.6 H 230.5 H 310.4 H 372.5 H 366.1 H 307.9 H 480 H 529.8 I 296 I 288.8 I 302.1 I 300.8 I 150.1 I 381.9 I 583.1 I 489.4 J 1.2 J 18.6 J 7.7 J 11.6 J 48.1 J 121.8 J 1284.1 J 944.7 L 0.5 L 44.4 L 80.9 L 15.3 L 80 L 379.9 L 940.6 L 829.3 M 323.6 M 401.5 M 162.1 M 136.5 M 139.4 M 363.3 M 280.7 M 356.5 N 197.6 N 245.9 N 221.5 N 224.3 N 185.4 N 265.3 N 304.8 N 351.9 O 189.9 O 237.3 O 247.1 O 230.4 O 272.2 O 155.1 O 270.7 O 315.2 P 48.4 P 15.5 P 53.1 P 72.8 P 74.8 P 132.3 P 550 P 478.7 On Fri, Jul 5, 2013 at 8:06 AM, Jose Iparraguirre jose.iparragui...@ageuk.org.uk wrote: Humber, Have a look at this: http://r.789695.n4.nabble.com/Multiple-Comparisons-Kruskal-Wallis-Test-kruskal-agricolae-and-kruskalmc-pgirmess-don-t-yield-the-sa-td4639004.html Hope it helps. Kind regards, José Prof. José Iparraguirre Chief Economist Age UK -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Humber Andrade Sent: 05 July 2013 11:38 To: r-help@r-project.org Subject: [R] kruskal.test followed by kruskalmc Hi all, After running kruskal.test I have got results (p0,005) pointing to reject the hypothesis that the samples were draw from the same population. Howerver when I run the kruskalmc there are no significant differences in any of the multiple comparisons. Is that possible? Some clarification? Thanks, Humber https://sites.google.com/site/humberandrade [[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. The Wireless from Age UK | Radio for grown-ups. www.ageuk.org.uk/thewireless If youre looking for a radio station that offers real variety, tune in to The Wireless from Age UK. Whether you choose to listen through the website at www.ageuk.org.uk/thewireless, on digital radio (currently available in London and Yorkshire) or through our TuneIn Radio app, you can look forward to an inspiring mix of music, conversation and useful information 24 hours a day. --- Age UK is a registered charity and company limited by guarantee, (registered charity number 1128267, registered company number 6825798). Registered office: Tavis House, 1-6 Tavistock Square, London WC1H 9NA. For the purposes of promoting Age UK Insurance, Age UK is an Appointed Representative of Age UK Enterprises Limited, Age UK is an Introducer Appointed Representative of JLT Benefit Solutions Limited and Simplyhealth Access for the purposes of introducing potential annuity and health cash plans customers respectively. Age UK Enterprises Limited, JLT Benefit Solutions Limited and Simplyhealth Access are all authorised and regulated by the Financial Services Authority. -- This email and any files transmitted with it are confidential and intended solely for the use of the individual or entity to whom they are addressed. If you receive a message in error, please advise the sender and delete immediately. Except where this email is sent in the usual course of our business, any opinions expressed in this email
Re: [R] ggplot2
Thank you Mr. Kane for your time, I finally achieve my objective with the package reshape2 n5dt-last(dezdiff,5) n5dt - as.data.frame(t(n5dt)) n5dt$id - c(1:35) subm5dt-melt(n5dt, id=id, c(2013-06-28, 2013-07-01, 2013-07-02, 2013-07-03, 2013-07-04) ) names(subm5dt) - c(Observations, variable,BasisPoints) ggplot(subm5dt, aes(Observations, BasisPoints, col=variable, group=2)) + geom_line() + facet_grid(variable~., scale=fixed) + geom_point() Have a nice day! -- View this message in context: http://r.789695.n4.nabble.com/ggplot2-tp4670852p4670931.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] repolr: multivariate repeated analysis possible?
Dear all, I'm trying to perform a non-parametric multivariate repeated measures analysis of 9 ordered variables (scale 0-3) at two time points. So basically a multivariate repeated measure GLM for ordered variables. While the package repolr can model one ordered variable over time, I have not found a way to do this with several DVs at the same time (I am interested in the time * DVs interaction term, so I need a multivariate model). Is this possible? Is there example syntax for this problem? Thank you for your help Torvon [[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] Error when building a custom package
Hi Ellison
Thanks for the reply. The test package build I did using the skeleton package
was in the same folder as the one I'm using now for the custom 'slo' package,
so that is not the problem.
I just investigated whether it was an issue with program rights regarding
downloading and installing new packages from RCMD, but even after double
checking that all packages were installed (and finding some that wasn't,. which
are now manually installed), it still doesn't work.
Regarding the new plot implementation, I don't have a clue whether that is
needed, but I fully trust the capabilities of the guy to made the source files
- and as written previously - he is able to build the package on his machine.
Venlig hilsen / Best regards
Christer P. Volk
Specialist
SenseLab
DELTA| Venlighedsvej 4 | 2970 Hørsholm | Denmark |
madebydelta.comhttp://madebydelta.com
c...@delta.dkmailto:c...@delta.dk | tel: +45 72 19 40 00 | direct: +45 72 19
43 45
[cid:image001.jpg@01CE7995.A881CD10]
We help ideas meet the real world
[cid:image002.jpg@01CE7995.A881CD10]http://www.madebydelta.com/delta/campaign/newsletter-signup.page?
From: S Ellison-2 [via R] [mailto:ml-node+s789695n4670933...@n4.nabble.com]
Sent: 5. juli 2013 15:07
To: Christer P. Volk
Subject: Re: Error when building a custom package
I only get an error message, when trying to build a binary
file (or runnning /R CMD check package name/), not the
standard tar.gz. Here is what the output looks like in the
command prompt:
...
I have gotten the impression that generally the error with
the Error in file... is a problem with a wrong path or
problems with permissions,
Not sure whether it will solve your problem, but specifying a writeable
location for R CMD INSTALL --build (as described in 1.3.3 Building binary
packages in Writing R extensions) should dodge a permission problem if that
is the issue. Personally I prefer to do that anyway to avoid overwriting an
existing (working) package installation.
As an aside - probably unrelated and apparently not a warning anyway - I don't
see why you would need to be be creating a new plot generic. Isn't the existing
plot generic sufficient?
S Ellison
***
This email and any attachments are confidential. Any use...{{dropped:8}}
__
[hidden email]/user/SendEmail.jtp?type=nodenode=4670933i=0 mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
If you reply to this email, your message will be added to the discussion below:
http://r.789695.n4.nabble.com/Error-when-building-a-custom-package-tp4670866p4670933.html
To unsubscribe from Error when building a custom package, click
herehttp://r.789695.n4.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_codenode=4670866code=Y3ZvQGRlbHRhLmRrfDQ2NzA4NjZ8LTE5MDc2NTQ5MTc=.
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
The information transmitted is intended only for the person or entity to which
it is addressed and may contain confidential and/or privileged material. Any
review, retransmission, dissemination or other use of, or taking of any action
in reliance upon, this information by persons or entities other than the
intended recipient is prohibited. If you received this in error, please contact
the sender and delete the material from any computer.
image001.jpg (37K)
http://r.789695.n4.nabble.com/attachment/4670938/0/image001.jpg
image002.jpg (38K)
http://r.789695.n4.nabble.com/attachment/4670938/1/image002.jpg
--
View this message in context:
http://r.789695.n4.nabble.com/Error-when-building-a-custom-package-tp4670866p4670938.html
Sent from the R help mailing list archive at Nabble.com.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] kruskal.test followed by kruskalmc
On Jul 5, 2013, at 15:00 , Humber Andrade wrote: Thank you Prof. José Iparraguirre. Maybe I am wrong but I think the issues are not the same. His data doesn't showed significant differences after kruskal.test(), that was not my case. Anyway follow below my the results I've got and the database. This can happen. It is a matter of probability theory, not of R. The following is a simplified paraphrase of what is going on: # 15 random normals, compare range test to variance test # Simulate everything for simplicity # Null distribution M0 - replicate(1, rnorm(15)) vv - apply(M0,2,var) rg - apply(M0,2,range) rr - apply(rg,2,diff) r.95 - quantile(rr, .95) v.95 - quantile(vv, .95) v.995 - quantile(vv, .995) # Distribution at quadrupled variance M - replicate(1, rnorm(15,sd=2)) vv - apply(M,2,var) rg - apply(M,2,range) rr - apply(rg,2,diff) plot(rr,vv) abline(h=c(v.95,v.995),v=r.95, col=red) Notice that the two statistics are correlated, but not equivalent. There are cases where one value is beyond the .95 level and the other is not. Since it is theoretically optimal to use the variance as the test statistic in this model, there are quite a few more cases where rr is below the cutoff and vv is above than the other way around. There are even a sizable number of cases where vv is beyond v.995 and rr does not reach r.95. (The theoretical optimality applies only because I use an increased variance alternative. For specific alternatives, the picture can change. Try it, for instance with a single mean substantially different from the others: M - replicate(1, rnorm(15,mean=rep(c(4,0),c(1,14 ) Thank you, # kruskal.test(data$resp,data$group) Kruskal-Wallis rank sum test data: data$resp and data$group Kruskal-Wallis chi-squared = 32.3546, df = 14, p-value = 0.003566 kruskalmc(data$resp,data$group) Multiple comparison test after Kruskal-Wallis p.value: 0.05 Comparisons obs.dif critical.dif difference A-B 4.8303571 62.37688 FALSE A-C 3.8928571 62.37688 FALSE A-D 0.4821429 62.37688 FALSE . . . M-P 14.250 60.26179 FALSE N-O 1.375 60.26179 FALSE N-P 6.125 60.26179 FALSE O-P 4.750 60.26179 FALSE # database group resp A 0.1 A 581.8 A 90.5 A 70.1 A 820.1 A 1159.2 A 2478.1 A 2475.3 B 351.8 B 370.1 B 326.1 B 751.9 B 931 B 588.2 B 70.1 B 1754.9 C 289.8 C 254.1 C 370.3 C 459.8 C 412.5 C 591.5 C 986.9 C 890 D 425.6 D 397.4 D 464 D 370.9 D 417.3 D 455 D 568.2 D 599.4 E 405.1 E 626.2 E 299 E 493.8 E 362.6 E 309.8 E 522.7 E 433.3 F 698.6 F 42.5 F 7.4 F 10.6 F 95.8 F 497.5 F 987.9 F 925.1 G 492.9 G 376 G 413 G 278.3 G 344.2 G 292.2 G 429.4 G 368 H 241.6 H 230.5 H 310.4 H 372.5 H 366.1 H 307.9 H 480 H 529.8 I 296 I 288.8 I 302.1 I 300.8 I 150.1 I 381.9 I 583.1 I 489.4 J 1.2 J 18.6 J 7.7 J 11.6 J 48.1 J 121.8 J 1284.1 J 944.7 L 0.5 L 44.4 L 80.9 L 15.3 L 80 L 379.9 L 940.6 L 829.3 M 323.6 M 401.5 M 162.1 M 136.5 M 139.4 M 363.3 M 280.7 M 356.5 N 197.6 N 245.9 N 221.5 N 224.3 N 185.4 N 265.3 N 304.8 N 351.9 O 189.9 O 237.3 O 247.1 O 230.4 O 272.2 O 155.1 O 270.7 O 315.2 P 48.4 P 15.5 P 53.1 P 72.8 P 74.8 P 132.3 P 550 P 478.7 On Fri, Jul 5, 2013 at 8:06 AM, Jose Iparraguirre jose.iparragui...@ageuk.org.uk wrote: Humber, Have a look at this: http://r.789695.n4.nabble.com/Multiple-Comparisons-Kruskal-Wallis-Test-kruskal-agricolae-and-kruskalmc-pgirmess-don-t-yield-the-sa-td4639004.html Hope it helps. Kind regards, José Prof. José Iparraguirre Chief Economist Age UK -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Humber Andrade Sent: 05 July 2013 11:38 To: r-help@r-project.org Subject: [R] kruskal.test followed by kruskalmc Hi all, After running kruskal.test I have got results (p0,005) pointing to reject the hypothesis that the samples were draw from the same population. Howerver when I run the kruskalmc there are no significant differences in any of the multiple comparisons. Is that possible? Some clarification? Thanks, Humber https://sites.google.com/site/humberandrade [[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. The Wireless from Age UK | Radio for grown-ups. www.ageuk.org.uk/thewireless If you‚re looking for a radio station
Re: [R] kruskal.test followed by kruskalmc
Thank you Dr. Dalgaard. I understood. I know that this list is not to discuss statistics but I would be very glad if you or someone else can give me some opinion on how to proceed. The kruskal.test says there are differences but the multiple comparisons do not point out what are the differences. Can you suggest a suitable way (maybe paired wilcoxon) to infer what are the differences? I am asking for the hint because I am sure the journal editor/reviewer will ask me to point out which groups differ from each other. Regards, Humber On Fri, Jul 5, 2013 at 11:04 AM, peter dalgaard pda...@gmail.com wrote: On Jul 5, 2013, at 15:00 , Humber Andrade wrote: Thank you Prof. José Iparraguirre. Maybe I am wrong but I think the issues are not the same. His data doesn't showed significant differences after kruskal.test(), that was not my case. Anyway follow below my the results I've got and the database. This can happen. It is a matter of probability theory, not of R. The following is a simplified paraphrase of what is going on: # 15 random normals, compare range test to variance test # Simulate everything for simplicity # Null distribution M0 - replicate(1, rnorm(15)) vv - apply(M0,2,var) rg - apply(M0,2,range) rr - apply(rg,2,diff) r.95 - quantile(rr, .95) v.95 - quantile(vv, .95) v.995 - quantile(vv, .995) # Distribution at quadrupled variance M - replicate(1, rnorm(15,sd=2)) vv - apply(M,2,var) rg - apply(M,2,range) rr - apply(rg,2,diff) plot(rr,vv) abline(h=c(v.95,v.995),v=r.95, col=red) Notice that the two statistics are correlated, but not equivalent. There are cases where one value is beyond the .95 level and the other is not. Since it is theoretically optimal to use the variance as the test statistic in this model, there are quite a few more cases where rr is below the cutoff and vv is above than the other way around. There are even a sizable number of cases where vv is beyond v.995 and rr does not reach r.95. (The theoretical optimality applies only because I use an increased variance alternative. For specific alternatives, the picture can change. Try it, for instance with a single mean substantially different from the others: M - replicate(1, rnorm(15,mean=rep(c(4,0),c(1,14 ) Thank you, # kruskal.test(data$resp,data$group) Kruskal-Wallis rank sum test data: data$resp and data$group Kruskal-Wallis chi-squared = 32.3546, df = 14, p-value = 0.003566 kruskalmc(data$resp,data$group) Multiple comparison test after Kruskal-Wallis p.value: 0.05 Comparisons obs.dif critical.dif difference A-B 4.8303571 62.37688 FALSE A-C 3.8928571 62.37688 FALSE A-D 0.4821429 62.37688 FALSE . . . M-P 14.250 60.26179 FALSE N-O 1.375 60.26179 FALSE N-P 6.125 60.26179 FALSE O-P 4.750 60.26179 FALSE # database group resp A 0.1 A 581.8 A 90.5 A 70.1 A 820.1 A 1159.2 A 2478.1 A 2475.3 B 351.8 B 370.1 B 326.1 B 751.9 B 931 B 588.2 B 70.1 B 1754.9 C 289.8 C 254.1 C 370.3 C 459.8 C 412.5 C 591.5 C 986.9 C 890 D 425.6 D 397.4 D 464 D 370.9 D 417.3 D 455 D 568.2 D 599.4 E 405.1 E 626.2 E 299 E 493.8 E 362.6 E 309.8 E 522.7 E 433.3 F 698.6 F 42.5 F 7.4 F 10.6 F 95.8 F 497.5 F 987.9 F 925.1 G 492.9 G 376 G 413 G 278.3 G 344.2 G 292.2 G 429.4 G 368 H 241.6 H 230.5 H 310.4 H 372.5 H 366.1 H 307.9 H 480 H 529.8 I 296 I 288.8 I 302.1 I 300.8 I 150.1 I 381.9 I 583.1 I 489.4 J 1.2 J 18.6 J 7.7 J 11.6 J 48.1 J 121.8 J 1284.1 J 944.7 L 0.5 L 44.4 L 80.9 L 15.3 L 80 L 379.9 L 940.6 L 829.3 M 323.6 M 401.5 M 162.1 M 136.5 M 139.4 M 363.3 M 280.7 M 356.5 N 197.6 N 245.9 N 221.5 N 224.3 N 185.4 N 265.3 N 304.8 N 351.9 O 189.9 O 237.3 O 247.1 O 230.4 O 272.2 O 155.1 O 270.7 O 315.2 P 48.4 P 15.5 P 53.1 P 72.8 P 74.8 P 132.3 P 550 P 478.7 On Fri, Jul 5, 2013 at 8:06 AM, Jose Iparraguirre jose.iparragui...@ageuk.org.uk wrote: Humber, Have a look at this: http://r.789695.n4.nabble.com/Multiple-Comparisons-Kruskal-Wallis-Test-kruskal-agricolae-and-kruskalmc-pgirmess-don-t-yield-the-sa-td4639004.html Hope it helps. Kind regards, José Prof. José Iparraguirre Chief Economist Age UK -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org ] On Behalf Of Humber Andrade Sent: 05 July 2013 11:38 To: r-help@r-project.org Subject: [R] kruskal.test followed by kruskalmc Hi all, After running kruskal.test I have got results (p0,005) pointing to reject the hypothesis that the
Re: [R] Error when building a custom package
Thanks for the reply. The test package build I did using the
skeleton package was in the same folder as the one I'm using
now for the custom 'slo' package, so that is not the problem.
Have you checked the detailed logs to see exactly which file is causing the
trouble?
***
This email and any attachments are confidential. Any use...{{dropped:8}}
__
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] Obtaining predicted values from a zero-inflated poisson regression model
I've got a data set of about 5,000 observations with a zero-inflated count response variable, two predictor variables and a variable which is effectively an area of opportunity, which I want to use as an offset term (all continuous). I want to explore the association between these variables, in particular I'm interested in the association of one particular variable net of the others and it would be useful when presenting the results of the analysis to hold my other variables constant at say their means or medians, and show what my model predicts would happen to my RV as it varied. In reality there are a bunch of methods issues with this that I don't want to explore in this post, (like would a hurdle model be better, and as it is a model with spatial variables, should I use something optimised for this, bootstrapping the CIs etc) as I want to get to an adequate level of understanding of how this particular model behaves before I go any further with this. I have incorporated the code to synthesise some data that look a bit like mine, I've then gone about making some predictions the best way I can given my current level of understanding of the package I've used (pscl), with some brief annotation. #packages required require(pscl) require(emdbook) #get some data set.seed(33) dsdata.frame(a=rzinbinom(500,mu=2,zprob=0.6,size=0.8),b=runif(500,3,70),c=runif(500,0.3,10),d=runif(500,0.03,0.8)) #run the model mymodel-zeroinfl(a~b+c,offset=d,data=ds,dist=negbin,EM=TRUE) summary(mymodel) #make a prediction data set, holding c and d at their means and varying b over the range of original values mypred-data.frame(b=3:70,c=5.18,d=0.39) #generate predicted values of response variable, add to the data frame pred_rv-predict.zeroinfl(mymodel,newdata=mypred,type=response,model=count,se=TRUE) #assemble the prediction data.frame predfinal-data.frame(mypred,pred_rv) colnames(predfinal)-c(b,c,d,predicted,lower,upper,se) predfinal My questions are: What exactly does this give me? I am trying to get counts in their original scale that my model predicts. Is what I have here the modelled counts from just the count component of the model? If I change the model term to zero from count I get the same result, why is this? Or indeed have i simply completely gone the wrong way about getting these estimates? Also, I can't seem to manually replicate any of these predicted values using the model results. I've tried looking at the help page of the predict.zeroinfl and predict commands, but still cannot get much further with this. Any help is much appreciated. Regards GavinR [[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] Substituting Greek symbols in some tick labels
On Jul 4, 2013, at 8:14 PM, Eric Archer - NOAA Federal wrote: I have a character vector that I'm using to label ticks in a dotchart. Some of the elements in the vector have an asterisk (*) where a Greek Delta needs to be placed when the plot is generated. Here's a simple example: x - 1:4 x.lab - c(a*a, bbb, c*c, ddd) dotchart(x, labels = x.lab) The first and third labels should be 'aDeltaa' and 'cDeltac'. I've tried things like, x.lab - strsplit(x.lab, [*]) x.lab - lapply(x.lab, function(y) expression(paste(y, sep = Delta))) but because 'y' is unevaluated, the resulting list elements won't work as tick labels. I've tried to modify it by using bquote and substitute, but couldn't get anything closer. Any suggestions? Thanks! Right. I would not throw away the * but rather change it to an inline Delta: x.lab - gsub(\\*, *Delta*, x.lab) # parse(text=x.lab) # returns: expression(a*Delta*a, bbb, c*Delta*c, ddd)# Then parse it: dotchart(x, labels = parse(text=x.lab) ) Cheers, eric -- Eric Archer, Ph.D. Southwest Fisheries Science Center NMFS, NOAA 8901 La Jolla Shores Drive La Jolla, CA 92037 USA 858-546-7121 (work) 858-546-7003 (FAX) Marine Mammal Genetics Group: swfsc.noaa.gov/prd-mmgenetics ETP Cetacean Assessment Program: swfsc.noaa.gov/prd-etp The universe doesn't care what you believe. The wonderful thing about science is that it doesn't ask for your faith, it just asks for your eyes. - Randall Munroe Lighthouses are more helpful than churches. - Benjamin Franklin ...but I'll take a GPS over either one. - John C. Craig George [[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. 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] Substituting Greek symbols in some tick labels
David, That's perfect! I just didn't think to use 'parse'. Thanks! Cheers, eric On Fri, Jul 5, 2013 at 8:20 AM, David Winsemius dwinsem...@comcast.netwrote: On Jul 4, 2013, at 8:14 PM, Eric Archer - NOAA Federal wrote: I have a character vector that I'm using to label ticks in a dotchart. Some of the elements in the vector have an asterisk (*) where a Greek Delta needs to be placed when the plot is generated. Here's a simple example: x - 1:4 x.lab - c(a*a, bbb, c*c, ddd) dotchart(x, labels = x.lab) The first and third labels should be 'aDeltaa' and 'cDeltac'. I've tried things like, x.lab - strsplit(x.lab, [*]) x.lab - lapply(x.lab, function(y) expression(paste(y, sep = Delta))) but because 'y' is unevaluated, the resulting list elements won't work as tick labels. I've tried to modify it by using bquote and substitute, but couldn't get anything closer. Any suggestions? Thanks! Right. I would not throw away the * but rather change it to an inline Delta: x.lab - gsub(\\*, *Delta*, x.lab) # parse(text=x.lab) # returns: expression(a*Delta*a, bbb, c*Delta*c, ddd)# Then parse it: dotchart(x, labels = parse(text=x.lab) ) Cheers, eric -- Eric Archer, Ph.D. Southwest Fisheries Science Center NMFS, NOAA 8901 La Jolla Shores Drive La Jolla, CA 92037 USA 858-546-7121 (work) 858-546-7003 (FAX) Marine Mammal Genetics Group: swfsc.noaa.gov/prd-mmgenetics ETP Cetacean Assessment Program: swfsc.noaa.gov/prd-etp The universe doesn't care what you believe. The wonderful thing about science is that it doesn't ask for your faith, it just asks for your eyes. - Randall Munroe Lighthouses are more helpful than churches. - Benjamin Franklin ...but I'll take a GPS over either one. - John C. Craig George [[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. David Winsemius Alameda, CA, USA -- Eric Archer, Ph.D. Southwest Fisheries Science Center NMFS, NOAA 8901 La Jolla Shores Drive La Jolla, CA 92037 USA 858-546-7121 (work) 858-546-7003 (FAX) Marine Mammal Genetics Group: swfsc.noaa.gov/prd-mmgenetics ETP Cetacean Assessment Program: swfsc.noaa.gov/prd-etp The universe doesn't care what you believe. The wonderful thing about science is that it doesn't ask for your faith, it just asks for your eyes. - Randall Munroe Lighthouses are more helpful than churches. - Benjamin Franklin ...but I'll take a GPS over either one. - John C. Craig George [[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] problem on reading many files
On Jul 5, 2013, at 1:04 AM, Pascal Oettli wrote: Hello, I'm not sure whether it is dur to the HTML version or not, but there is a problem with quotation marks in your first line. I'm guessing the OP was using Word as an editor, a practice that is doomed to frustration. -- David. Regards, Pascal 2013/7/5 RODRIGUEZ MORENO VICTOR MANUEL rodriguez.vic...@inifap.gob.mx Hi, I have to run almost 120 stations files of temperatura (mx and min), separated, and rain. I am following the instructions on RHtests tutorial as on http://www.cmc.org.ve/mediawiki/index.php?title=Preparando_los_datoslink. Bit I have no success on running on multiple files. I have the .ls file which include the names of the file station, which I include on this msg for your consideration. The instruction that I wrote on the console is listatmn - readLines („tmaxcopia.ls‰, warn=false) for(ifile in listatmn) FindU(paste(./,ifile,sep=),paste(salidas/,ifile,sep=),-999.9) Of course the tmaxcopia.ls is the list of stations for tmax data. But once running on the console, this is the error msg I got listatmn - readLines (Œtmaxcopia.ls‚, warn=false) Error: inesperado entrada in listatmn - readLines (Œ for(ifile in listatmn) FindU(paste(./,ifile,sep=),paste(salidas/,ifile,sep=),-999.9) Error: objeto 'ifnames' no encontrado What I'd like to do is to make easier and fast running the RHtest Thank's a lot in advance for your help. Dr Victor M Rodríguez M Doctor en Ciencias, en Ciencias de la Tierra / Geociencias Ambientales INIFAP. CEPAB. (465) 9580167 y 86 Ext 108 y 220 / Este msg y su(s) atado(s) tiene carácter de CONFIDENCIAL y solo es para los fines descritos en su contenido. Queda expresamente prohibida bajo la ley de protección de datos del Gobierno de los Estados Unidos Mexicanos y para su difusion o extensión a terceros, se requiere del consentimiento expresso del titular de la cuenta de correo electrónico y/o del administrador del sitio en extenso de sus responsabilidades . / [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. 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] R and MatLab implementations of the same model differs
Thank you very much to all of you for your help. This model now works as it
should (I believe). This is the final code:
rm(list=ls())
library(deSolve)
ST - function(time, init, parms) {
with(as.list(c(init, parms)),{
#functions to calculate activation m and inactivation h of the currents
#Axon
mNax - function(v) 1/(1+exp(-(v+24.7)/5.29))
taumNa - function(v) 1.32-(1.26/(1+exp(-(v+120)/25)))
hNax - function(v) 1/(1+exp((v+48.9)/5.18))
tauhNa - function(v)
0.67/(1+exp(-(v+62.9)/10))*(1.5+(1/(1+exp((v+34.9)/3.6
mKx - function(v) 1/(1+exp(-(v+14.2)/11.8))
taumK - function(v) 7.2-(6.4/(1+exp(-(v+28.3)/19.2)))
# Currents as product of maximal conducatance(g), activation(m) and
inactivation(h)
# Driving force (v-E) where E is the reversal potential of the
particular ion
# AB axon
iNa_axon_AB -
gNa_axon_AB*mNa_axon_AB^3*hNa_axon_AB*(v_axon_AB-ENa_axon_AB)
iK_axon_AB - gK_axon_AB*mK_axon_AB^4*(v_axon_AB-EK_axon_AB)
iLeak_axon_AB - gLeak_axon_AB*(v_axon_AB-ELeak_axon_AB)
# AB Axon equations
dv_axon_AB - (0 -iNa_axon_AB -iK_axon_AB -iLeak_axon_AB)/C_axon_AB
dmNa_axon_AB - (mNax(v_axon_AB)-mNa_axon_AB)/taumNa(v_axon_AB)
dhNa_axon_AB - (hNax(v_axon_AB)-hNa_axon_AB)/tauhNa(v_axon_AB)
dmK_axon_AB - (mKx(v_axon_AB)-mK_axon_AB)/taumK(v_axon_AB)
list(c(dv_axon_AB,dmNa_axon_AB, dhNa_axon_AB, dmK_axon_AB
))
})}
## Set initial state
init = c(v_axon_AB=-55,mNa_axon_AB=0,hNa_axon_AB=1,mK_axon_AB=0)
## Set parameters
# AB
gNa_axon_AB=300e-3
gK_axon_AB=52.5e-3
gLeak_axon_AB=0.0018e-3
# AB Axon Reversal potentials
ENa_axon_AB=50
EK_axon_AB=-80
ELeak_axon_AB=-60
C_axon_AB=1.5e-3
I=0
parms =
c(ENa_axon_AB,EK_axon_AB,ELeak_axon_AB,gNa_axon_AB,gK_axon_AB,gLeak_axon_AB,C_axon_AB)
## Set integrations times
times = seq(from=0, to=5000, by = 0.05);
out-ode(y=init, times=times, func=ST, parms=parms, method=ode45)
o-data.frame(out)
plot(o$v_axon_AB~o$time,type='l')
My MatLab expert came to my rescue and both model now produce similar
results :-)
Regards
Jannetta
On 5 July 2013 12:11, Jannetta Steyn janne...@henning.org wrote:
On 5 July 2013 09:44, Berend Hasselman b...@xs4all.nl wrote:
On 05-07-2013, at 09:53, Jannetta Steyn janne...@henning.org wrote:
I don't quite know how to explain the doesn't work in more detail
without
any visual aid.
You said that R got into an indefinite loop, whatever that maybe.
When I change the solver to ode45 the script never stops running. I
have even left it over night at one stage but the next day it was still
busy. Is there a way to see what it is doing and to determine why it seems
to be in an
infinite loop?
The script just ran using ode45!! For the first time ever.
Please keep the conversation on R-help.
Don't reply to me personally.
Apologies for that. It was meant to go to the list.
Which script?
With the corrections suggested by me and Peter? (gK_axon_AB=52.5-3 ==
gK_axon_AB=52.5e-3)
The R script. With the correction of AB=52.5e-3
For what it's worth: lsdoda seems quicker.
Yes, lsoda is quicker but I want to use ode45 because that is what they
used in the paper I am referencing and I am trying to get the same results
as they did.
Variable v_acon_AB now converges to -60 (the equilibrium state?)
I'm trying to determine now whether this is correct.
Regards
Jannetta
--
===
Web site: http://www.jannetta.com
Email: janne...@henning.org
===
--
===
Web site: http://www.jannetta.com
Email: janne...@henning.org
===
[[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] R and MatLab implementations of the same model differs
On 05-07-2013, at 17:37, Jannetta Steyn janne...@henning.org wrote:
Thank you very much to all of you for your help. This model now works as it
should (I believe). This is the final code:
rm(list=ls())
library(deSolve)
ST - function(time, init, parms) {
with(as.list(c(init, parms)),{
#functions to calculate activation m and inactivation h of the currents
#Axon
mNax - function(v) 1/(1+exp(-(v+24.7)/5.29))
taumNa - function(v) 1.32-(1.26/(1+exp(-(v+120)/25)))
hNax - function(v) 1/(1+exp((v+48.9)/5.18))
tauhNa - function(v)
0.67/(1+exp(-(v+62.9)/10))*(1.5+(1/(1+exp((v+34.9)/3.6
mKx - function(v) 1/(1+exp(-(v+14.2)/11.8))
taumK - function(v) 7.2-(6.4/(1+exp(-(v+28.3)/19.2)))
# Currents as product of maximal conducatance(g), activation(m) and
inactivation(h)
# Driving force (v-E) where E is the reversal potential of the
particular ion
# AB axon
iNa_axon_AB -
gNa_axon_AB*mNa_axon_AB^3*hNa_axon_AB*(v_axon_AB-ENa_axon_AB)
iK_axon_AB - gK_axon_AB*mK_axon_AB^4*(v_axon_AB-EK_axon_AB)
iLeak_axon_AB - gLeak_axon_AB*(v_axon_AB-ELeak_axon_AB)
# AB Axon equations
dv_axon_AB - (0 -iNa_axon_AB -iK_axon_AB -iLeak_axon_AB)/C_axon_AB
dmNa_axon_AB - (mNax(v_axon_AB)-mNa_axon_AB)/taumNa(v_axon_AB)
dhNa_axon_AB - (hNax(v_axon_AB)-hNa_axon_AB)/tauhNa(v_axon_AB)
dmK_axon_AB - (mKx(v_axon_AB)-mK_axon_AB)/taumK(v_axon_AB)
list(c(dv_axon_AB,dmNa_axon_AB, dhNa_axon_AB, dmK_axon_AB
))
})}
## Set initial state
init = c(v_axon_AB=-55,mNa_axon_AB=0,hNa_axon_AB=1,mK_axon_AB=0)
## Set parameters
# AB
gNa_axon_AB=300e-3
gK_axon_AB=52.5e-3
gLeak_axon_AB=0.0018e-3
# AB Axon Reversal potentials
ENa_axon_AB=50
EK_axon_AB=-80
ELeak_axon_AB=-60
C_axon_AB=1.5e-3
I=0
parms =
c(ENa_axon_AB,EK_axon_AB,ELeak_axon_AB,gNa_axon_AB,gK_axon_AB,gLeak_axon_AB,C_axon_AB)
## Set integrations times
times = seq(from=0, to=5000, by = 0.05);
out-ode(y=init, times=times, func=ST, parms=parms, method=ode45)
o-data.frame(out)
plot(o$v_axon_AB~o$time,type='l')
You must make the parms vector a named vector because the with(.. in your
function ST is not using the global values for those entities available in the
global environment.
You can check this by inserting this line after the line with times = and
before the line out - ode(…
rm(list=c(ENa_axon_AB,EK_axon_AB,ELeak_axon_AB,gNa_axon_AB,gK_axon_AB,gLeak_axon_AB,C_axon_AB))
You will get an error.
So create the parms vector like this
parms =
c(ENa_axon_AB=ENa_axon_AB,EK_axon_AB=EK_axon_AB,ELeak_axon_AB=ELeak_axon_AB,gNa_axon_AB=gNa_axon_AB,
gK_axon_AB=gK_axon_AB,gLeak_axon_AB=gLeak_axon_AB,C_axon_AB=C_axon_AB)
or better
parms - c(ENa_axon_AB=50,EK_axon_AB=-80,ELeak_axon_AB=-60,gNa_axon_AB=300e-3,
gK_axon_AB=52.5e-3,gLeak_axon_AB=0.0018e-3,C_axon_AB=1.5e-3)
and remove the expressions with the separate assignments to the variables.
And use - !!
Berend
__
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 and MatLab implementations of the same model differs
Aah. Thanks. I remember someone mentioned a named vector and I meant to ask
what that was, but I forgot. I have now fixed it.
On 5 July 2013 17:18, Berend Hasselman b...@xs4all.nl wrote:
On 05-07-2013, at 17:37, Jannetta Steyn janne...@henning.org wrote:
Thank you very much to all of you for your help. This model now works as
it
should (I believe). This is the final code:
rm(list=ls())
library(deSolve)
ST - function(time, init, parms) {
with(as.list(c(init, parms)),{
#functions to calculate activation m and inactivation h of the
currents
#Axon
mNax - function(v) 1/(1+exp(-(v+24.7)/5.29))
taumNa - function(v) 1.32-(1.26/(1+exp(-(v+120)/25)))
hNax - function(v) 1/(1+exp((v+48.9)/5.18))
tauhNa - function(v)
0.67/(1+exp(-(v+62.9)/10))*(1.5+(1/(1+exp((v+34.9)/3.6
mKx - function(v) 1/(1+exp(-(v+14.2)/11.8))
taumK - function(v) 7.2-(6.4/(1+exp(-(v+28.3)/19.2)))
# Currents as product of maximal conducatance(g), activation(m) and
inactivation(h)
# Driving force (v-E) where E is the reversal potential of the
particular ion
# AB axon
iNa_axon_AB -
gNa_axon_AB*mNa_axon_AB^3*hNa_axon_AB*(v_axon_AB-ENa_axon_AB)
iK_axon_AB - gK_axon_AB*mK_axon_AB^4*(v_axon_AB-EK_axon_AB)
iLeak_axon_AB - gLeak_axon_AB*(v_axon_AB-ELeak_axon_AB)
# AB Axon equations
dv_axon_AB - (0 -iNa_axon_AB -iK_axon_AB -iLeak_axon_AB)/C_axon_AB
dmNa_axon_AB - (mNax(v_axon_AB)-mNa_axon_AB)/taumNa(v_axon_AB)
dhNa_axon_AB - (hNax(v_axon_AB)-hNa_axon_AB)/tauhNa(v_axon_AB)
dmK_axon_AB - (mKx(v_axon_AB)-mK_axon_AB)/taumK(v_axon_AB)
list(c(dv_axon_AB,dmNa_axon_AB, dhNa_axon_AB, dmK_axon_AB
))
})}
## Set initial state
init = c(v_axon_AB=-55,mNa_axon_AB=0,hNa_axon_AB=1,mK_axon_AB=0)
## Set parameters
# AB
gNa_axon_AB=300e-3
gK_axon_AB=52.5e-3
gLeak_axon_AB=0.0018e-3
# AB Axon Reversal potentials
ENa_axon_AB=50
EK_axon_AB=-80
ELeak_axon_AB=-60
C_axon_AB=1.5e-3
I=0
parms =
c(ENa_axon_AB,EK_axon_AB,ELeak_axon_AB,gNa_axon_AB,gK_axon_AB,gLeak_axon_AB,C_axon_AB)
## Set integrations times
times = seq(from=0, to=5000, by = 0.05);
out-ode(y=init, times=times, func=ST, parms=parms, method=ode45)
o-data.frame(out)
plot(o$v_axon_AB~o$time,type='l')
You must make the parms vector a named vector because the with(.. in your
function ST is not using the global values for those entities available in
the global environment.
You can check this by inserting this line after the line with times = and
before the line out - ode(
rm(list=c(ENa_axon_AB,EK_axon_AB,ELeak_axon_AB,gNa_axon_AB,gK_axon_AB,gLeak_axon_AB,C_axon_AB))
You will get an error.
So create the parms vector like this
parms =
c(ENa_axon_AB=ENa_axon_AB,EK_axon_AB=EK_axon_AB,ELeak_axon_AB=ELeak_axon_AB,gNa_axon_AB=gNa_axon_AB,
gK_axon_AB=gK_axon_AB,gLeak_axon_AB=gLeak_axon_AB,C_axon_AB=C_axon_AB)
or better
parms -
c(ENa_axon_AB=50,EK_axon_AB=-80,ELeak_axon_AB=-60,gNa_axon_AB=300e-3,
gK_axon_AB=52.5e-3,gLeak_axon_AB=0.0018e-3,C_axon_AB=1.5e-3)
and remove the expressions with the separate assignments to the variables.
And use - !!
Berend
--
===
Web site: http://www.jannetta.com
Email: janne...@henning.org
===
[[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] multcomp on significant interaction in coxme model
Dear R community I currently try to get post hoc multiple comparisons with the package multcomp from a cox mixed-effects model, where the survival is explained by two variables (cover with levels: nocover and cover; treatment with levels: tx, uv, meta), whose interaction is significant. I read Hothorn, T. 2011: Additional multcomp Examples and there is an example involving a two-way ANOVA with significant interaction. Using this guideline I wanted to compare the levels of treatment in my data within the levels of cover. My model is as follows: m.lapf-coxme(Surv(day,status)~cover*treatment+(1|r.brood)+(1|r.worker)+(1|species/brood),data=removal.lapf) First I constructed a contrast matrix based on Tukey-contrasts as follows: Tukey-contrMat(table(removal.lapf$treatment), Tukey) Tukey K1-cbind(Tukey,matrix(0,nrow=nrow(Tukey),ncol=ncol(Tukey))) K1 rownames(K1)-paste(levels(removal.lapf$cover)[1],rownames(K1),sep=:) K1 K2-cbind(matrix(0,nrow=nrow(Tukey),ncol=ncol(Tukey)),Tukey) K2 rownames(K2)-paste(levels(removal.lapf$cover)[2],rownames(K1),sep=:) K2 K-rbind(K1,K2) colnames(K)-c(colnames(Tukey),colnames(Tukey)) K This gives the following matrix meta tx uv meta tx uv cocooned:tx - meta -1 1 0 0 0 0 cocooned:uv - meta -10 1 0 0 0 cocooned:uv - tx0 -1 1 0 0 0 naked:cocooned:tx - meta00 0 -11 0 naked:cocooned:uv - meta0 0 0 -1 0 1 naked:cocooned:uv - tx 00 0 0-1 1 then I constructed a secon matrix tmp-expand.grid(treatment=unique(removal.lapf$treatment),cover=unique(removal.lapf$cover)) X-model.matrix(~cover*treatment,data=tmp) X which gives the following (Intercept) covernaked treatmenttx treatmentuv covernaked:treatmenttx covernaked:treatmentuv 1 1 1 1 0 1 0 2 1 1 0 1 0 1 3 1 1 0 0 0 0 4 1 0 1 0 0 0 5 1 0 0 1 0 0 6 1 0 0 0 0 0 attr(,assign) [1] 0 1 2 2 3 3 attr(,contrasts) attr(,contrasts)$cover [1] contr.treatment attr(,contrasts)$treatment [1] contr.treatment however when performing the tests via summary(glht(m.lapf,linfct=K%*%X)) I get the following error message Error in summary(glht(m.lapf, linfct = K %*% X)) : error in evaluating the argument 'object' in selecting a method for function 'summary': Error in glht.matrix(m.lapf, linfct = K %*% X) : 'ncol(linfct)' is not equal to 'length(coef(model))' Performing exactly the same routine with the same data on a logistic model with family=binomial does not give this error message. So my question is, what am I missing here? Thanks, for any possible input -- Simon Tragust Animal Ecology I University Bayreuth D-95440 Bayreuth +49 921 552464 [[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] kruskal.test followed by kruskalmc
On Jul 5, 2013, at 16:34 , Humber Andrade wrote: Thank you Dr. Dalgaard. I understood. I know that this list is not to discuss statistics but I would be very glad if you or someone else can give me some opinion on how to proceed. The kruskal.test says there are differences but the multiple comparisons do not point out what are the differences. Can you suggest a suitable way (maybe paired wilcoxon) to infer what are the differences? I am asking for the hint because I am sure the journal editor/reviewer will ask me to point out which groups differ from each other. I'm afraid it can't be done. You really can be in a situation where you reject the global null hypothesis that all groups are the same, yet cannot point out any two groups that differ from eachother. -pd Regards, Humber On Fri, Jul 5, 2013 at 11:04 AM, peter dalgaard pda...@gmail.com wrote: On Jul 5, 2013, at 15:00 , Humber Andrade wrote: Thank you Prof. José Iparraguirre. Maybe I am wrong but I think the issues are not the same. His data doesn't showed significant differences after kruskal.test(), that was not my case. Anyway follow below my the results I've got and the database. This can happen. It is a matter of probability theory, not of R. The following is a simplified paraphrase of what is going on: # 15 random normals, compare range test to variance test # Simulate everything for simplicity # Null distribution M0 - replicate(1, rnorm(15)) vv - apply(M0,2,var) rg - apply(M0,2,range) rr - apply(rg,2,diff) r.95 - quantile(rr, .95) v.95 - quantile(vv, .95) v.995 - quantile(vv, .995) # Distribution at quadrupled variance M - replicate(1, rnorm(15,sd=2)) vv - apply(M,2,var) rg - apply(M,2,range) rr - apply(rg,2,diff) plot(rr,vv) abline(h=c(v.95,v.995),v=r.95, col=red) Notice that the two statistics are correlated, but not equivalent. There are cases where one value is beyond the .95 level and the other is not. Since it is theoretically optimal to use the variance as the test statistic in this model, there are quite a few more cases where rr is below the cutoff and vv is above than the other way around. There are even a sizable number of cases where vv is beyond v.995 and rr does not reach r.95. (The theoretical optimality applies only because I use an increased variance alternative. For specific alternatives, the picture can change. Try it, for instance with a single mean substantially different from the others: M - replicate(1, rnorm(15,mean=rep(c(4,0),c(1,14 ) Thank you, # kruskal.test(data$resp,data$group) Kruskal-Wallis rank sum test data: data$resp and data$group Kruskal-Wallis chi-squared = 32.3546, df = 14, p-value = 0.003566 kruskalmc(data$resp,data$group) Multiple comparison test after Kruskal-Wallis p.value: 0.05 Comparisons obs.dif critical.dif difference A-B 4.8303571 62.37688 FALSE A-C 3.8928571 62.37688 FALSE A-D 0.4821429 62.37688 FALSE . . . M-P 14.250 60.26179 FALSE N-O 1.375 60.26179 FALSE N-P 6.125 60.26179 FALSE O-P 4.750 60.26179 FALSE # database group resp A 0.1 A 581.8 A 90.5 A 70.1 A 820.1 A 1159.2 A 2478.1 A 2475.3 B 351.8 B 370.1 B 326.1 B 751.9 B 931 B 588.2 B 70.1 B 1754.9 C 289.8 C 254.1 C 370.3 C 459.8 C 412.5 C 591.5 C 986.9 C 890 D 425.6 D 397.4 D 464 D 370.9 D 417.3 D 455 D 568.2 D 599.4 E 405.1 E 626.2 E 299 E 493.8 E 362.6 E 309.8 E 522.7 E 433.3 F 698.6 F 42.5 F 7.4 F 10.6 F 95.8 F 497.5 F 987.9 F 925.1 G 492.9 G 376 G 413 G 278.3 G 344.2 G 292.2 G 429.4 G 368 H 241.6 H 230.5 H 310.4 H 372.5 H 366.1 H 307.9 H 480 H 529.8 I 296 I 288.8 I 302.1 I 300.8 I 150.1 I 381.9 I 583.1 I 489.4 J 1.2 J 18.6 J 7.7 J 11.6 J 48.1 J 121.8 J 1284.1 J 944.7 L 0.5 L 44.4 L 80.9 L 15.3 L 80 L 379.9 L 940.6 L 829.3 M 323.6 M 401.5 M 162.1 M 136.5 M 139.4 M 363.3 M 280.7 M 356.5 N 197.6 N 245.9 N 221.5 N 224.3 N 185.4 N 265.3 N 304.8 N 351.9 O 189.9 O 237.3 O 247.1 O 230.4 O 272.2 O 155.1 O 270.7 O 315.2 P 48.4 P 15.5 P 53.1 P 72.8 P 74.8 P 132.3 P 550 P 478.7 On Fri, Jul 5, 2013 at 8:06 AM, Jose Iparraguirre jose.iparragui...@ageuk.org.uk wrote: Humber, Have a look at this: http://r.789695.n4.nabble.com/Multiple-Comparisons-Kruskal-Wallis-Test-kruskal-agricolae-and-kruskalmc-pgirmess-don-t-yield-the-sa-td4639004.html Hope it helps. Kind regards, José Prof. José Iparraguirre Chief Economist Age UK -Original Message- From: r-help-boun...@r-project.org
Re: [R] Lattice barchart with error bars
Hi Bert, Dennis, I'll agree that using a barchart was a poor choice. I was in fact using a notched bwplot to show the median and confidence interval of the median. In this case it's the median and confidence interval that I want to highlight, and I find that the visual noise of the box and whiskers is detracting from the focus, and those wee notches are not much to focus on. So, I'd like to draw a stripplot with error bars, preferably in Lattice. Let's call this a TIE fighter plot. Any suggestions? Cheers, Shaun On 4 July 2013 18:00, Dennis Murphy djmu...@gmail.com wrote: If you consult the lattice package help, you'll discover there is no panel_errorbar() function, which would imply the package developers have a distaste for that type of graphic. If you fish around the R-help archives, though, you might be able to find someone who wrote a function to do error bars in lattice. (Use a searchable archive such as Nabble to hunt for it.) Error bar plots are easier to do in the ggplot2 package, since there is a specific function to generate the error bar 'geometry' (geom_errorbar). See http://docs.ggplot2.org/current/ for an expanded version of the package help pages, which include the graphs generated by the code. I believe there's also a base graphics version that you can get from the gplots package, but I don't know a lot about it. Dennis On Thu, Jul 4, 2013 at 2:53 PM, Shaun Jackman sjack...@gmail.com wrote: Hi, I'd like to draw a lattice barchart of means with error bars to show the standard deviation. I have the barchart, how do I add the error bars? require(datasets) require(lattice) x - aggregate(weight ~ Diet, ChickWeight, function(x) c(mean=mean(x), sd=sd(x))) barchart(weight[,'mean'] ~ Diet, x) Thanks, Shaun __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Lattice barchart with error bars
On Jul 5, 2013, at 11:15 AM, Shaun Jackman wrote: Hi Bert, Dennis, I'll agree that using a barchart was a poor choice. I was in fact using a notched bwplot to show the median and confidence interval of the median. In this case it's the median and confidence interval that I want to highlight, and I find that the visual noise of the box and whiskers is detracting from the focus, and those wee notches are not much to focus on. So, I'd like to draw a stripplot with error bars, preferably in Lattice. Let's call this a TIE fighter plot. Any suggestions? I like the TIE fighter label. Try this: library(latticeExtra) data(USCancerRates) segplot(reorder(factor(county), rate.male) ~ LCL95.male + UCL95.male, data = subset(USCancerRates, state == Washington), draw.bands = FALSE, centers = rate.male, segments.fun = panel.arrows, ends = both, angle = 90, length = 1, unit = mm) It's what Sarkar has recommended in the past when this request has been posted. -- David Cheers, Shaun On 4 July 2013 18:00, Dennis Murphy djmu...@gmail.com wrote: If you consult the lattice package help, you'll discover there is no panel_errorbar() function, which would imply the package developers have a distaste for that type of graphic. If you fish around the R-help archives, though, you might be able to find someone who wrote a function to do error bars in lattice. (Use a searchable archive such as Nabble to hunt for it.) Error bar plots are easier to do in the ggplot2 package, since there is a specific function to generate the error bar 'geometry' (geom_errorbar). See http://docs.ggplot2.org/current/ for an expanded version of the package help pages, which include the graphs generated by the code. I believe there's also a base graphics version that you can get from the gplots package, but I don't know a lot about it. Dennis On Thu, Jul 4, 2013 at 2:53 PM, Shaun Jackman sjack...@gmail.com wrote: Hi, I'd like to draw a lattice barchart of means with error bars to show the standard deviation. I have the barchart, how do I add the error bars? require(datasets) require(lattice) x - aggregate(weight ~ Diet, ChickWeight, function(x) c(mean=mean(x), sd=sd(x))) barchart(weight[,'mean'] ~ Diet, x) Thanks, Shaun __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. 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] Lattice barchart with error bars
Be careful! You are talking about 2 different varieties of apples here. As I read it, the CI's in the cancer data, which I know is just for example purposes, are CI's for the **individual means**; the notches in boxplots are nonparametric and for 2 groups with roughly equal sample sizes, The idea appears to be to give roughly a 95% confidence interval for the **difference** in two medians. (from ?boxplot.stats). So I'm not sure which you want, but they are certainly different (by a factor of around sqrt(2),right?), even if both are for the mean or both are for the median. Cheers, Bert On Fri, Jul 5, 2013 at 11:28 AM, David Winsemius dwinsem...@comcast.net wrote: On Jul 5, 2013, at 11:15 AM, Shaun Jackman wrote: Hi Bert, Dennis, I'll agree that using a barchart was a poor choice. I was in fact using a notched bwplot to show the median and confidence interval of the median. In this case it's the median and confidence interval that I want to highlight, and I find that the visual noise of the box and whiskers is detracting from the focus, and those wee notches are not much to focus on. So, I'd like to draw a stripplot with error bars, preferably in Lattice. Let's call this a TIE fighter plot. Any suggestions? I like the TIE fighter label. Try this: library(latticeExtra) data(USCancerRates) segplot(reorder(factor(county), rate.male) ~ LCL95.male + UCL95.male, data = subset(USCancerRates, state == Washington), draw.bands = FALSE, centers = rate.male, segments.fun = panel.arrows, ends = both, angle = 90, length = 1, unit = mm) It's what Sarkar has recommended in the past when this request has been posted. -- David Cheers, Shaun On 4 July 2013 18:00, Dennis Murphy djmu...@gmail.com wrote: If you consult the lattice package help, you'll discover there is no panel_errorbar() function, which would imply the package developers have a distaste for that type of graphic. If you fish around the R-help archives, though, you might be able to find someone who wrote a function to do error bars in lattice. (Use a searchable archive such as Nabble to hunt for it.) Error bar plots are easier to do in the ggplot2 package, since there is a specific function to generate the error bar 'geometry' (geom_errorbar). See http://docs.ggplot2.org/current/ for an expanded version of the package help pages, which include the graphs generated by the code. I believe there's also a base graphics version that you can get from the gplots package, but I don't know a lot about it. Dennis On Thu, Jul 4, 2013 at 2:53 PM, Shaun Jackman sjack...@gmail.com wrote: Hi, I'd like to draw a lattice barchart of means with error bars to show the standard deviation. I have the barchart, how do I add the error bars? require(datasets) require(lattice) x - aggregate(weight ~ Diet, ChickWeight, function(x) c(mean=mean(x), sd=sd(x))) barchart(weight[,'mean'] ~ Diet, x) Thanks, Shaun __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius Alameda, CA, USA -- 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] Lattice barchart with error bars
Yes! Thank you, David. That's exactly what I'm I'm looking for. For the record, here's a couple pages leading to this answer: http://www.hep.by/gnu/r-patched/r-faq/R-FAQ_89.html http://latticeextra.r-forge.r-project.org/man/segplot.html http://rgm3.lab.nig.ac.jp/RGM/r_function?p=latticeExtraf=segplot For a related question, what's the tidiest way to calculate the medians and confidence intervals? Currently I'm using `boxplot`: require(datasets) ci - with(boxplot(weight ~ Diet, ChickWeight), data.frame( Diet = names, median = stats[3,], lower = conf[1,], upper = conf[2,])) Cheers, Shaun On 5 July 2013 11:28, David Winsemius dwinsem...@comcast.net wrote: On Jul 5, 2013, at 11:15 AM, Shaun Jackman wrote: Hi Bert, Dennis, I'll agree that using a barchart was a poor choice. I was in fact using a notched bwplot to show the median and confidence interval of the median. In this case it's the median and confidence interval that I want to highlight, and I find that the visual noise of the box and whiskers is detracting from the focus, and those wee notches are not much to focus on. So, I'd like to draw a stripplot with error bars, preferably in Lattice. Let's call this a TIE fighter plot. Any suggestions? I like the TIE fighter label. Try this: library(latticeExtra) data(USCancerRates) segplot(reorder(factor(county), rate.male) ~ LCL95.male + UCL95.male, data = subset(USCancerRates, state == Washington), draw.bands = FALSE, centers = rate.male, segments.fun = panel.arrows, ends = both, angle = 90, length = 1, unit = mm) It's what Sarkar has recommended in the past when this request has been posted. -- David Cheers, Shaun On 4 July 2013 18:00, Dennis Murphy djmu...@gmail.com wrote: If you consult the lattice package help, you'll discover there is no panel_errorbar() function, which would imply the package developers have a distaste for that type of graphic. If you fish around the R-help archives, though, you might be able to find someone who wrote a function to do error bars in lattice. (Use a searchable archive such as Nabble to hunt for it.) Error bar plots are easier to do in the ggplot2 package, since there is a specific function to generate the error bar 'geometry' (geom_errorbar). See http://docs.ggplot2.org/current/ for an expanded version of the package help pages, which include the graphs generated by the code. I believe there's also a base graphics version that you can get from the gplots package, but I don't know a lot about it. Dennis On Thu, Jul 4, 2013 at 2:53 PM, Shaun Jackman sjack...@gmail.com wrote: Hi, I'd like to draw a lattice barchart of means with error bars to show the standard deviation. I have the barchart, how do I add the error bars? require(datasets) require(lattice) x - aggregate(weight ~ Diet, ChickWeight, function(x) c(mean=mean(x), sd=sd(x))) barchart(weight[,'mean'] ~ Diet, x) Thanks, Shaun __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. 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] Lattice barchart with error bars
Hmm. Interesting point, Bert. I don't know whether the notches show the 95% confidence interval or the median, or the 95% confidence interval that two non-overlapping notches have different medians. You're saying it's the latter? Anyone know what the 95% confidence interval of the median would be? Cheers, Shaun The notches (if requested) extend to +/-1.58 IQR/sqrt(n). This seems to be based on the same calculations as the formula with 1.57 in Chambers et al. (1983, p. 62), given in McGill et al. (1978, p. 16). They are based on asymptotic normality of the median and roughly equal sample sizes for the two medians being compared, and are said to be rather insensitive to the underlying distributions of the samples. The idea appears to be to give roughly a 95% confidence interval for the difference in two medians. On 5 July 2013 11:48, Bert Gunter gunter.ber...@gene.com wrote: Be careful! You are talking about 2 different varieties of apples here. As I read it, the CI's in the cancer data, which I know is just for example purposes, are CI's for the **individual means**; the notches in boxplots are nonparametric and for 2 groups with roughly equal sample sizes, The idea appears to be to give roughly a 95% confidence interval for the **difference** in two medians. (from ?boxplot.stats). So I'm not sure which you want, but they are certainly different (by a factor of around sqrt(2),right?), even if both are for the mean or both are for the median. Cheers, Bert On Fri, Jul 5, 2013 at 11:28 AM, David Winsemius dwinsem...@comcast.net wrote: On Jul 5, 2013, at 11:15 AM, Shaun Jackman wrote: Hi Bert, Dennis, I'll agree that using a barchart was a poor choice. I was in fact using a notched bwplot to show the median and confidence interval of the median. In this case it's the median and confidence interval that I want to highlight, and I find that the visual noise of the box and whiskers is detracting from the focus, and those wee notches are not much to focus on. So, I'd like to draw a stripplot with error bars, preferably in Lattice. Let's call this a TIE fighter plot. Any suggestions? I like the TIE fighter label. Try this: library(latticeExtra) data(USCancerRates) segplot(reorder(factor(county), rate.male) ~ LCL95.male + UCL95.male, data = subset(USCancerRates, state == Washington), draw.bands = FALSE, centers = rate.male, segments.fun = panel.arrows, ends = both, angle = 90, length = 1, unit = mm) It's what Sarkar has recommended in the past when this request has been posted. -- David Cheers, Shaun On 4 July 2013 18:00, Dennis Murphy djmu...@gmail.com wrote: If you consult the lattice package help, you'll discover there is no panel_errorbar() function, which would imply the package developers have a distaste for that type of graphic. If you fish around the R-help archives, though, you might be able to find someone who wrote a function to do error bars in lattice. (Use a searchable archive such as Nabble to hunt for it.) Error bar plots are easier to do in the ggplot2 package, since there is a specific function to generate the error bar 'geometry' (geom_errorbar). See http://docs.ggplot2.org/current/ for an expanded version of the package help pages, which include the graphs generated by the code. I believe there's also a base graphics version that you can get from the gplots package, but I don't know a lot about it. Dennis On Thu, Jul 4, 2013 at 2:53 PM, Shaun Jackman sjack...@gmail.com wrote: Hi, I'd like to draw a lattice barchart of means with error bars to show the standard deviation. I have the barchart, how do I add the error bars? require(datasets) require(lattice) x - aggregate(weight ~ Diet, ChickWeight, function(x) c(mean=mean(x), sd=sd(x))) barchart(weight[,'mean'] ~ Diet, x) Thanks, Shaun __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius Alameda, CA, USA -- 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,
Re: [R] Filter Dataframe for Alarm for particular column(s).
Hi,
May be this helps:
If you had showed your solution, it would be easier to compare.
res-data.frame(lapply(sapply(MyDF[,c(2,4)],function(x)
{x1-which(c(0,diff(x))0);x1[length(x1)==0]-0;x1}),`[`,1))
res
# TNH BIX
#1 3 9
#Speed
set.seed(24)
MyDFNew-
data.frame(TNH=sample(0:1,1e6,replace=TRUE),BIX=sample(0:1,1e6,replace=TRUE))
system.time(res1-data.frame(lapply(sapply(MyDFNew,function(x)
{x1-which(c(0,diff(x))0);x1[length(x1)==0]-0;x1}),`[`,1)))
# user system elapsed
# 0.364 0.000 0.363
res1
# TNH BIX
#1 7 2
MyDFNew[1:10,]
# TNH BIX
#1 0 1
#2 0 0
#3 1 1
#4 1 1
#5 1 0
#6 1 0
#7 0 1
#8 1 1
#9 1 1
#10 0 0
A.K.
Hi,
Hi here i have a dataframe called MyDF.
a-c(1,1,1,1,1,0,0,0,1,1)
b-c(1,1,0,1,1,0,0,0,1,1)
c-c(1,1,1,1,1,1,1,0,1,1)
d-c(1,1,1,1,1,1,1,1,0,1)
MyDF-data.frame(DWATT=a,TNH=b,CSGV=c,BIX=d)
My requirement is, here i need a function - to get for a
particular row number(s), when particular column(s) value change from
one-to-zero (for the first change). Suppose there is no change is
happening then it should return Zero
For example, Using MyDF,
DWATT TNH CSGV BIX
1 1 1 1
1 1 1 1
1 0 1 1
1 1 1 1
1 1 1 1
0 0 1 1
0 0 1 1
0 0 0 1
1 1 1 0
1 1 1 1
Here i want to know, the row number where TNH-column and BIX-column values
change happening from one-to-zero for the first time.
Note:- Suppose there is no change is happening then it should return Zero
Answer should be a dataframe with single row.
So here answer should return a dataframe like this.
TNH BIX
--
3 9
i used some ways to get a solution using loops. But there is a bulk files with
bulk rows to process.
So performace is most important. Could someone please suggest better ideas ?
Thanks,
Antony.
__
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] converting a list of loglin terms to a model formula
On 7/3/2013 6:38 PM, Henrique Dallazuanna wrote:
Try this:
as.formula(sprintf( ~ %s, do.call(paste, c(lapply(mutual(3), paste,
collapse = :), sep = +
Thanks for this. I encapsulated this as a function, loglin2formula()
and a related function,
loglin2string() to give a character string representation.
They seem to work when used from the command line, but don't give
what I need when I use it in another function. I think it has something
to do with the environment
for the formula or how I pass it to MASS::loglm in my function
test_loglm (below).
I'll demonstrate the problem first, then give the functions test cases
I'm using as
runnable code.
joint(3, table=HairEyeColor)
$term1
[1] Hair Eye
$term2
[1] Sex
loglin2formula(joint(3, table=HairEyeColor))
~Hair:Eye + Sex
environment: 0x0884a640
loglin2string(joint(3, table=HairEyeColor))
[1] [Hair,Eye] [Sex]
These look like what I want, more or less; but, when I use these in the
function test_loglm
(also below) , the formula doesn't get resolved when I call loglm (it
appears as
formula = form), and the result of loglin2string gets garbled. The
numerical results are
correct; I'm concerned about the labeling in the computed loglm object.
test_loglm(HairEyeColor, type='joint')
formula:
~Hair:Eye + Sex
environment: 0x08842de0
model.string:
[1] [~] [+,Hair:Eye,Sex] ### I want [Hair,Eye]
[Sex] here
model:
Call:
loglm(formula = form, data = x) ### I want formula = ~Hair:Eye +
Sex here
Statistics:
X^2 df P( X^2)
Likelihood Ratio 19.85656 15 0.1775045
Pearson 19.56712 15 0.1891745
## --- functions and test cases, should be runnable (mod line folding)
--- ###
#' convert a loglin model to a model formula for loglm
#' @param x a list of terms in a loglinear model, such as returned by
\code{joint}, \code{conditional}, \dots
#' @source Code from Henrique Dallazuanna, www...@gmail.com, R-help
7-4-2013
loglin2formula - function(x) {
terms - lapply(x, paste, collapse = :)
as.formula(sprintf( ~ %s, do.call(paste, c(terms, sep = +
}
#' convert a loglin model to a string, using bracket notation for the
high-order terms
#' @param x a list of terms in a loglinear model, such as returned by
\code{joint}, \code{conditional}, \dots
#' @param brackets characters to use to surround model terms.
#' @param sep characters used to separate factor names within a term
#' @param collapse characters used to separate terms
#' @param abbrev
loglin2string - function(x, brackets = c('[', ']'), sep=',', collapse='
', abbrev) {
if (length(brackets)==1 (nchar(brackets)1)) brackets -
unlist(strsplit(brackets, ))
terms - lapply(x, paste, collapse=sep)
terms - paste(brackets[1], terms, brackets[2], sep='')
paste(terms, collapse= ' ')
}
#' models of joint independence, of some factors wrt one or more other
factors
#' @param nf number of factors for which to generate model
#' @param table a contingency table used for factor names, typically the
output from \code{\link[base]{table}}
#' @param factors names of factors used in the model when \code{table}
is not specified
#' @param withindices of the factors against which others are
considered jointly independent
#' @export
joint - function(nf, table=NULL, factors=1:nf, with=nf) {
if (!is.null(table)) factors - names(dimnames(table))
if (nf == 1) return (list(term1=factors[1]))
if (nf == 2) return (list(term1=factors[1], term2=factors[2]))
others - setdiff(1:nf, with)
result - list(term1=factors[others], term2=factors[with])
result
}
### Test using these in another function
test_loglm - function(
x, type = c(joint, conditional),
...)
{
nf - length(dim(x))
factors - names(dimnames(x))
type - match.arg(type)
margins - switch(type,
joint = joint(nf, factors=factors),
conditional = conditional(nf, factors=factors)
)
form - loglin2formula(margins);
cat(formula:\n); print(form)
model.string - loglin2string(form)
cat(model.string:\n); print(model.string)
mod - loglm(formula=form, data=x)
cat(model:\n); print(mod)
invisible(list(form=form, model.string=model.string, mod=mod))
}
## tests
loglin2formula(joint(3, table=HairEyeColor))
loglin2string(joint(3, table=HairEyeColor))
test_loglm(HairEyeColor, type='joint')
On Wed, Jul 3, 2013 at 6:55 PM, Michael Friendly frien...@yorku.ca
mailto:frien...@yorku.ca wrote:
I'm developing some functions to create symbolic specifications
for loglinear models of different types.
I don't really know how to 'compute' with model formulas, so I've
done this in the notation
for stats::loglin(), which is a list of high-order terms in the model.
What I'd like is a function to turn the results of these into a
model formula, suitable for
MASS::loglm. That's the reverse of what loglm does.
For example, the simplest versions of
[R] Subset and order
Hello, I have a data frame with several columns. I'd like to select some subset *and* order by another field at the same time. Example: a b c 1 2 3 3 3 4 2 4 5 1 3 4 etc I want to select all rows where b=3 and then order by a. To subset is easy: x[x$b==3,] To order is easy: x[order(x$a),] Is there a way to do both in a single efficient statement? Thanks, -- Noah Silverman, M.S., C.Phil UCLA Department of Statistics 8117 Math Sciences Building Los Angeles, CA 90095 [[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] Subset and order
Hello, Maybe like this? subset(x[order(x$a), ], b == 3) Hope this helps, Rui Barradas Em 05-07-2013 20:33, Noah Silverman escreveu: Hello, I have a data frame with several columns. I'd like to select some subset *and* order by another field at the same time. Example: a b c 1 2 3 3 3 4 2 4 5 1 3 4 etc… I want to select all rows where b=3 and then order by a. To subset is easy: x[x$b==3,] To order is easy: x[order(x$a),] Is there a way to do both in a single efficient statement? Thanks, -- Noah Silverman, M.S., C.Phil UCLA Department of Statistics 8117 Math Sciences Building Los Angeles, CA 90095 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Subset and order
That would work, but is painfully slow. It forces a new sort of the data with every query. I have 200,000 rows and need almost a hundred queries. Thanks, -N On Jul 5, 2013, at 12:43 PM, Rui Barradas ruipbarra...@sapo.pt wrote: Hello, Maybe like this? subset(x[order(x$a), ], b == 3) Hope this helps, Rui Barradas Em 05-07-2013 20:33, Noah Silverman escreveu: Hello, I have a data frame with several columns. I'd like to select some subset *and* order by another field at the same time. Example: ab c 12 3 33 4 24 5 13 4 etc… I want to select all rows where b=3 and then order by a. To subset is easy: x[x$b==3,] To order is easy: x[order(x$a),] Is there a way to do both in a single efficient statement? Thanks, -- Noah Silverman, M.S., C.Phil UCLA Department of Statistics 8117 Math Sciences Building Los Angeles, CA 90095 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Subset and order
Hello, If time is one of the problems, precompute an ordered index, and use it every time you want the df sorted. But that would mean you can't do it in a single operation. iord - order(x$a) subset(x[iord, ], b == 3) Rui Barradas Em 05-07-2013 20:47, Noah Silverman escreveu: That would work, but is painfully slow. It forces a new sort of the data with every query. I have 200,000 rows and need almost a hundred queries. Thanks, -N On Jul 5, 2013, at 12:43 PM, Rui Barradas ruipbarra...@sapo.pt wrote: Hello, Maybe like this? subset(x[order(x$a), ], b == 3) Hope this helps, Rui Barradas Em 05-07-2013 20:33, Noah Silverman escreveu: Hello, I have a data frame with several columns. I'd like to select some subset *and* order by another field at the same time. Example: a b c 1 2 3 3 3 4 2 4 5 1 3 4 etc… I want to select all rows where b=3 and then order by a. To subset is easy: x[x$b==3,] To order is easy: x[order(x$a),] Is there a way to do both in a single efficient statement? Thanks, -- Noah Silverman, M.S., C.Phil UCLA Department of Statistics 8117 Math Sciences Building Los Angeles, CA 90095 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Subset and order
It may be that single and efficient are opposing goals. Two steps lets you create the subset and then just order each query. Alternatively, if the data do not change often, create an ordered version and query that. David Carlson Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Noah Silverman Sent: Friday, July 5, 2013 2:47 PM To: Rui Barradas Cc: R-help@r-project.org Subject: Re: [R] Subset and order That would work, but is painfully slow. It forces a new sort of the data with every query. I have 200,000 rows and need almost a hundred queries. Thanks, -N On Jul 5, 2013, at 12:43 PM, Rui Barradas ruipbarra...@sapo.pt wrote: Hello, Maybe like this? subset(x[order(x$a), ], b == 3) Hope this helps, Rui Barradas Em 05-07-2013 20:33, Noah Silverman escreveu: Hello, I have a data frame with several columns. I'd like to select some subset *and* order by another field at the same time. Example: ab c 12 3 33 4 24 5 13 4 etc. I want to select all rows where b=3 and then order by a. To subset is easy: x[x$b==3,] To order is easy: x[order(x$a),] Is there a way to do both in a single efficient statement? Thanks, -- Noah Silverman, M.S., C.Phil UCLA Department of Statistics 8117 Math Sciences Building Los Angeles, CA 90095 [[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. __ 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] converting a list of loglin terms to a model formula
Call:
loglm(formula = form, data = x) ### I want formula = ~Hair:Eye + Sex
here
Since your function made the call
loglm(form, data=x)
the 'call' component of output is going to show 'form', not '~ Hair:Eye + Sex'.
You can use bquote to pre-evaluate the formula=form argument to get the call
to look nicer, as in:
form - mpg ~ wt + hp
eval(bquote(lm(.(form), data=mtcars)))
instead of
lm(form, data=mtcars)
In your loglin2formula, I would make by environment of the generated formula
the environment of the caller of loglin2formula (or somewhere else if the user
wishes)
by adding the argument
env = parent.frame()
and replacing
as.formula( sprintf(...) )
with
formula( sprintf(...), env=env)
(I don't think you've run into any problems related to having an irrelevant
environment attached to the formula, but they will happen if the formula
involves any variable names that happen to be in loglin2formula.)
I would also change the loglin2formula so it worked with non-syntactic names.
Wrapping them with backquotes would probably do it, but I may have missed
something in the back and forth between character strings and formula.
As for loglin2string, you complain that it works when given a list of character
vectors but not when is given a formula. That is not surprising. Did you mean
for test_loglm to pass it 'margins' instead of 'form'?
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 Michael Friendly
Sent: Friday, July 05, 2013 12:21 PM
To: Henrique Dallazuanna
Cc: R-help
Subject: Re: [R] converting a list of loglin terms to a model formula
On 7/3/2013 6:38 PM, Henrique Dallazuanna wrote:
Try this:
as.formula(sprintf( ~ %s, do.call(paste, c(lapply(mutual(3), paste,
collapse = :), sep = +
Thanks for this. I encapsulated this as a function, loglin2formula()
and a related function,
loglin2string() to give a character string representation.
They seem to work when used from the command line, but don't give
what I need when I use it in another function. I think it has something
to do with the environment
for the formula or how I pass it to MASS::loglm in my function
test_loglm (below).
I'll demonstrate the problem first, then give the functions test cases
I'm using as
runnable code.
joint(3, table=HairEyeColor)
$term1
[1] Hair Eye
$term2
[1] Sex
loglin2formula(joint(3, table=HairEyeColor))
~Hair:Eye + Sex
environment: 0x0884a640
loglin2string(joint(3, table=HairEyeColor))
[1] [Hair,Eye] [Sex]
These look like what I want, more or less; but, when I use these in the
function test_loglm
(also below) , the formula doesn't get resolved when I call loglm (it
appears as
formula = form), and the result of loglin2string gets garbled. The
numerical results are
correct; I'm concerned about the labeling in the computed loglm object.
test_loglm(HairEyeColor, type='joint')
formula:
~Hair:Eye + Sex
environment: 0x08842de0
model.string:
[1] [~] [+,Hair:Eye,Sex] ### I want [Hair,Eye]
[Sex] here
model:
Call:
loglm(formula = form, data = x) ### I want formula = ~Hair:Eye +
Sex here
Statistics:
X^2 df P( X^2)
Likelihood Ratio 19.85656 15 0.1775045
Pearson 19.56712 15 0.1891745
## --- functions and test cases, should be runnable (mod line folding)
--- ###
#' convert a loglin model to a model formula for loglm
#' @param x a list of terms in a loglinear model, such as returned by
\code{joint}, \code{conditional}, \dots
#' @source Code from Henrique Dallazuanna, www...@gmail.com, R-help
7-4-2013
loglin2formula - function(x) {
terms - lapply(x, paste, collapse = :)
as.formula(sprintf( ~ %s, do.call(paste, c(terms, sep = +
}
#' convert a loglin model to a string, using bracket notation for the
high-order terms
#' @param x a list of terms in a loglinear model, such as returned by
\code{joint}, \code{conditional}, \dots
#' @param brackets characters to use to surround model terms.
#' @param sep characters used to separate factor names within a term
#' @param collapse characters used to separate terms
#' @param abbrev
loglin2string - function(x, brackets = c('[', ']'), sep=',', collapse='
', abbrev) {
if (length(brackets)==1 (nchar(brackets)1)) brackets -
unlist(strsplit(brackets, ))
terms - lapply(x, paste, collapse=sep)
terms - paste(brackets[1], terms, brackets[2], sep='')
paste(terms, collapse= ' ')
}
#' models of joint independence, of some factors wrt one or more other
factors
#' @param nf number of factors for which to generate model
#' @param table a contingency table used for factor names, typically the
output from \code{\link[base]{table}}
#' @param factors names of factors used in the model when \code{table}
is not specified
#' @param with
Re: [R] Lattice barchart with error bars
This is not an R question. Read the references. Bert Sent from my iPhone -- please excuse typos. On Jul 5, 2013, at 12:15 PM, Shaun Jackman sjack...@gmail.com wrote: Hmm. Interesting point, Bert. I don't know whether the notches show the 95% confidence interval or the median, or the 95% confidence interval that two non-overlapping notches have different medians. You're saying it's the latter? Anyone know what the 95% confidence interval of the median would be? Cheers, Shaun The notches (if requested) extend to +/-1.58 IQR/sqrt(n). This seems to be based on the same calculations as the formula with 1.57 in Chambers et al. (1983, p. 62), given in McGill et al. (1978, p. 16). They are based on asymptotic normality of the median and roughly equal sample sizes for the two medians being compared, and are said to be rather insensitive to the underlying distributions of the samples. The idea appears to be to give roughly a 95% confidence interval for the difference in two medians. On 5 July 2013 11:48, Bert Gunter gunter.ber...@gene.com wrote: Be careful! You are talking about 2 different varieties of apples here. As I read it, the CI's in the cancer data, which I know is just for example purposes, are CI's for the **individual means**; the notches in boxplots are nonparametric and for 2 groups with roughly equal sample sizes, The idea appears to be to give roughly a 95% confidence interval for the **difference** in two medians. (from ?boxplot.stats). So I'm not sure which you want, but they are certainly different (by a factor of around sqrt(2),right?), even if both are for the mean or both are for the median. Cheers, Bert On Fri, Jul 5, 2013 at 11:28 AM, David Winsemius dwinsem...@comcast.net wrote: On Jul 5, 2013, at 11:15 AM, Shaun Jackman wrote: Hi Bert, Dennis, I'll agree that using a barchart was a poor choice. I was in fact using a notched bwplot to show the median and confidence interval of the median. In this case it's the median and confidence interval that I want to highlight, and I find that the visual noise of the box and whiskers is detracting from the focus, and those wee notches are not much to focus on. So, I'd like to draw a stripplot with error bars, preferably in Lattice. Let's call this a TIE fighter plot. Any suggestions? I like the TIE fighter label. Try this: library(latticeExtra) data(USCancerRates) segplot(reorder(factor(county), rate.male) ~ LCL95.male + UCL95.male, data = subset(USCancerRates, state == Washington), draw.bands = FALSE, centers = rate.male, segments.fun = panel.arrows, ends = both, angle = 90, length = 1, unit = mm) It's what Sarkar has recommended in the past when this request has been posted. -- David Cheers, Shaun On 4 July 2013 18:00, Dennis Murphy djmu...@gmail.com wrote: If you consult the lattice package help, you'll discover there is no panel_errorbar() function, which would imply the package developers have a distaste for that type of graphic. If you fish around the R-help archives, though, you might be able to find someone who wrote a function to do error bars in lattice. (Use a searchable archive such as Nabble to hunt for it.) Error bar plots are easier to do in the ggplot2 package, since there is a specific function to generate the error bar 'geometry' (geom_errorbar). See http://docs.ggplot2.org/current/ for an expanded version of the package help pages, which include the graphs generated by the code. I believe there's also a base graphics version that you can get from the gplots package, but I don't know a lot about it. Dennis On Thu, Jul 4, 2013 at 2:53 PM, Shaun Jackman sjack...@gmail.com wrote: Hi, I'd like to draw a lattice barchart of means with error bars to show the standard deviation. I have the barchart, how do I add the error bars? require(datasets) require(lattice) x - aggregate(weight ~ Diet, ChickWeight, function(x) c(mean=mean(x), sd=sd(x))) barchart(weight[,'mean'] ~ Diet, x) Thanks, Shaun __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius Alameda, CA, USA -- 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
Re: [R] Data Package Query
Hello, When I run the below syntax: Trial-read.table(Trial.txt,header=TRUE) Trial save.image(file=Trial.RData) load(Trial.RData) fit-logistf(data=Trial, y~x1+x2) summary(fit) AIC(fit) I am getting the below error: AIC(fit) Error in UseMethod(logLik) : no applicable method for 'logLik' applied to an object of class logistf Can you please help with that? Regards, Yasmine Date: Sat, 29 Jun 2013 05:05:28 -0800 From: jrkrid...@inbox.com Subject: Re: [R] Data Package Query To: y_re...@hotmail.com; rolf.tur...@xtra.co.nz; jdnew...@dcn.davis.ca.us CC: r-help@r-project.org I don't know what you think data(Trial) is doing but what it in fact is doing is trying to load a stored data set called Trial and it does not exist. Have a look at ?data to see what I mean. In your program data(Trial) is redundant, well actually closer to meaningless. Trial is already loaded since you created it in the read statement John Kane Kingston ON Canada -Original Message- From: y_re...@hotmail.com Sent: Fri, 28 Jun 2013 12:31:11 + To: rolf.tur...@xtra.co.nz, jdnew...@dcn.davis.ca.us Subject: Re: [R] Data Package Query hello, please advice what is wrong at the below syntax: Trial-read.table(Trial.txt,header=TRUE) Trial save.image(file=Trial.RData) data(Trial) fit-logistf(data=Trial, y~x1+x2) and here is the error I get: Warning message: In data(Trial) : data set ?Trial? not found regards, yasmine Date: Fri, 28 Jun 2013 10:29:21 +1200 From: rolf.tur...@xtra.co.nz To: jdnew...@dcn.davis.ca.us CC: y_re...@hotmail.com; r-help@r-project.org Subject: Re: [R] Data Package Query On 28/06/13 04:47, Jeff Newmiller wrote: SNIP A common error by beginners (which may or may not be your problem in this case) is to create a variable called data. Unfortunately this hides the function named data and from that time forward that R session doesn't work when you type example code that uses the data function. SNIP This is simply not true. I believe it *used* to be true, sometime wa back, but hasn't been true for years. The R language is much cleverer now. If there is a function melvin() somewhere on the search path and also a data object melvin (earlier on the search path) then doing melvin(whatever) will correctly call the function melvin() with no complaints. The R language can tell by the parentheses that you mean the *function* melvin and not the data object melvin. E.g. data - 42 require(akima) akima Error: object 'akima' not found data(akima) # No error message, nor nothin'! akima # The data set akima is displayed. All that being said it is ***BAD PRACTICE***, just in terms of comprehensibility and avoiding confusion, to give a data set set the same name as a function (either built in, or one of your own). fortune(dog) is relevant. cheers, Rolf Turner [[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. FREE 3D MARINE AQUARIUM SCREENSAVER - Watch dolphins, sharks orcas on your desktop! Check it out at http://www.inbox.com/marineaquarium [[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] save rds as text
I created a table like this: Analysis of Variance Table Response: dati Df Sum Sq Mean Sq F value Pr(F) groups 2114 57.00 76 4.134e-11 *** Residuals 24 180.75 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 and saved it to a variable called k. When I tried to save the variable to txt file using dput(k, file = ..\\Readouts\\anova_spike_number.txt) or lapply(k, write, ..\\Readouts\\anova_spike_number.txt, append=TRUE) it destroys the formatting. Then I saved to rds using saveRDS(k, ..\\Readouts\\anova_spike_number.rds) and the format was preserved. I would like to convert the rds to txt preserving the table format. Or save it as excel file or pdf would also be fine. (If there is a way to avoid passing through rds file better, but still OK using rds). Thank you very much. Best regards, Alberto -- View this message in context: http://r.789695.n4.nabble.com/save-rds-as-text-tp4670925.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Non-linear modelling with several variables including a categorical variable
Off list I sent the OP a note that wrapnls() from nlmrt calls nls after nlxb finishes. This is not, of course, super-efficient, but returns the nls-structured answer. JN On 13-07-05 06:00 AM, r-help-requ...@r-project.org wrote: Message: 49 Date: Fri, 5 Jul 2013 08:30:39 +0700 From: Robbie Weteringsrobbie.weteri...@gmail.com To: Prof J C Nash (U30A)nas...@uottawa.ca,r-help@r-project.org Subject: Re: [R] Non-linear modelling with several variables including a categorical variable Message-ID: CAFe5dHZdXFbFtwKmTE1_QPi1rqNGsd+=82tproyfs6mg6zm...@mail.gmail.com Content-Type: text/plain Dear Prof. Nash, I tried to run nls with the nlxb function and as you mention it is fairly slower in terms of running the code. However, if I would have used this function earlier I would have saved a lot of time trying to find the start values. The output looks a little bit sloppy but I think it is very usefull to use in combination with nls. Thanks Robbie __ 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] geeglm
How to extract the Std.err and the alpha estimated value from the geeglm function in R. -- View this message in context: http://r.789695.n4.nabble.com/geeglm-tp4670936.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Filter Dataframe for Alarm for each column.
Hi here i have a dataframe called MyDF. a-c(1,1,1,1,1,0,0,0,1,1) b-c(1,1,0,1,1,0,0,0,1,1) c-c(1,1,1,1,1,1,1,0,1,1) d-c(1,1,1,1,1,1,1,1,0,1) MyDF-data.frame(DWATT=a,TNH=b,CSGV=c,BIX=d) My requirement is, here i need a function - to get for a particular row number(s), when particular column(s) value change from zero to one (for the first change). For example, Using MyDF, DWATT TNH CSGV BIX 1 11 1 1 11 1 1 01 1 1 11 1 1 11 1 0 01 1 0 01 1 0 00 1 1 11 0 1 11 1 Here i want to know, the row number where TNH-column and BIX-column values change happening from one-to-zero for the first time. Answer should be a dataframe with single row. So here answer should return a dataframe like this. TNH BIX -- 3 9 i used some ways to get a solution using loops. But there is a bulk files with bulk rows to process. So performace is most important. Could someone please suggest better ideas ? Thanks, Antony. -- View this message in context: http://r.789695.n4.nabble.com/Filter-Dataframe-for-Alarm-for-each-column-tp4670950.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] list construction with automatic names
I would suggest seeing if you can use the `lapply` or `sapply` function,
that way most of the details of creating the object are taken care of for
you. If you want the list named then you can use the `names` or `setNames`
function, or the USE.NAMES argument to `sapply`.
On Thu, Jul 4, 2013 at 6:43 AM, Alex van der Spek do...@xs4all.nl wrote:
I often find myself (wanting t)o constructing lists or
data.frames
Apologies; previous post should have said Read R inforno on 'Growing
Objects' and should have added the URL:
http://www.burns-stat.com/pages/Tutor/R_inferno.pdf
S Ellison
***
This email and any attachments are confidential. Any use...{{dropped:8}}
__
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.
Thank you! That was very helpful. The R-Inferno pdf makes nice reading too!
__
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.
--
Gregory (Greg) L. Snow Ph.D.
538...@gmail.com
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Transferring commas in character vector to expression
I'm trying to format a given character vector as an expression with Greek symbols to be used in labeling axis ticks. Thanks to some help from David Winsemius, I've learned how to make the substitution and place the Greek symbols in, however I've run into another problem: Some of my labels have commas in them, so when the parse command is executed, there is an unexpected symbol error. For example: x - c(aa, aaa,aa, bb*Delta*b, *Delta*cc,c) parse(text = x) Error in parse(text = x) : text:2:4: unexpected ',' 1: aa 2: aaa, ^ I've tried various iterations of wrapping the commas in interior quotes (aaa\,\aa), but then the error shifts to the quote. I see in plotmath that 'list(a,b,c)' gives me comma separated values, but I haven't been able to work out how to get this construction for elements that have a comma. Is this possible? -- Eric Archer, Ph.D. Southwest Fisheries Science Center NMFS, NOAA 8901 La Jolla Shores Drive La Jolla, CA 92037 USA 858-546-7121 (work) 858-546-7003 (FAX) Marine Mammal Genetics Group: swfsc.noaa.gov/prd-mmgenetics ETP Cetacean Assessment Program: swfsc.noaa.gov/prd-etp The universe doesn't care what you believe. The wonderful thing about science is that it doesn't ask for your faith, it just asks for your eyes. - Randall Munroe Lighthouses are more helpful than churches. - Benjamin Franklin ...but I'll take a GPS over either one. - John C. Craig George [[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] Subset and order
David Carlson dcarlson at tamu.edu writes: It may be that single and efficient are opposing goals. Two steps lets you create the subset and then just order each query. Alternatively, if the data do not change often, create an ordered version and query that. I don't know the data.table package well, but it seems as though this might be an appropriate job for it. __ 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] Operations on a big data frame
Hi, May be this helps: dat1- read.table(text= P1_prom Nom 1 -6.17 Pt_00187 2 -6.17 Pt_00187 3 -6.17 Pt_00187 4 -6.17 Pt_00187 5 -6.17 Pt_00187 6 -6.17 Pt_01418 7 -5.77 Pt_01418 8 -5.37 Pt_01418 9 -4.97 Pt_01418 10 -4.57 Pt_01418 ,sep=,header=TRUE,stringsAsFactors=FALSE) library(zoo) dat1$PT_promMean-rollmean(dat1$P1_prom,5,fill=NA,align=left) dat1 # P1_prom Nom PT_promMean #1 -6.17 Pt_00187 -6.17 #2 -6.17 Pt_00187 -6.17 #3 -6.17 Pt_00187 -6.09 #4 -6.17 Pt_00187 -5.93 #5 -6.17 Pt_00187 -5.69 #6 -6.17 Pt_01418 -5.37 #7 -5.77 Pt_01418 NA #8 -5.37 Pt_01418 NA #9 -4.97 Pt_01418 NA #10 -4.57 Pt_01418 NA A.K. Hello all, I have a big data frame that looks like this: P1_prom Nom 1 -6.17 Pt_00187 2 -6.17 Pt_00187 3 -6.17 Pt_00187 4 -6.17 Pt_00187 5 -6.17 Pt_00187 6 -6.17 Pt_01418 7 -5.77 Pt_01418 8 -5.37 Pt_01418 9 -4.97 Pt_01418 10 -4.57 Pt_01418 - - - 25000 where Nom represents a point in a map, and P1_prom represents the value of an operation we perfomed on each point (note that we performed 5 repetitions for each point, hence, each point has 5 values). What I am trying to do, with no success, is to create a new column, in which each row corresponds to the mean value of P1_prom for each point. So basically what I need the program to do is to write in the first row of the new column the average of the first five values of P1_prom, in the second row the average of the next five values, and so on. Could anybody guide me on how to do this. Thank you very much, Veronica __ 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] IF function
Hi, May be this helps. dat1- read.table(text= Col1,Col2 High value,9 Low value,0 High value,7 Low value,0 Low value,0 No data,0 High value,8 No data,0 ,sep=,,header=TRUE,stringsAsFactors=FALSE) dat1$Col2[dat1$Col1==No data]- NA dat1 # Col1 Col2 #1 High value 9 #2 Low value 0 #3 High value 7 #4 Low value 0 #5 Low value 0 #6 No data NA #7 High value 8 #8 No data NA A.K. Hello, I am an R novice so excuse me if this is woefully straight forward, but I have tried the help files to no avail. I am trying to identify cells in 1 column with the value of 'No data', so I can change the values in the next column to 'Null'. Currently I am struggling with the data set, as it assigns both 'No data' and 'Low values' as zero which skews my analysis. I've tried a number of different attempts but just get the error unexpected symbol ? __ 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 BootCV for coxph in pec after feature selection with glmnet (lasso)
Hi, I am attempting to evaluate the prediction error of a coxph model that was built after feature selection with glmnet. In the preprocessing stage I used na.omit (dataset) to remove NAs. I reconstructed all my factor variables into binary variables with dummies (using model.matrix) I then used glmnet lasso to fit a cox model and select the best performing features. Then I fit a coxph model using only these feature. When I try to evaluate the model using pec and a bootstrap I get an error that the prediction matrix has wrong dimensions. Suddenly the cox object has 318 variables instead of 356 variables in the dataset. I don't know why this is happening. The cox object I assign to pec and the dataframe are both of the same size. However, once pec refits the model its size changes (356 - 318 variables). Apparently something is happening during the bootstrap sampling that removes some variables. As mentioned, I used na.omit in the preprocessing so there should not be any NAs. Here are some details from my workspace: Reformat_dataSet: 368 obs. Of 356 variables print (glmnet.cox) 354 df, p=1 n= 368, number of events= 288 (354 df = 354 variables + time and status = 356 variables) Here is the pec function and the error: pec.f - as.formula(Hist(time,status) ~ 1) brierGlmnet - pec(list(glmnet.cox), data = reformat_Dataset, splitMethod=BootCV, B=50, formula = pec.f) Error in predictSurvProb.coxph(object = structure(list(coefficients = structure(c(-4.27787223119601, : Prediction matrix has wrong dimensions: 368 rows and 318 columns. But requested are predicted probabilities for 118 subjects (rows) in newdata and 356 time points (columns) This may happen when some covariate values are missing in newdata!? Here are the relevant sections of the code: trainSet - na.omit (dataset) #creat Y (survival matrix) for glmnet surv_obj - Surv(trainSet$time,trainSet$status) ## tranform categorical variables into binary variables with dummy for trainSet predict_matrix - model.matrix(~ ., data=trainSet, contrasts.arg = lapply (trainSet[,sapply(trainSet, is.factor)], contrasts, contrasts=FALSE)) ## remove the statu/time variables from the predictor matrix (x) for glmnet predict_matrix - subset (predict_matrix, select=c(-time,-status)) ## create a glmnet cox object using lasso regularization glmnet.obj - glmnet (predict_matrix, surv_obj, family=cox) # find lambda for which dev.ratio is max max.dev.index - which.max(glmnet.obj$dev.ratio) optimal.lambda - glmnet.obj$lambda[max.dev.index] # take beta for optimal lambda optimal.beta - glmnet.obj$beta[,max.dev.index] # find non zero beta coef nonzero.coef - abs(optimal.beta)0 selectedBeta - optimal.beta[nonzero.coef] # take only covariates for which beta is not zero selectedVar - predict_matrix[,nonzero.coef] # create a dataframe for trainSet with time, status and selected variables in binary representation for evaluation in pec reformat_dataSet - as.data.frame(cbind(surv_obj,selectedVar)) # create coxph object with pre-defined coefficients glmnet.cox - coxph(surv_obj ~ selectedVar, init=selectedBeta,iter=0) ## Brier score for the cox-glmnet model brierGlmnet - pec(list(glmnet.cox), data = reformat_Dataset, splitMethod=BootCV, B=50, formula = pec.f) 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.
