Re: [R] gplot heatmaps: clustering according to rowsidecolors + key.xtickfun
Oops, I forgot to mention that an bug was preventing RowSideColors from
working properly. It is fixed in version 2.14.2 of gplots which I've just
uploaded to CRAN and am attaching to this email.
-Greg
On Wed, Sep 17, 2014 at 5:12 PM, Gregory R. Warnes g...@warnes.net wrote:
Hello Tim,
Sorry about the slow response, I just found this message.
On Sep 4, 2014, at 9:23 AM, Tim Richter-Heitmann trich...@uni-bremen.de
wrote:
Hi there,
I have two questions about heatmap.2 in gplot.
My input is a simple square matrix with numeric values between 75 and 100
(it is a similarity matrix based on bacterial DNA sequences).
1. I can sort my input matrix into categories with rowsidecolors (in this
case, very conveniently by bacterial taxa). I do a clustering and
reordering of my matrix by Rowv=TRUE (and Colv=Rowv).
The question is now, can i combine the two features that the
clustering/reordering is done only for submatrices defined by the vectors
given in rowsidecolors (so, in this case, that the original ordering by
bacterial taxa is preserved)?
That would be very amazing.
Hmm.To get the individual species clustered within taxa would require
doing the hierarchical clustering first separately, then combining the
dendrograms. This should do the trick:
set.seed(1234567)
## Dummy Distances
x - matrix( rnorm(400, mean=87.5, sd=12.5/4), ncol=20)
## Dummy Taxa
taxa - sample(letters[1:4], 20, replace=T)
taxa - as.factor(taxa)
# sort the data by taxa
ord - order(taxa)
x - x[ord, ord]
taxa - taxa[ord]
rownames(x) - 1:nrow(x)
# stats:::merge.dendrogram is broken. This is the corrected version.
# See R BUG 15648
# (https://bugs.r-project.org/bugzilla3/show_bug.cgi?id=15648) for
# details
merge.dendrogram - function(x, y, ..., height) {
stopifnot(inherits(x,dendrogram), inherits(y,dendrogram))
### FIX
inx.add - function(inx, add) {
if(is.leaf(inx)) {
inx - inx + add
}
return(inx)
}
y - dendrapply(y, inx.add, add=max(unlist(x)))
### FIX
r - list(x,y)
if(length(xtr - list(...))) {
if(!all(is.d - vapply(xtr, inherits, NA, what=dendrogram))) {
xpr - substitute(c(...))
nms - sapply(xpr[-1][!is.d], deparse, nlines = 1L)
## do not simplify: xgettext needs this form
msg - ngettext(length(nms),
extra argument %s is not of class \%s\,
extra arguments %s are not of class \%s\s)
stop(sprintf(msg, paste(nms, collapse=, ), dendrogram),
domain = NA)
}
## GRW
for(i in 1:length(xtr))
{
add - max(c(unlist(r), unlist(xtr)))
print(add)
xtr[[i]] - dendrapply(xtr[[i]], inx.add, add=add)
}
## /GRW
r - c(r, xtr)
}
attr(r, members) - sum(vapply(r, attr, 0L, which=members))
h.max - max(vapply(r, attr, 0., which=height))
if(missing(height) || is.null(height))
height - 1.1 * h.max
else if(height h.max) {
msg - gettextf('height' must be at least %g, the maximal height of
its components, h.max)
stop(msg, domain = NA)
}
attr(r, height) - height
class(r) - dendrogram
midcache.dendrogram(r, quiet=TRUE)
}
## Compute dendrograms within each taxum, then merge into a combined
dendrogram
dendList - list()
for( taxon in levels(taxa) )
{
items - which(taxon==taxa)
submatrix - x[ items, items]
dend - as.dendrogram(hclust(dist(submatrix)))
dendList[[taxon]] - dend
}
names(dendList) - NULL
dends - do.call(merge.dendrogram, dendList)
## Now generate the heatmap
heatmap.2(x,
Rowv=dends,
Colv=dends,
symm=TRUE,
RowSideColors=c(red,blue,green,black)[as.numeric(taxa)],
ColSideColors=c(red,blue,green,black)[as.numeric(taxa)],
trace=none
)
2. I have set my own coloring rules with:
mypal - c(grey,blue, green,yellow,orange,red)
col_breaks = c(seq(0,74.9), seq(75.0,78.4), seq(78.5,81.9),
seq(82.0,86.4), seq(86.5, 94.5), seq(94.5,100.0))
Is it possible to pass this sequential ordering to key.xtickfun? May i ask
for an example code?
Use the ‘breaks’ and ‘col’ arguements:
## Custom color key
mypal - c(grey,blue, green,yellow,orange,red)
col_breaks - c(0,75.0,78.5,82.0,86.5,94.5,100.0)
heatmap.2(x,
Rowv=dends,
Colv=dends,
symm=TRUE,
RowSideColors=c(red,blue,green,black)[as.numeric(taxa)],
ColSideColors=c(red,blue,green,black)[as.numeric(taxa)],
trace=none,
breaks=col_breaks,
col=mypal
)
-Greg
--
Whereas true religion and good morals are the only solid foundations of
public liberty and happiness . . . it is hereby earnestly recommended to
the several States to take the most effectual measures for the
encouragement thereof. Continental Congress, 1778
Re: [R] issue with ... in write.fwf in gdata
Hi Jan, The issue isn't that the ... arguments aren't passed on. Rather, the problem is that in the current implementation the ... arguments are passed to format(), which doesn't understand the eol argument. The solution is to modify write.fwf() to explicitly accept all of the appropriate the arguments for write.table() and to only pass the ... arguments to format() and format.info(). I've just modified gdata to make this change, and have submitted the new version to CRAN as gdata version 2.8.1. -Greg On Fri, Nov 12, 2010 at 7:08 AM, Jan Wijffels janwijff...@hotmail.comwrote: Dear R-list This is just message to inform that the there is an issue with write.fwf in the gdata library (from version 2.5.0 on). It does not seem to accept further arguments to write.table like eol as the help file indicates as it stops when executing tmp - lapply(x, format.info, ...). Great package though - I use it a lot except for this function :) See example below. require(gdata) saveto - tempfile(pattern = test.txt, tmpdir = tempdir()) write.fwf(x = data.frame(a=1:length(LETTERS), b=LETTERS), file=saveto, eol=\r\n) Error in FUN(X[[1L]], ...) : unused argument(s) (eol = \r\n) sessionInfo() R version 2.12.0 (2010-10-15) Platform: x86_64-pc-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] gdata_2.8.0 loaded via a namespace (and not attached): [1] gtools_2.6.2 kind regards, Jan [[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.
[R] Correction: RStat version 2.7.*1* now available
Hello Everyone, I mis-typed the version number in the annoucement that just went out. Here is the corrected message: On Fri, Jul 18, 2008 at 6:50 AM, Gregory R. Warnes [EMAIL PROTECTED] wrote: Random Technologies LLC is pleased to announce immediate availability of RStat version 2.7.1. http://2.7.2. This version of RStat provides the features of the latest version of R from the R-project with with enterprise-level validation, documentation, software support, and consulting services from Random Technologies LLC. RStat 2.7.1 is available in Personal, Professional, and Enterprise editions. Visit the RStat(R) product page http://random-technologies-llc.com/products/products/RStat/rstat for features and pricing. Gregory R. Warnes, Ph.D. Chief Scientist Random Technologies LLC http://random-technologies-llc.com email: [EMAIL PROTECTED] tel:(585) 419-6853 fax:(585) 672-5085 [[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] Messge ``no visible binding''.
Simply add a definition for this variable to the source for your package,
e.g.;
.Last.make.date - NULL
And the warning should go away.
-g
On 6/30/08 11:44PM , Rolf Turner [EMAIL PROTECTED] wrote:
I have written a function make.fun(), which I keep in my personal
``miscellaneous'' package. It searches the current working directory
for files named ``*.R'' and sources them if their modification date
is later than the date stored in ``.Last.make.date''. (The variable
.Last.make.date then gets updated.) The function checks whether
.Last.make.date exists before comparing its value with the modification
dates of the *.R files.
When I do R CMD check on my ``miscellaneous'' package, I get a ``NOTE''
(not even a warning!) to the effect:
``no visible binding for global variable '.Last.make.date'
Is this something that I should worry about at all? Am I running any
risk by this procedure? Is there a better way to check whether code file
have modified sufficiently recently that they need to be re-sourced?
Thanks for any tips. (But no asparagus, please! :-) )
cheers,
Rolf Turner
##
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
__
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 R. Warnes, Ph.D
Program Director
Center for Computational Arts, Sciences, and Engineering
University of Rochester
Tel: 585-273-2794
Fax: 585-276-2097
Email: [EMAIL PROTECTED]
__
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] heatmap and continuous variable
Hello Hans, Do you need *both* a categorical *and* a continuous variable, or just a continuous variable? -Greg On 6/24/08 4:51PM , Hans-Ulrich Klein [EMAIL PROTECTED] wrote: Dear All, I want to plot a heat map with annotated columns. Both functions heatmap (stats) and heatmap.2 (gplots) can plot a horizontal side bar that can be used to visualize a categorical variable. In addition to a categorical variable, I would like to visualize a continuous variable. This could be done by small bars, a curve or simply numbers above the columns. (The Sample names are already located at the bottom.) Does someone have example code for this or a similar problem? Is this so special that I have to get down to the lattice package? Thank you, Hans-Ulrich -- Gregory R. Warnes, Ph.D Program Director Center for Computational Arts, Sciences, and Engineering University of Rochester Tel: 585-273-2794 Fax: 585-276-2097 Email: [EMAIL PROTECTED] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Survey: Commercial R companies
Hello Everyone, In preparation for an upcoming talk, I would like to assemble a list of companies that provide consulting, services, products, or training for R. I am already aware of a number of such companies including (in alphabetical order): BlueReference http://inference.us Insightful http://www.insightful.com Metrum Institute http://www.metruminstitute.org Metrum Research Group, LLC http://www.metrumrg.com Random Technologies LLC http://random-technologies-llc.com REvolution Computing Inc http://revolution-computing.com Statistics.com http://www.statistics.com XLSolutions Corporation http://www.xlsolutions-corp.com If you are aware of a company that provides consulting, services, products, or training for R, please complete the form below and forward it directly to me at [EMAIL PROTECTED] by Wednesday June 11. I will collate all of the responses and post a complete list back to R-help. Thanks, -Greg Gregory R. Warnes, Ph.D. Associate Professor Center for Biodefence Immune Modeling and Department of Biostatistics and Computational Biology University of Rochester === Cut Here, Send to [EMAIL PROTECTED] == Company Name: Address (including country): Web Site: Email Address: Phone Number: Type of Products and Services: [ ] Commercial versions of R Product Name: [ ] General Training (e.g. Introduction to R, Programming in R, etc.) Courses Taught: [ ] Domain-specific Training (e.g. R for PK/PD Modeling, Bayesian Modeling using R, etc.) Domain: [ ] Add-on packages (Off the shelf, not custom) Package Names: [ ] Software/package development (custom) [ ] Parallel Computing Tools [ ] Statistical Consulting utilizing R Area of expertise: [ ] GUI environments Product Name: [ ] Qualification/Validation Services [ ] Other: [[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] SQL INSERT using RMySQL
Hi All,
I figured out my problem. There was a combination of lack of
understanding on my part, and a bit of missing functionality. I made
a small patch to the rmysqlWriteTable() function passes the field
names to MySQL corresponding to the data columns passed in:
diff -ru RMySQL.orig/R/MySQLSupport.R RMySQL/R/MySQLSupport.R
--- RMySQL.orig/R/MySQLSupport.R2007-05-31 22:36:02.0 -0400
+++ RMySQL/R/MySQLSupport.R 2008-04-11 17:50:29.0 -0400
@@ -616,7 +616,9 @@
on.exit(unlink(fn), add = TRUE)
sql4 - paste(LOAD DATA LOCAL INFILE ', fn, ',
INTO TABLE , name,
- LINES TERMINATED BY '\n' , sep=)
+ LINES TERMINATED BY '\n' ,
+ ( , paste(names(field.types), collapse=, ),
);,
+ sep=)
rs - try(dbSendQuery(new.con, sql4))
if(inherits(rs, ErrorClass)){
warning(could not load data into table)
I also defined a useful function for describing the structure of an
existing table:
setGeneric(
dbDescribeTable,
function(conn, name, ...)
standardGeneric(dbDescribeTable),
valueClass = character
)
setMethod(
dbDescribeTable,
signature(conn=MySQLConnection, name=character),
def = function(conn, name, ...){
rs - dbGetQuery(conn, paste(describe, name))
fields - rs$Type
names(fields) - rs$Field
if(length(fields)==0)
fields - character()
fields
},
valueClass = character
)
And I now have working code:
## Columns in the table
dbDescribeTable(con, past_purchases)
idcustomer_id item_upc
int(10) unsigned int(11) bigint(12)
suggested quantity total
tinyint(1) int(11) int(11)
on_sale actual_price featured
tinyint(1) double tinyint(1)
date
date
## columns in my data (note the absence of the primary key 'id')
head(fulldata)
customer_iditem_upc suggested quantity total on_sale
1 3 632 FALSE1 1 FALSE
2 3 733 FALSE1 1 FALSE
3 3 1116095 FALSE1 1 FALSE
4 3 1117164 FALSE1 1 FALSE
5 3 1117210 FALSE1 1 FALSE
6 3 1119092 FALSE1 1 FALSE
actual_price featured date
110.49FALSE 2008-03-22
2 4.99FALSE 2008-03-22
3 5.49FALSE 2008-03-22
4 9.99FALSE 2008-03-22
5 4.19FALSE 2008-03-22
6 3.99FALSE 2008-03-22
dim(fulldata)
[1] 75 9
## Size of the table before adding my data
dbGetQuery(con, SELECT COUNT(ID) FROM past_purchases)[1,1]
[1] 675
## Insert the data
dbWriteTable(
+ con,
+ past_purchases,
+ value=fulldata,
+ overwrite=FALSE,
+ append=TRUE,
+ row.names=FALSE #,
+ #field.types=field.types
+ )
[1] TRUE
## Size of the table after adding my data
dbGetQuery(con, SELECT COUNT(ID) FROM past_purchases)[1,1]
[1] 750
-Greg
On Apr 11, 2008, at 10:57PM , Chris Stubben wrote:
Greg,
If you have a MySQL table with an auto_increment field, you could just
insert a NULL value into that column and the database will
increment the key
(it may not work in SQL STRICT mode, I'm not sure). I don't think
there's
any way to specify which columns you want to load data into using
dbWriteTable yet, but that would be a nice feature since LOAD data now
allows that (and SET syntax and other options).
Try this code below - I used cbind(NA, x) to insert a null into the
first
column.
Chris
dbSendQuery(con, create table tmp (id int not null auto_increment
primary
key, a char(1), b char(1)))
MySQLResult:(369,1,67)
x-data.frame( a=letters[1:3], b=letters[4:6])
x
a b
1 a d
2 b e
3 c f
dbWriteTable(con, tmp, cbind(NA, x), row.name=FALSE, append=TRUE)
[1] TRUE
dbWriteTable(con, tmp, cbind(NA, x), row.name=FALSE, append=TRUE)
[1] TRUE
dbReadTable(con, tmp)
id a b
1 1 a d
2 2 b e
3 3 c f
4 4 a d
5 5 b e
6 6 c f
Gregory. R. Warnes wrote:
Hi All,
I've finally gotten around to database access using R. I'm happily
extracting rows from a MySQL database using RMySQL, but am having
problems appending rows to an existing table.
What I *want* to do is to append rows to the table, allowing the
database to automatically generate primary key values. I've only
managed to add rows by using
dbWriteTable( con, past_purchases, newRecords, overwrite=FALSE,
append=TRUE, ...)
And this only appears to properly append rows (as opposed to
overwriting them) IFF
Re: [R] prediction intervals from a mixed-effects models?
On Apr 13, 2008, at 1:41PM , Dieter Menne wrote: Spencer Graves spencer.graves at pdf.com writes: How can I get prediction intervals from a mixed-effects model? Consider the following example: library(nlme) fm3 - lme(distance ~ age*Sex, data = Orthodont, random = ~ 1) df3.1 - with(Orthodont, data.frame(age=seq(5, 20, 5), Subject=rep(Subject[1], 4), Sex=rep(Sex[1], 4))) predict(fm3, df3.1, interval='prediction') # M01 M01 M01 M01 # 22.69012 26.61199 30.53387 34.45574 # NOTE: The 'interval' argument to the 'predict' function was ignored. # It works works for an 'lm' object, but not an 'lme' object. In theory, ci from gmodels should work library(gmodels) ci(df3.1) However, I get an error message. I will forward this to Greg to let him check if I did something stupid. gmodels::ci() will give confidence intervals for the fixed effects via ci(fm3) ci() won't generate prediction intervals for individual parameters, and according to ?predict.lme it won't either. -Greg Dieter __ 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] lme and confidence intervals
FWIW, the ci() function in the gmodels package supports generating confidence intervals for the fixed effects of lme objects. -G On Apr 8, 2008, at 1:34PM , Dieter Menne wrote: Cristian Carranza cristiancarranza_1 at hotmail.com writes: After fitting a mixed effects model to repeated measurements data set, and after several unsuccessful atempts to make a simple plot of the confidence interval for the fitted model, I gave up and now I am asking for help in this useful list. Could anyone be so kind to give me some code lines in order to make a plot of the fitted equation and the correspondent confidence interval? The fitted equation is simple: use predict in its variations. For the confidence interval, most people on the list will ask back what it means. But I know my colleagues, they or their reviewers insist on having error bar in the graphics, and don't care about the interpretation. Dimitri has given an approach that could be used to produce some reviewer-satisfaction plots limited borborygmy. http://finzi.psych.upenn.edu/R/Rhelp02a/archive/118174.html Note that his method is meant to estimate well-defined confidence intervals for the coefficients, but you can use the profiled fits and compute the confidence intervals at the data points. Alternatively, you could use mcmc to get a range of curves for averaging. Dieter __ 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] Reading microsoft .xls format and openoffice OpenDocument files
Hello Ajay, I'm the author of the gdata package. If you send me a copy of an .XLS file that doesn't work with read.xls(), I'll see about fixing the code. -Greg On Mar 7, 2008, at 6:17AM , Ajay Shah wrote: 1. I have used gdata::read.xls() with much happiness. But every now and then it breaks. I have not, as yet, been able to construct a mental model about the class of .xls files for which it works. Does someone have a simple rule for predicting the circumstances under which it will work? 2. Just like there is a read.xls(), it'd be great if we have a read.ods() which directly reads files from openoffice. This should be easier than grokking Microsoft formats given that openoffice is gpl. I hunted a bit and couldn't find any. Does someone know how we might approach this? Am I correct in thinking that our goal is reading OpenDocument files (http://en.wikipedia.org/wiki/OpenDocument) ? -- Ajay Shah http://www.mayin.org/ ajayshah [EMAIL PROTECTED] http:// ajayshahblog.blogspot.com *(:-? - wizard who doesn't know the answer. __ 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] testing for significantly different slopes
What's the problem? The problem is that I would like to do a pair-wise comparison between the multiple slopes. For example with this model: lm1 - lm (Sepal.Length ~ Species/Sepal.Width -1, data=iris) # truncated output from summary(lm1) # just the slope terms Speciessetosa:Sepal.Width 0.6905 0.1657 4.166 5.31e-05 *** Speciesversicolor:Sepal.Width 0.8651 0.2002 4.321 2.88e-05 *** Speciesvirginica:Sepal.Width0.9015 0.1948 4.628 8.16e-06 *** If I wanted to test the hypothesis that Speciessetosa:Sepal.Width was different than Speciesversicolor:Sepal.Width, what course of action should I take? I have found an example in the gmodels package, using the estimable() function, but the documentation is not clear to me. The gmodels::estimable function will allow you to test any hypothesis constructed from the model coefficients. Each row of the contrast matrix 'cm' is applied to the vector of coefficients 'beta' individually and tested (by default) against the null: cm[i] %*% beta == 0 This is accomplished using t-tests using the appropriate transformation of the model parameter covariance matrix to determine the correct standard deviation. Here is the example: library(gmodels) # example from manual lm1 - lm (Sepal.Length ~ Sepal.Width + Species + Sepal.Width:Species, data=iris) cm - rbind( 'Setosa vs. Versicolor' = c(0, 0, 1, 0, 1, 0), 'Setosa vs. Virginica'= c(0, 0, 0, 1, 0, 1), 'Versicolor vs. Virginica'= c(0, 0, 1,-1, 1,-1) ) estimable(lm1, cm) This *appears* to test what I am after, but I am not clear on how the 'cm' argument works. This example code tests whether the intercept *and* slope are equal across species, with a little extra complexity because the species=Setosa is represeted as the 'baseline' group that is included in the intercept term. To test whether just the slope Speciesversicolor:Sepal.Width is equal to the slope Speciesvirginica:Sepal.Width, you can do a single contrast row: estimable(lm1, c(0, 0, 0, 0, 1, -1) ) Estimate Std. Errort value DF Pr(|t|) (0 0 0 0 1 -1) -0.03645676 0.2793277 -0.1305161 144 0.8963403 For the one-contrast case, one can use the shortcut: estimable(lm1, c(Speciesversicolor:Sepal.Width=1, Speciesvirginica:Sepal.Width=-1) ) Estimate Std. Errort value DF Pr(|t|) (0 0 0 0 1 -1) -0.03645676 0.2793277 -0.1305161 144 0.8963403 Now, to each of the pairwise tests in a single call, construct a matrix with one row per contrast, and one column per model parameter: cm - rbind( + 'Setosa vs. Versicolor Slope' = c(0, 0, 0, 1, -1, 0), + 'Setosa vs. Virginica Slope'= c(0, 0, 0, 1, 0, -1), + 'Versicolor vs. Virginica Slope'= c(0, 0, 0, 0, 1, -1) + ) colnames(cm) - names(coef(lm1)) #print(cm) # omitted here for brevity, but instructive estimable(lm1, cm) Estimate Std. Errort value DF Pr(|t|) Setosa vs. Versicolor Slope-0.17458800 0.2598919 -0.6717717 144 0.5028054 Setosa vs. Virginica Slope -0.21104476 0.2557557 -0.8251811 144 0.4106337 Versicolor vs. Virginica Slope -0.03645676 0.2793277 -0.1305161 144 0.8963403 I hope this helps. -Greg __ 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] compress data on read, decompress on write
You might look at storing the data using R's raw data type... -G On Feb 28, 2008, at 5:38PM , Ramon Diaz-Uriarte wrote: Dear Christos, Thanks for your reply. Actually, I should have been more careful with language: its not really a sparse matrix, but rather a ragged array that results from a more compact representation we though of for the hidden states in a Hidden Markov Model in many runs of MCMC. However, it might make sense for us to check sparseMatrix and see how its done there. Thanks, R On Thu, Feb 28, 2008 at 7:49 PM, Christos Hatzis [EMAIL PROTECTED] wrote: Ramon, If you are looking for a solution to your specific application (as opposed to a general compression/ decompression mechanism), it might be worth checking out the Matrix package, which has facilities for storing and manipulating sparse matrices. The sparseMatrix class stores matrices in the triplet representation (i.e. only indices and values of the non-zero elements) and this affords great compression ratios, depending on the size and degree of sparseness of the matrix. -Christos -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Ramon Diaz- Uriarte Sent: Thursday, February 28, 2008 1:18 PM To: [EMAIL PROTECTED] Subject: [R] compress data on read, decompress on write Dear All, I'd like to be able to have R store (in a list component) a compressed data set, and then write it out uncompressed. gzcon and gzfile work in exactly the opposite direction. What would be a good way to handle this? Details: -- We have a package that uses C; part of the C output is a large sparse matrix. This is never manipulated directly by R, but always by the C code. However, we need to store that data somewhere (inside an R object) for further calls to the functions in our package. We'd like to store that matrix as part of the R object (say, as an element of a list). Ideally, it would be stored in as compressed a way as possible. Then, when we need to use that information, it would be decompressed and passed to the C function. I guess one way to do it is to have C deal with the compression and uncompression (e.g., using zlib or the bzip2 libraries) and then use readBin, etc, from R. But, if I can, I'd like to avoid our C code having to call zlib, etc, so as to make our package easily portable. Thanks, R. -- Ramon Diaz-Uriarte Statistical Computing Team Structural Biology and Biocomputing Programme Spanish National Cancer Centre (CNIO) http://ligarto.org/rdiaz __ 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. -- Ramon Diaz-Uriarte Statistical Computing Team Structural Biology and Biocomputing Programme Spanish National Cancer Centre (CNIO) http://ligarto.org/rdiaz __ 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] Loading Data to R
On Microsoft Windows systems, it may be more convenient to install and use the XLSReadWRite packge. For non-windows systems, the gdata package provides this function, but requires perl to be present. -Greg (Maintainer of gdata) On Feb 9, 2008, at 1:09PM , Henrique Dallazuanna wrote: You need library(gdata) before On 08/02/2008, Wensui Liu [EMAIL PROTECTED] wrote: # READ DATA FROM XLS FILE # xls - read.xls(file = C:/projects/Rintro/Part01/export.xls, sheet = 3, type = data.frame, from = 1, colNames = TRUE) On Feb 8, 2008 3:49 PM, Christine Lynn [EMAIL PROTECTED] wrote: This is the most basic question ever...I haven't used R in a couple years since college so I forget and haven't been able to find what I'm looking for in any of the manuals. I just need to figure out how to load a dataset into the program from excel! Thanks! CL [[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. -- === WenSui Liu ChoicePoint Precision Marketing Phone: 678-893-9457 Email : [EMAIL PROTECTED] Blog : statcompute.spaces.live.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. -- Henrique Dallazuanna Curitiba-Paraná-Brasil 25° 25' 40 S 49° 16' 22 O __ 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] Boxplots illustrating the fixed effects in a lme object
As with many things, I suspect that someone has created plots of this
sort, but no-one seems to have contributed a function to do so. If
you find/create one, send it to me and I'll be glad to include it in
either gmodes or gplots as appropriate.
One place to check before writing one is Frank Harrell's Hmisc or
design packages, which do have functions for plotting estimates and
standard errors for (at least) standard glm's. These might work for
lme's as well...
-G
On Nov 30, 2007, at 10:09AM , CG Pettersson wrote:
Hello all,
I posted a similar question recently but I suspect that it was not
well
enough formulated to trigger any answers. So I try again:
Is there a way to produce boxplots (or something similar) that uses
the
estimated fixed effects of an lme{nlme} object? When I want to know
the
mean, I use estimable{gmodels} together with an appropriate matrix.
I did
look for a similar tool to estimable in gplots but did not find any.
As nearly everything is possible to do in R, the reason I did not
get any
respons ought to be that the question is too trivial. So, please
give me a
hint of how to proceed.
Cheers
/CG
--
CG Pettersson, Ph.D.
Swedish University of Agricultural Sciences (SLU)
Dept. of Crop Production Ecology. Box 7043.
SE-750 07 Uppsala, Sweden
[EMAIL PROTECTED]
__
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] Heatmap problem
Hi Allen, The dotted line is an 'extra' feature to help quantify the difference between adjacent color boxes. It can be turned off by using the argument heatmap.2(..., trace=none). For more details do library(gplots) ?heatmap.2 -G On Nov 22, 2007, at 12:03PM , affy snp wrote: Hi Marco, Thanks! I actually tried heatmap.2 earlier but I observed dot line along each column generated as well. Allen On 11/22/07, Boks, M.P.M. [EMAIL PROTECTED] wrote: Try Heatmap.2 in the gplots package. BW, Marco -Oorspronkelijk bericht- Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] project.org] Namens Jim Lemon Verzonden: donderdag 22 november 2007 11:21 Aan: affy snp CC: r-help@r-project.org Onderwerp: Re: [R] Heatmap problem affy snp wrote: Hi friends, I used heatmap(as.matrix(y2),col=rainbow(256),scale = column) to generate the heatmap. But it did not show the code that which color correspond the value. Is there any parameter for this in heatmap()? Hi Allen, If you want colors corresponding to the values with a color legend, have a look at color2D.matplot in the plotrix package. Jim __ 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-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] sequence of vectors
Well, this is a natural thing to program up using 3 nested 'for',
loops. Alternatively, one could use something like:
combn - function( ..., l=list(...) )
+ {
+ lens - sapply( args, length)
+ ncomb - prod(lens)
+ retval - matrix(ncol=length(args), nrow=ncomb)
+ for(i in 1:length(args))
+ {
+ retval[,i] - rep(args[[i]], length=ncomb)
+ }
+ retval
+ }
combn(a=1:3, b=4:5, c=6:8)
[,1] [,2] [,3]
[1,]146
[2,]257
[3,]348
[4,]156
[5,]247
[6,]358
[7,]146
[8,]257
[9,]348
[10,]156
[11,]247
[12,]358
[13,]146
[14,]257
[15,]348
[16,]156
[17,]247
[18,]358
On Nov 19, 2007, at 9:59AM , Marc Bernard wrote:
Dear All,
I wonder if there is any R function to generate a sequence of
vectors from existing ones. For example:
x 1- c(1,2,3)
x2 - c(4,5)
x3 - c(6,7,8)
The desired output is a list of all 3*2*3 = 18 possible
combinations of elements of x1,x2 and x3. One element for example
is (1,4,6).
Many thanks in advance,
Bernard
-
[[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] R as a programming language
(2) More process and I/O facilities, specifically I'd like
forking and
something like a functionconnection which works like a
textconnection but obtains input from / feeds output to a
function.
This would allow running an external process that receives input
*and* provides output to R -- currently, pipes allow either one
or the other, but not both.
Both of these features are provided by the 'fork' package, at least
on *nix systems. From the example section of fork::fork,
# Use a socket to communicate between two processes. Information
# typed on the console, which is read by the initial process,
will be send
# to the child process for display.
###
## Not run:
send - function()
{
pid - getpid()
con1 - socketConnection(Sys.info()[nodename], port = 6011)
i - 1
while(TRUE)
{
cat([,pid,] ,i,: ,sep=)
data - readLines(stdin(), n=1)
writeLines(data, con=con1)
if( length(grep(quit, data))0 )
break;
i - i+1
}
close(con1)
}
recieve - function()
{
pid - getpid()
con2 - socketConnection(port = 6011, block=TRUE,
server=TRUE)
i - 1
while(TRUE)
{
data - readLines(con2, n=1)
cat([,pid,] ,i,: ,sep=)
writeLines(data, stdout())
if( length(grep(quit, data))0 )
break;
i - i+1
}
close(con2)
}
## Not run:
# Important: if we aren't going to explicitly wait() for the child
# process to exit, we need to tell the OS that we are ignoring
child
# process signals so it will let us die cleanly.
signal(SIGCHLD,ignore)
## End(Not run)
# fork off the process
pid - fork(recieve)
# start sending...
send()
## End(Not run)
-Greg
With these two, I'd be quite inclined to use R as my primary scripting
language. (I tried to get there a while ago and wrote a filterpipe
patch to achieve (2), currently I work use a mix of Python and R).
Best regards, Jan
--
+- Jan T. Kim
---+
| email:
[EMAIL PROTECTED] |
| WWW: http://www.cmp.uea.ac.uk/people/
jtk |
*-= hierarchical systems are for files, not for humans
=-*
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-
guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to test combined effects?
Ooops. One typo in the estimable command:
estimable(ModelFit, c('IQ:age'=1, 'IQ:I(age^2)'= 1, 'IQ:I(age^3)' =
1))
(Remove a trailing space in the second string.)
-G
On Nov 2, 2007, at 11:51AM , Gregory Warnes wrote:
Hello Gang,
First, if you would like to performa an overall test of whether the
IQ interactions are necessary, you may find it most useful to use
anova to compare a full and reduced model. Something like:
ModelFit.full -lme(mct~ IQ*age+IQ*I(age^2)+IQ*I(age^3), MyData,
random=~1|ID)
ModelFit.reduced -lme(mct~ IQ + age+I(age^2)+I(age^3), MyData,
random=~1|ID)
anova(ModelFit.full, ModelFit.reduced, test=F)
Second, you don't have the syntax right for estimable(). As
described and shown by example in the manual page. The correct
syntax is:
library(gmodels)
estimable(ModelFit, c('IQ:age'=1, 'IQ:I(age^2) '= 1, 'IQ:I(age^3)' =
1))
Note the pattern of quoted name, followed by '=', and then the value
1 (not zero). This will perform a single joint test whether these
three coefficients are zero.
-G
On Oct 30, 2007, at 5:26PM , Gang Chen wrote:
Dieter,
Thank you very much for the help!
I tried both glht() in multcomp and estimable() in gmodels, but
couldn't get them work as shown below. Basically I have trouble
specifying those continuous variables. Any suggestions?
Also it seems both glht() and estimable() would give multiple t
tests. Is there a way to obtain sort of partial F test?
ModelFit-lme(mct~ IQ*age+IQ*I(age^2)+IQ*I(age^3), MyData,
random=~1|ID)
anova(ModelFit)
mDF denDF F-value p-value
(Intercept) 1 257 54393.04 .0001
IQ 1 215 3.02 0.0839
age 1 25746.06 .0001
I(age^2)1 257 8.80 0.0033
I(age^3)1 25721.30 .0001
IQ:age 1 257 1.18 0.2776
IQ:I(age^2) 1 257 0.50 0.4798
IQ:I(age^3) 1 257 0.23 0.6284
library(multcomp)
glht(ModelFit, linfct = c(IQ:age = 0, IQ:I(age^2) = 0, IQ:I
(age^3) = 0))
Error in coefs(ex[[3]]) :
cannot interpret expression 'I''age^2' as linear function
library(gmodels)
estimable(ModelFit, rbind('IQ:age'=0, 'IQ:I(age^2) = 0', 'IQ:I
(age^3) = 0'))
Error in FUN(newX[, i], ...) :
`param' has no names and does not match number of coefficients of
model. Unable to construct coefficient vector
Thanks,
Gang
On Oct 30, 2007, at 9:08 AM, Dieter Menne wrote:
Gang Chen gangchen at mail.nih.gov writes:
Suppose I have a mixed-effects model where yij is the jth sample
for
the ith subject:
yij= beta0 + beta1(age) + beta2(age^2) + beta3(age^3) + beta4(IQ) +
beta5(IQ^2) + beta6(age*IQ) + beta7(age^2*IQ) + beta8(age^3 *IQ)
+random intercepti + eij
In R how can I get an F test against the null hypothesis of
beta6=beta7=beta8=0? In SAS I can run something like contrast
age*IQ
1, age^2*IQ 1, age^3 *IQ 1, but is there anything similar in R?
Check packages multcomp and gmodels for contrast tests that work
with lme.
Dieter
__
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] textplot() in gplots causes problems (0x9)
Hi Jonas,
You are correct that manually setting cex won't entirely remove the
warning messages (note: these are *warnings* not *errors*). It will
reduce them greatly, from, perhaps 40, down to 3 or so.
Since the tab character won't be properly displayed anyway, you can
use gsub or similar to replace it with a space. EG:
pdf(file=C:/..., paper=a4, width=8, height=12)
.model - lm(.model.formula, data=database)
text - capture.output(summary(.model))
text - gsub('\t',' ', text)
textplot(, valign=top, halign=left)
dev.off()
As the inclusion of tab characters seems likely to be a problem for
mony people, I've just added code the gplots::textplot() that will
replace all occurences of tab with an appropriate number of spaces
before attempting to compute text heigh or width. This will be part
of the next release of gplots, which should show up on CRAN shortly.
-Greg
On Nov 1, 2007, at 2:18PM , Jonas Malmros wrote:
Dear Gregory,
How can I avoid using tab character when all I want to do is to print
a model summary on my pdf device using textplot()?
How do I set the font size? If you mean using cex inside textplot,
then it does not work. Whether cex is 1 or 0.2 I get the same result,
exemplified here:
Call:
lm(formula =...)
Resuduals:
...
Coefficients Estimate
(intercept)1.32
......
......
Std.Error t-value
(Intercept)0.2 0.1
.........
.........
Pr(|t|)
(intercept) 0.01
......
......
−−−
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.748...
Is there no solution to this problem?
I am using Vista, R2.6.0 patched, RWinEdt.
textplot(capture.output(summary(.model)), valign=top,
halign=left, cex=0.5)
Thanks in advance,
Jonas
On 11/1/07, Gregory Warnes [EMAIL PROTECTED] wrote:
Hi Jonas,
By default, textplot() attempts to automatically select a font size
that is 'just big enough, but not too big'. It does this by a
binary-
search approach where it sets a font size, then asks R to compute the
actual width of the text to be displayed (without actually displaying
it), then increases or decreases the font size appropriately until it
finds the largest font that doesn't extend beyond the plot region
vertically or horizontally. It appears that on your system, R
doesn't know how wide a tab character is. This isn't particularly
surprising since tab characters vary in width depending on the
context.
There are two simple solutions. First, avoid using characters R
can't figure out sizes for (i.e. tab), or manually specify the font
size so textplot() doesn't attempt to optimize it.
I personally choose the former, avoid tab characters, since the
appropriate font size varies greatly by device.
-Greg
On Oct 31, 2007, at 3:22PM , Jonas Malmros wrote:
Hello,
I am using textplot function in gplots package to put some model
output inside a PDF file, but it does not seem to work properly with
PDF.
I am doing follwing:
pdf(file=C:/..., paper=a4, width=8, height=12)
.model - lm(.model.formula, data=database)
textplot(capture.output(summary(.model)), valign=top,
halign=left)
I am getting these error messages:
Warning messages:
1: In FUN(c(C, a, l, l, :, l, m, (, f, o,
r, :
font width unknown for character 0x9
2: In strwidth(object, cex = cex) : font width unknown for
character 0x9
3: In FUN(c(C, a, l, l, :, l, m, (, f, o,
r, :
font width unknown for character 0x9
4: In strwidth(object, cex = cex) : font width unknown for
character 0x9
5: In FUN(c(C, a, l, l, :, l, m, (, f, o,
r, :
font width unknown for character 0x9
6: In strwidth(object, cex = cex) : font width unknown for
character 0x9
7: In text.default(x = xpos, y = ypos, labels = object, adj = c
(0, :
font width unknown for character 0x9
8: In text.default(x = xpos, y = ypos, labels = object, adj = c
(0, :
font width unknown for character 0x9
This is a tab character that causes problems, I guess. Is there any
way to solve this?
Thank you in advance
--
Jonas Malmros
Stockholm University
Stockholm, Sweden
__
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.
--
Jonas Malmros
Stockholm University
Stockholm, Sweden
mime-attachment.txt
__
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] [R-pkgs] SASxport v. 1.2.0
SASxport Version 1.2.0 is now available --- The SASxport package provides R with full support for reading and writing SAS xport format files. Version 1.2.0 corrects a critical issues with storage of negative numbers, as well as adding additional improvements to the handling of SAS display and input format specifications. With these enhancements, both reading and writing SAS transport files is almost lossless*. SASxport 1.2.0 is available _now_ at http://random-technologies-llc.com/products/SASxport/ and will become available on your favorite CRAN repository shortly. For additional information and commercial support packages, please visit the Random Technologies web site: http://random-technologies-llc.com/products/SASxport/ -Greg Gregory R. Warnes, Ph.D Chief Scientist Random Technologies, LLC. (*) Only information that cannot be stored in a SAS xport format file is lost: variable names are uppercased and truncated at 8 characters, missing values for character string variables are converted to empty strings, and numeric values outside of the SAS xport format numeric floating point range are converted to missing values. ___ R-packages mailing list [EMAIL PROTECTED] https://stat.ethz.ch/mailman/listinfo/r-packages __ 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] fit.contrast and interaction terms
Hello Berta, gmodels::fit.contrasts() simply performs a single-variable contrast using the model you have specified. To perform the more involved contrasts that you are describing, there are two approaches: 1) use the estimable() function in the gmodels package. gmodels::estimable() allows you to provides an arbitrary function to be applied to the estimated model parameters, allowing you to perform any appropriate (or inappropriate) contrast calculation. 2) Use the contrast functions provided by Frank Harrell's Hmisc package. Frank's functions allow you to specify the desired value or level of each parameter for the contrast, and handle unspecified parameters. I hope this helps, -Greg On Oct 9, 2007, at 6:10AM , Berta wrote: Dear R-users, I want to fit a linear model with Y as response variable and X a categorical variable (with 4 categories), with the aim of comparing the basal category of X (category=1) with category 4. Unfortunately, there is another categorical variable with 2 categories which interact with x and I have to include it, so my model is s reg3: Y=x*x3. Using fit.contrast to make the contrast (category 1 vs category 4) with options(contrasts=c (contr.treatment, contr.poly)), it makes the contrast but just for the basal category of x3, (coincident with the corresponding result of summary(reg3)), so that it is not what I am looking for, and it seems that when I write (contrasts=c(contr.sum, contr.poly)) before fit.contrast, it adjust for x3. Here I send a SMALL EXAMPLE that tries to imiate my problem. library(gmodels) set.seed(100) options(contrasts=c(contr.treatment, contr.poly)) y - rnorm(100) x - cut(rnorm(100, mean=y, sd=0.25),c(-4,-1.5,0,1.5,4)) # 4 CATEGORIES x3 - cut(rnorm(100, mean=y, sd=8),c(-50,0,50)) # 2 CATEGORIES reg3-lm(y~ x * x3 ) summary(reg3) # category 1 vs category 4 estimate: 3.84776 , for basal category of x3 fit.contrast(reg3,x,c(-1,0,0,1)) # g4 vs g1 # category 1 vs category 4 estimate: 3.84776 , options(contrasts=c(contr.sum, contr.poly)) fit.contrast(reg3,x,c(-1,0,0,1)) # g4 vs g1 # category 1 vs category 4 estimate: 4.010439 I deduce that the global comparison of category 1 vs category 4 is 4.01, and not 3.84. Unfortunately, the real variable that interact with x is not categorical but continuous in my real case, and the new model is Y=x*x2. Again, with the contr.treatment option, the comparison of category 1 vs category 4 is done for the intercept, that is, for the numerical variable assumed to be 0. As i am interested in comparing category 1 vs category 4 adjusting for x2 in general terms, I use contr.sum as before, but it does not seem to produce any modification: x2 - rnorm(100,mean=y,sd=0.5) # NUMERIC reg2 - lm(y ~ x * x2 ) summary(reg2) # category 1 vs category 4 estimate: 3.067346 for x2=0 options(contrasts=c(contr.treatment, contr.poly)) fit.contrast(reg2,x,c(-1,0,0,1)) # g4 vs g1 estimate: 3.067346 for x2=0 options(contrasts=c(contr.sum, contr.poly)) fit.contrast(reg2,x,c(-1,0,0,1)) # g4 vs g1 estimate: 3.067346 for x2=0 The question is: How could I contrast simulaneously in global terms (not intercept and slope separately) if there are differences among category 1 vs category 4 but adjusting for a continuous variable? Or it is not possible to do so, and I have to contrast both (difference of intercepts and difference of slopes separately) and obtain a global conclussion? Thanks a lot in advance, Berta [[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] R and FDA trials
Hi All, I'm excited to see that R-Core has released the document R: Regulatory Compliance and Validation Issues A Guidance Document for the Use of R in Regulated Clinical Trial Environments (http://www.r-project.org/doc/R-FDA.pdf). I know it represents a great deal of effort on the part of R core, and it has taken many months to achieve. Random Technologies, LLC (http://random-technologies-llc.com) has been extended this work to cover additional R packages, as well as offering additional documentation, products, and tools that facilitate use of R for mission-critical and regulated contexts. (Lots of this work falls into the Make your *boss* happy context, grin). These efforts are supported by the sale of support contracts and related services. We are very interested in working with R users and developers to identify additional 1) packages, 2) features, and 3) documents which will (with appropriate documentation and testing, of course) make R easier to use (both functionally and politically) for FDA sumbissions. Please send us your suggestions and comments at [EMAIL PROTECTED] technologies-llc.com -G Gregory R. Warnes, Ph.D. Chief Scientist Random Technologies, LLC. http://random-technologies-llc.com On Oct 8, 2007, at 8:54AM , Marc Schwartz wrote: On Sun, 2007-10-07 at 21:25 -0500, Frank E Harrell Jr wrote: Ricardo Pietrobon wrote: Yesterday I just noticed the new document on R and regulatory aspects for biomedical research posted at http://www.r-project.org/doc/R-FDA.pdf Coming from an institution that performs a large number of clinical trials for FDA and being an advocate of R myself, I have found that the following issues usually come up when discussing the use of R for FDA trials: 1. Most FDA submissions come down to a series of r x k tables, and it is hard to claim that one system is better than another for that. 2. Data is to be submitted to the FDA in SAS (considered by many as the industry standard) or CDISC XML formats (http://www.cdisc.org/); there are pretty good SAS tools for that; does R have comparable? 3. Some packages in R provide acknowledgedly better functionality than their SAS-equivalent, but an entire FDA validation would have to occur each time an enhancement is made to the R package because often an enhancement breaks something else or the syntax would change from one release to another. Your item 3. is up to the company's policy. FDA does not require it and the word validation is not well defined in this context. No package does a complete validation any time any piece of the package is enhanced. Frank Just to augment Frank's response: With respect to point number 2, I believe the only R based package that offers the ability to write SAS XPORT format files at this point is the SASxport CRAN package by Greg Warnes et al. There are several, including Frank's Hmisc package and the standard foreign package that will read XPORT files. Keep in mind that it is the SAS XPORT files, not the proprietary SAS files, that the FDA wants, since the XPORT format is openly documented. There are other third party packages, such as DBMS/Copy and Stat/Transfer, that will also write XPORT format files. At some point, CDISC based exports will likely need to be looked at by one or more folks motivated to do so. In the mean time, you would need to consider third party tools for that. My recollection of Mat Soukup's comments from our panel at useR back in August, is that the FDA itself is not yet at the point where CDISC based file exchanges are without issue. Finally, bear in mind that the document to which you refer covers specifically listed packages released via the R Foundation. See Section 2 of the document. Other CRAN/non-CRAN R packages are not covered. In either case, you would need to engage in appropriate risk mitigation processes when using those and as Frank notes, these would be internal procedures and policies as deemed appropriate within the GxP framework. HTH, Marc Schwartz __ 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.
