Re: [R] coxph means not equal to means of model matrix
I looked at the nocenter and it says (-1,0,1) values but it seems that any
three-level factor is included in that (represented as 1,2,3 in R) .
Also, is the baseline curve now showing the reference level and not the
fictional .428 sex? If I predict the risk for a new row, should I multiply the
coefficient shown in the output by 1 for a sex=1? It used to be (1-.428)*coef.
Thanks for clarifying.
SB
From: Therneau, Terry M., Ph.D.
Sent: Friday, 3 September, 2021 12:37
To: Bond, Stephen
Cc: R-help
Subject: Re: coxph means not equal to means of model matrix
[EXTERNAL]
See ?coxph, in particular the new "nocenter" option.
Basically, the "mean" component is used to center later computations. This can
be critical for continuous variables, avoiding overflow in the exp function,
but is not necessary for 0/1 covariates. The fact that the default survival
curve would be for a sex of .453, say, was off-putting to many.
Terry T.
On 9/3/21 11:01 AM, Bond, Stephen wrote:
Hi,
Please, help me understand what is happening with the means of a Cox model?
I have:
R version 4.0.2 (2020-06-22) -- "Taking Off Again"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)
getOption("contrasts")
unordered ordered
"contr.treatment" "contr.poly"
According to the help coxph.object has a component holding the means of the X
(model.matrix). This does not hold any more.
```
library(survival)
test1 <- list(time=c(4,3,1,1,2,2,3),
status=c(1,1,1,0,1,1,0),
x=c(0,2,1,1,1,0,0),
sex=factor(c(0,0,0,0,1,1,1)))
m1 <- coxph(Surv(time, status) ~ x + sex, test1)
m1$means
##x sex1
## 0.7142857 0.000
colMeans(model.matrix(m1))
## x sex1
## 0.7142857 0.4285714
```
Will new observations be scored using the zero mean from the object?? Is this
just a reporting change where $means shows the reference level and no longer
the mean of the model matrix??
Thanks everybody
ATTENTION : This email originated outside your organization. Exercise caution
before clicking links, opening attachments, or responding with personal
information.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] coxph means not equal to means of model matrix
Hi,
Please, help me understand what is happening with the means of a Cox model?
I have:
R version 4.0.2 (2020-06-22) -- "Taking Off Again"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)
getOption("contrasts")
unordered ordered
"contr.treatment" "contr.poly"
According to the help coxph.object has a component holding the means of the X
(model.matrix). This does not hold any more.
```
library(survival)
test1 <- list(time=c(4,3,1,1,2,2,3),
status=c(1,1,1,0,1,1,0),
x=c(0,2,1,1,1,0,0),
sex=factor(c(0,0,0,0,1,1,1)))
m1 <- coxph(Surv(time, status) ~ x + sex, test1)
m1$means
##x sex1
## 0.7142857 0.000
colMeans(model.matrix(m1))
## x sex1
## 0.7142857 0.4285714
```
Will new observations be scored using the zero mean from the object?? Is this
just a reporting change where $means shows the reference level and no longer
the mean of the model matrix??
Thanks everybody
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[ESS] ESS adds symbols at beginning of line
Hello ESS friends,
ESS started to put "> " after the initial "> " at the beginning of the line.
Sometimes it works even though it looks strange, but sometimes it stops working
and I have to copy from the source buffer and paste in my R buffer to execute
lines, which is bad. Any ideas how to get rid of the unwanted symbols? Looks
like this:
> > value.date <- as.Date('2020-11-30')
I send lines by C-c C-n from the source buffer.
My ESS buffer has:
[ess-site:] ess-lisp-directory =
'c:/emacs/emacs-26.3/site-lisp/ess-18.10.2/lisp'[ess-site.el]:
ess-customize-alist=nil
[ess-site.el _2_]: ess-customize-alist=nil
Creating global Emacs toolbar[ess-site:] *very* end ...(ess-mode-1):
ess-language=S, ess-dialect=R buf=scratch.R
(ess-mode-1.5): ;;very long output omitted
(R): ess-dialect=R, buf=scratch.R, start-arg=nil
current-prefix-arg=nil
(inf-ess 1): lang=S, dialect=R, tmp-dialect=R, buf=scratch.R
(inf-ess 1.1): procname=R temp-dialect=R, buf-name=*R*
(inf-ess 2.0) Method #3 start=c:/emacs/New buf=*R*
(inf-ess 2.1): ess-language=S, ess-dialect=R buf=*R*
(i-ess 1): buf=*R*, lang=S, comint..echo=t,
comint..sender=inferior-ess-input-sender,
(i-ess end): buf=*R*, lang=S, comint..echo=t,
comint..sender=inferior-ess-input-sender,
(inf-ess 3.0): prog=Rterm, start-args=--ess , echoes=t
Making Process...Buf *R*, :Proc R, :Prog Rterm
Thanks everybody.
[[alternative HTML version deleted]]
__
ESS-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/ess-help
Re: [R] how to add a child to a child in XML
Just to clarify and hopefully catch the attention of the maintainer:
The newXMLNode function is not mentioned in:
https://cran.r-project.org/web/packages/XML/XML.pdf
which supposedly describes all functions in the package.
Stephen
From: Ben Tupper [mailto:btup...@bigelow.org]
Sent: Thursday, March 22, 2018 10:40 AM
To: Bond, Stephen
Cc: r-help
Subject: Re: [R] how to add a child to a child in XML
Hi,
It's a reasonable question. The answer is that it actually is included, but
there are many instances across packages where multiple functions are
documented on a single help page. The following brings up such a page... (for
XML_3.98-1.9)
> library(XML)
> ?newXMLNode
You can see the same on line...
https://www.rdocumentation.org/packages/XML/versions/3.98-1.9/topics/newXMLDoc
You have dig in to find it.
If you are just starting out with XML, you might want to spend some time
comparison shopping with the xml2 package.
https://www.rdocumentation.org/packages/xml2/versions/1.2.0 I like each one,
and I use both XML and xml2 (not at the same time). I have been slowly
migrating toward xml2 as I use more of the tidyverse stuff.
Cheers,
Ben
On Mar 22, 2018, at 9:19 AM, Bond, Stephen
<stephen.b...@cibc.com<mailto:stephen.b...@cibc.com>> wrote:
Big thanks. newXMLNode works great. Wonder why it is not included in the
documentation.
There is newXMLDoc and newXMLNamespace, but no mention of newXMLNode.
Stephen
From: Ben Tupper [mailto:btup...@bigelow.org]
Sent: Wednesday, March 21, 2018 6:18 PM
To: Bond, Stephen
Cc: r-help
Subject: Re: [R] how to add a child to a child in XML
Hi,
XML doesn't use the `$` to access child nodes. Instead use either `[name]` to
get a list of children of that name or `[[name]]` to get the just the first
child of that name encountered in the genealogy. Thus for your example...
> root$child1
NULL
> root[['child1']]
On the other hand, you might consider using newXMLNode() instead of xmlNode()
as it accepts a "parent = " argument. The alternative using newXMLNode() would
look like...
atts_root <- c("val1","val2","val3")
names(atts_root) <- c("att1","att2","att3")
root <- newXMLNode("root", attrs = atts_root)
atts_child <- LETTERS[1:3]
names(atts_child) <- paste("name",1:3,sep="")
child <- newXMLNode("child",attrs = atts_child, parent = root)
atts_grandchild <- letters[1:3]
names(atts_grandchild) <- paste("name",4:6,sep="")
grandchild <- newXMLNode("grandchild",attrs = atts_grandchild, parent = child)
root
#
#
#
#
#
Cheers,
Ben
On Mar 21, 2018, at 4:25 PM, Bond, Stephen
<stephen.b...@cibc.com<mailto:stephen.b...@cibc.com>> wrote:
I am trying to add a child to a child using XML package in R. the following
fails
library(XML)
node1 <- c("val1","val2","val3")
names(node1) <- c("att1","att2","att3")
root <- xmlNode("root", attrs=node1)
node2 <- LETTERS[1:3]
names(node2) <- paste("name",1:3,sep="")
root <- addChildren(root,xmlNode("child1",attrs=node2))
node3 <- letters[1:3]
names(node3) <- paste("name",4:6,sep="")
root <- addChildren(root$child1,xmlNode("child2",attrs=node3))
Error in UseMethod("addChildren") : no applicable method for 'addChildren'
applied to an object of class "NULL"
Stephen B
[[alternative HTML version deleted]]
__
R-help@r-project.org<mailto:R-help@r-project.org> mailing list -- To
UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html<http://www.r-project.org/posting-guide.html>
and provide commented, minimal, self-contained, reproducible code.
Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org<http://www.bigelow.org/>
Tick Forecasting: https://eco.bigelow.org/
Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org
Tick Forecasting: https://eco.bigelow.org/
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to add a child to a child in XML
Big thanks. newXMLNode works great. Wonder why it is not included in the
documentation.
There is newXMLDoc and newXMLNamespace, but no mention of newXMLNode.
Stephen
From: Ben Tupper [mailto:btup...@bigelow.org]
Sent: Wednesday, March 21, 2018 6:18 PM
To: Bond, Stephen
Cc: r-help
Subject: Re: [R] how to add a child to a child in XML
Hi,
XML doesn't use the `$` to access child nodes. Instead use either `[name]` to
get a list of children of that name or `[[name]]` to get the just the first
child of that name encountered in the genealogy. Thus for your example...
> root$child1
NULL
> root[['child1']]
On the other hand, you might consider using newXMLNode() instead of xmlNode()
as it accepts a "parent = " argument. The alternative using newXMLNode() would
look like...
atts_root <- c("val1","val2","val3")
names(atts_root) <- c("att1","att2","att3")
root <- newXMLNode("root", attrs = atts_root)
atts_child <- LETTERS[1:3]
names(atts_child) <- paste("name",1:3,sep="")
child <- newXMLNode("child",attrs = atts_child, parent = root)
atts_grandchild <- letters[1:3]
names(atts_grandchild) <- paste("name",4:6,sep="")
grandchild <- newXMLNode("grandchild",attrs = atts_grandchild, parent = child)
root
#
#
#
#
#
Cheers,
Ben
On Mar 21, 2018, at 4:25 PM, Bond, Stephen
<stephen.b...@cibc.com<mailto:stephen.b...@cibc.com>> wrote:
I am trying to add a child to a child using XML package in R. the following
fails
library(XML)
node1 <- c("val1","val2","val3")
names(node1) <- c("att1","att2","att3")
root <- xmlNode("root", attrs=node1)
node2 <- LETTERS[1:3]
names(node2) <- paste("name",1:3,sep="")
root <- addChildren(root,xmlNode("child1",attrs=node2))
node3 <- letters[1:3]
names(node3) <- paste("name",4:6,sep="")
root <- addChildren(root$child1,xmlNode("child2",attrs=node3))
Error in UseMethod("addChildren") : no applicable method for 'addChildren'
applied to an object of class "NULL"
Stephen B
[[alternative HTML version deleted]]
__
R-help@r-project.org<mailto:R-help@r-project.org> mailing list -- To
UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org
Tick Forecasting: https://eco.bigelow.org/
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] how to add a child to a child in XML
I am trying to add a child to a child using XML package in R. the following
fails
library(XML)
node1 <- c("val1","val2","val3")
names(node1) <- c("att1","att2","att3")
root <- xmlNode("root", attrs=node1)
node2 <- LETTERS[1:3]
names(node2) <- paste("name",1:3,sep="")
root <- addChildren(root,xmlNode("child1",attrs=node2))
node3 <- letters[1:3]
names(node3) <- paste("name",4:6,sep="")
root <- addChildren(root$child1,xmlNode("child2",attrs=node3))
Error in UseMethod("addChildren") : no applicable method for 'addChildren'
applied to an object of class "NULL"
Stephen B
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] offset with a factor
Knowledgeable useRs, Please, advise how to use offset with a factor. I estimate monthly effects from a much bigger data set as monthly effects seem to be stable, and other variables are estimated from a small, but recent data set as there is variation in those non-seasonal coefficients. How can I use the seasonality estimates from the big data set as an offset provided to the small data set. I know an offset is supposed to be quantitative, but this is such a practical and sensible scenario, I feel compelled. Assume I have 11 coefs estimated with contr.sum. Thanks everybody Stephen [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [ESS] connection to Oracle does not work in ESS
ESS or emacs stores sth between sessions! I tried quitting the R session and
connecting from a new session in the same emacs window and it did not work.
Today I closed the whole emacs window and reopened the same R code and it
worked.
Stephen
-Original Message-
From: Vitalie Spinu [mailto:spinu...@gmail.com]
Sent: Wednesday, December 13, 2017 5:34 PM
To: Bond, Stephen
Cc: 'ess-help@r-project.org'
Subject: Re: [ESS] connection to Oracle does not work in ESS
Hi Stephen,
I doubt it's problem with ESS. Running R from emacs shell and with ESS should
be equivalent. Make sure that it's indeed the same R instance and RODBC what
are called in your experiments. Do you have a custom
inferior-ess-r-program-name setup?
Do you issue Rterm or R in the emacs shell?
Vitalie
>> On Wed, Dec 13 2017 22:23, Bond, Stephen wrote:
> Further info:
> When I use Rterm in a cmd shell the connection works.
> When I use Rterm in a shell within Emacs the connection works.
> Something must be going wrong in the way ESS does things.
> Stephen
> From: Bond, Stephen
> Sent: Wednesday, December 13, 2017 12:53 PM
> To: 'ess-help@r-project.org'
> Subject: connection to Oracle does not work in ESS
> On 64-bit win7:
> I have a DSN which tests successfully in control panel. I can connect from
> RGui with:
> con=odbcConnect("DSNname",uid="myuid",pwd="mypass")
> it does not work in ESS:
> R version 3.4.2 (2017-09-28) -- "Short Summer"
> Copyright (C) 2017 The R Foundation for Statistical Computing
> Platform: x86_64-w64-mingw32/x64 (64-bit)
>> > options(chmhelp=FALSE, help_type="text")
>> options(STERM='iESS', str.dendrogram.last="'",
>> editor='emacsclient.exe', show.error.locations=TRUE)
>> library(RODBC)
>> con=odbcConnect("DSNname",uid="myuid",pwd="mypass")
> Warning messages:
> 1: In RODBC::odbcDriverConnect() :
> [RODBC] ERROR: state 08004, code 12154, message
> [Oracle][ODBC][Ora]ORA-12154: TNS:could not resolve the connect
> identifier specified
> 2: In RODBC::odbcDriverConnect() :
> ODBC connection failed
> I have all the paths setup and the driver is installed correctly. Tnsping
> finds the db.
> Let me know if anybody can help.
> Stephen
> [[alternative HTML version deleted]]
> __
> ESS-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/ess-help
__
ESS-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/ess-help
[ESS] connection to Oracle does not work in ESS
On 64-bit win7:
I have a DSN which tests successfully in control panel. I can connect from RGui
with:
con=odbcConnect("DSNname",uid="myuid",pwd="mypass")
it does not work in ESS:
R version 3.4.2 (2017-09-28) -- "Short Summer"
Copyright (C) 2017 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)
> > options(chmhelp=FALSE, help_type="text")
> options(STERM='iESS', str.dendrogram.last="'", editor='emacsclient.exe',
> show.error.locations=TRUE)
> library(RODBC)
> con=odbcConnect("DSNname",uid="myuid",pwd="mypass")
Warning messages:
1: In RODBC::odbcDriverConnect() :
[RODBC] ERROR: state 08004, code 12154, message [Oracle][ODBC][Ora]ORA-12154:
TNS:could not resolve the connect identifier specified
2: In RODBC::odbcDriverConnect() :
ODBC connection failed
I have all the paths setup and the driver is installed correctly. Tnsping finds
the db.
Let me know if anybody can help.
Stephen
[[alternative HTML version deleted]]
__
ESS-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/ess-help
Re: [R] update R version in windows
Thanks Caitlin and Richard MH. Works great.
Stephen
From: Caitlin [mailto:bioprogram...@gmail.com]
Sent: Friday, November 10, 2017 12:33 PM
To: Bond, Stephen
Subject: Re: [R] update R version in windows
install.packages("installr")
updateR()
rather...
On Friday, November 10, 2017, Bond, Stephen
<stephen.b...@cibc.com<mailto:stephen.b...@cibc.com>> wrote:
Is there a utility which will allow me to upgrade my R version and update all
packages from the old version?
If I manually upgrade, then I have to manually re-install 50 packages.
Thank you.
Stephen B
[[alternative HTML version deleted]]
__
R-help@r-project.org<javascript:;> mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] update R version in windows
This issue does not exist on Linux. My Ubuntu updates both R and all packages. Stephen B -Original Message- From: J C Nash [mailto:profjcn...@gmail.com] Sent: Friday, November 10, 2017 1:19 PM To: r-help; RICHARD M. HEIBERGER; Bond, Stephen Subject: Re: [R] update R version in windows However, trying this on Linux Mint gave package ‘installr’ is not available (for R version 3.4.2) Has the package not been updated yet? JN Try the installr package. It was designed for this purpose. On Fri, Nov 10, 2017 at 11:49 AM, Bond, Stephen wrote: > Is there a utility which will allow me to upgrade my R version and update all > packages from the old version? > If I manually upgrade, then I have to manually re-install 50 packages. > Thank you. > > Stephen B > > __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] update R version in windows
Is there a utility which will allow me to upgrade my R version and update all packages from the old version? If I manually upgrade, then I have to manually re-install 50 packages. Thank you. Stephen B [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] fedback from foreach
Hi useRs,
I am running a foreach loop and hoped to get a small message when it hits a
multiple of 1000, but it does not work.
p <- foreach(i=1:1, .combine='c') %dopar% {
if(i%%1000==0) print(i)
sqrt(i)
}
What is the proper way to do it.
Thanks everybody.
Stephen B
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] minimal attributes to get se.fit
Does anybody know what are the attributes of a glm fit object that will allow the "predict call" to produce an se.fit? I am deleting most of the attributes as the size of the final object is 5Gb and I want to reduce it to under 20Mb, but that causes as error when I ask for an se.fit . mod.b$fitted.values <- 1:10 mod.b$prior.weights <- 1:10 mod.b$data <-mod.b$data[1:10,] mod.b$residuals <- 1:10 mod.b$linear.predictors <- 1:10 mod.b$qr$qr <- mod.b$qr$qr[1:10,] mod.b$effects <- mod.b$effects[1:100] mod.b$weights <- mod.b$weights[1:100] mod.b$model <- mod.b$model[1:10,] mod.b$y <- mod.b$y[1:10] p1 <- predict(mod.b,new=newdata,type="link",se.fit=T) Error in Qr$qr[p1, p1, drop = FALSE] : subscript out of bounds I believe the covariance matrix of the coefficients is all that should be needed and that is quite small. However, the covariance matrix is not an attribute of the model object. Thanks everybody. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] lattice add a fit
Folks,
This is just misunderstanding. I did not want a panel function for locfit. In
my email I say:
Instead, I want to put a fit from lm (but
not a simple straight line) and the fit has to be done for each panel
separately, not one fit for the full data set, so sth like an lm
equivalent of panel.locfit (there is no panel.lmfit) Thank you.
Bert Gunter provided the answer to my question. Maybe I should have sent a
thank you note to the list to finalize.
Kind regards
Stephen Bond
-Original Message-
From: David Winsemius [mailto:dwinsem...@comcast.net]
Sent: Wednesday, October 08, 2014 12:30 PM
To: Duncan Mackay
Cc: R; Bond, Stephen
Subject: Re: [R] lattice add a fit
On Oct 7, 2014, at 9:15 PM, Duncan Mackay wrote:
I'm a tad puzzled by the comments about needing to build a panel function for
locfit. The various plot.locfit functions are actually lattice calls.
locfit:::panel.locfit # already exists, and even has versions for 1d, 2d and
3d purposes.
And there is a llines.locfit function that will add locfit smooths to existing
lattice plots.
It's a very simple function and could easily be modified to any regression
method that has a predict functions:
locfit:::llines.locfit
function (x, m = 100, tr = x$trans, ...) {
newx - lfmarg(x, m = m)[[1]] # probably need to modify to your purposes
y - predict(x, newx, tr = tr)
llines(newx, y, ...)
}
environment: namespace:locfit
--
David
You will have to make your own panel function for locfit if you want
to use it I have done it in the past - read the help for
library(locfit)
?plot.locfit
and the links
?lattice::prepanel
Regards
Duncan
Duncan Mackay
Department of Agronomy and Soil Science University of New England
Armidale NSW 2351
Email: home: mac...@northnet.com.au
-Original Message-
From: r-help-boun...@r-project.org
[mailto:r-help-boun...@r-project.org] On Behalf Of Bond, Stephen
Sent: Tuesday, 7 October 2014 23:02
To: r-help@R-project.org
Subject: [R] lattice add a fit
What is the way to add an arbitrary fit from a model to a lattice
conditioning plot ?
For example
xyplot(v1 ~v2 | v3,data=mydata,
panel=function(...){
panel.xyplot(...)
panel.loess(...,col.line=red)
}
)
Will add a loess smoother. Instead, I want to put a fit from lm (but
not a simple straight line) and the fit has to be done for each panel
separately, not one fit for the full data set, so sth like an lm
equivalent of panel.locfit (there is no panel.lmfit) Thank you.
Stephen B
[[alternative HTML version deleted]]
David Winsemius
Alameda, CA, USA
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] lattice add a fit
What is the way to add an arbitrary fit from a model to a lattice conditioning
plot ?
For example
xyplot(v1 ~v2 | v3,data=mydata,
panel=function(...){
panel.xyplot(...)
panel.loess(...,col.line=red)
}
)
Will add a loess smoother. Instead, I want to put a fit from lm (but not a
simple straight line) and the fit has to be done for each panel separately, not
one fit for the full data set, so sth like an lm equivalent of panel.locfit
(there is no panel.lmfit)
Thank you.
Stephen B
[[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 add a fit
Bert,
Can you provide an example how to pass the conditioned data set?
xyplot( x~y | z, data=mydata,
panel=function(...){
mod=lm(x~ y+w +q, data=??)
panel.lines(fitted(mod))
}
If I use mydata in place of ?? I get a global fit, not a fit for each level of
z, which is what I want.
Stephen B
-Original Message-
From: Bert Gunter [mailto:gunter.ber...@gene.com]
Sent: Tuesday, October 07, 2014 9:30 AM
To: Bond, Stephen
Cc: r-help@R-project.org
Subject: Re: [R] lattice add a fit
Fit your model in the panel function using lm and plot the fits using
?panel.points, ?panel.lines, etc.
-- Bert
Bert Gunter
Genentech Nonclinical Biostatistics
(650) 467-7374
Data is not information. Information is not knowledge. And knowledge is
certainly not wisdom.
Clifford Stoll
On Tue, Oct 7, 2014 at 6:01 AM, Bond, Stephen stephen.b...@cibc.com wrote:
What is the way to add an arbitrary fit from a model to a lattice
conditioning plot ?
For example
xyplot(v1 ~v2 | v3,data=mydata,
panel=function(...){
panel.xyplot(...)
panel.loess(...,col.line=red)
}
)
Will add a loess smoother. Instead, I want to put a fit from lm (but
not a simple straight line) and the fit has to be done for each panel
separately, not one fit for the full data set, so sth like an lm equivalent
of panel.locfit (there is no panel.lmfit) Thank you.
Stephen B
[[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] Rscript fails where Rterm works
Hadley, You are a genius. Stephen B -Original Message- From: Hadley Wickham [mailto:h.wick...@gmail.com] Sent: Thursday, June 12, 2014 5:18 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] Rscript fails where Rterm works Explicitly load the methods package: library(methods) Hadley On Thu, Jun 12, 2014 at 2:22 PM, Bond, Stephen stephen.b...@cibc.com wrote: I have a script which loads library(XLConnect) wb - loadWorkbook(wbname) the code works without errors when run from ESS which uses R version 3.0.1 (2013-05-16) -- Good Sport Copyright (C) 2013 The R Foundation for Statistical Computing Platform: i386-w64-mingw32/i386 (32-bit) But fails when run from Rscript which Rscript C:\Program Files\R\R-3.0.1\bin\i386\Rscript.EXE Rscript updatevols.R Warning message: package 'RODBC' was built under R version 3.0.2 XLConnect 0.2-7 by Mirai Solutions GmbH http://www.mirai-solutions.com , http://miraisolutions.wordpress.com Warning message: package 'XLConnect' was built under R version 3.0.3 Error in .jfield(x, Ljava/lang/Class;, TYPE) : could not find function getClass Calls: loadWorkbook ... ._java_class_list - lapply - FUN - .jfield - .Call Execution halted Warning message: closing unused RODBC handle 1 Any clues? Stephen B [[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. -- http://had.co.nz/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Rscript fails where Rterm works
I have a script which loads library(XLConnect) wb - loadWorkbook(wbname) the code works without errors when run from ESS which uses R version 3.0.1 (2013-05-16) -- Good Sport Copyright (C) 2013 The R Foundation for Statistical Computing Platform: i386-w64-mingw32/i386 (32-bit) But fails when run from Rscript which Rscript C:\Program Files\R\R-3.0.1\bin\i386\Rscript.EXE Rscript updatevols.R Warning message: package 'RODBC' was built under R version 3.0.2 XLConnect 0.2-7 by Mirai Solutions GmbH http://www.mirai-solutions.com , http://miraisolutions.wordpress.com Warning message: package 'XLConnect' was built under R version 3.0.3 Error in .jfield(x, Ljava/lang/Class;, TYPE) : could not find function getClass Calls: loadWorkbook ... ._java_class_list - lapply - FUN - .jfield - .Call Execution halted Warning message: closing unused RODBC handle 1 Any clues? Stephen B [[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] rJava fail
R version 3.1.0 (2014-04-10) -- Spring Dance Copyright (C) 2014 The R Foundation for Statistical Computing Platform: i386-w64-mingw32/i386 (32-bit) library(rJava) Error : .onLoad failed in loadNamespace() for 'rJava', details: call: dirname(this$RuntimeLib) error: a character vector argument expected Error: package or namespace load failed for 'rJava' Things used to work on R 3.0.1 but suddenly stopped. I installed the new R and new packages. Then started downgrading Java. Went from Java 7 to Java 6 update 16 and still no luck. Please, advise which Java I need and if any paths need to be modified. Thank you. Stephen B [[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] windowing
Is there a package or a command that does window aggregation like select sum(col1) over (partition by col2, col3 order by col4 rows between unbounded preceding and current row) as sum1 from table1 ; the above is Netezza syntax, but Postgre has same capability. Stephen B [[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] odfWeave on WinXP
Hello useRs, Having trouble getting odfWeave to work here. Crashing at unzipping like noted in many other posts. The unzip utility referred in the doc no longer exists (page last updated in 2004 and links are not found.) If anybody knows how to make it work on XP pls advise. I can manually unzip from the command line and tried setting that option in the control, but still: odfWeave(examples.odt.zip,outFile,control =odfWeaveControl(zipCmd = c(pkzipc -add $$file$$ ., pkzipc -extract $$file$$))) Creating C:\DOCUME~1\BondStep\LOCALS~1\Temp\RtmpMV18bU/odfWeave04122034507 Copying examples.odt.zip Setting wd to C:\Documents and Settings\BondStep\Local Settings\Temp\RtmpMV18bU\odfWeave04122034507 Unzipping ODF file using pkzipc -extract examples.odt.zip Error in odfWeave(examples.odt.zip, outFile, control = odfc) : Error unzipping file From cygwin: $ pkzipc -extract examples.odt.zip PKZIP(R) Version 8 ZIP Compression Utility for Windows Copyright (C) 1989-2005 PKWARE, Inc. All Rights Reserved. Registered Version PKZIP Reg. U.S. Pat. and Tm. Off. Patent No. 5,051,745 Patent Pending Extracting files from .ZIP: X:\examples.odt.zip Unstoring: Configurations2/floater OK Unstoring: Configurations2/images/Bitmaps OK Unstoring: Configurations2/menubar OK Unstoring: Configurations2/popupmenu OK Unstoring: Configurations2/progressbar OK Unstoring: Configurations2/statusbar OK Unstoring: Configurations2/toolbar OK Inflating: Configurations2/accelerator/current.xml OK Inflating: META-INF/manifest.xml OK Unstoring: Pictures/115F0163AB6E6424.png OK Inflating: Thumbnails/thumbnail.png OK Inflating: content.xml OK Inflating: layout-cache OK Unstoring: meta.xml OK Unstoring: mimetype OK Inflating: settings.xml OK Inflating: styles.xmlOK I put a .zip extension as pkzip chokes on anything which is not .zip Have a great day everybody. Stephen B [[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] unbalanced design
Please, help with a formula for dealing with unbalanced design: To see the counts: aggregate(dfa$CertId,by=list(type=dfa$ComType,stat=dfa$StatusCodeId),length) type stat x 1C1 6571 2C3 28957 3C8 12390 4C 11 12415 5E 13 9 6R 1351 7E 15 2079 8R 15 6692 I would like to have a slope for statuses 1,3,8,11,13 and two slopes for status 15 one for type E and one for type R. I tried nesting, but it assumes that all levels exist for each factor and complains about singular model matrix. Is there a theoretically proper way to deal with this or I should just relabel status 15 and make it 16 for type R and regress on status alone?? Thanks everybody Stephen B __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] obtainl survival curves for single strata
Dr. Therneau,
You are correct about the fit:
(cph.approve -
coxph(Surv(fundterm,resp)~I(CommitAmt/1e5)+res+CommitedRate+dayflag+mth+strata(termfac),
data=dfa,
subset=(HedgeDate2012-11-15 p1!=FixedOpen), method=efron))
However the output from survfit has individuals in each column and the strata
are stacked so, in order to use that I have to subset the exact rows for the
strata I need even though the strata is provided in the newdata argument.
Based on the help file it seems like the function should be able to use the
strata info and return a single strata per subject. I am pasting my current
code, which is not fast due to calls to reshape. If sth similar can be achieved
by calling the compiled code it should run much faster allowing 100k subjects
to be done in 2-3 min.
To use survfit as is, I would need to achieve subsetting in a for loop (going
through columns), which is even slower than reshape. In my old code I processed
subjects one by one in a for loop and that was much slower than grouping all
subjects from the same strata together as in the code below.
surv.approve - survfit(cph.approve)
b1 - c(1,cumsum(surv.approve$strata)+1)
b2 - cumsum(surv.approve$strata)
b1 - b1[-length(b1)]
stratbins -
data.frame(termfac=as.integer(substring(names(b2),9,9)),start=b1,finish=b2)
stratbins
termfac start finish
1 1 1 93
2 294187
3 3 188282
4 4 283372
5 5 373462
6 6 463553
7 7 554633
8 8 634695
9 9 696789
strats - sort(unique(newapp$termfac))
for (jj in strats){
cat('strata ',jj,'\n')
block - newapp[newapp$termfac==jj,]
surv -
surv.approve$surv[stratbins[stratbins$termfac==jj,start]:stratbins[stratbins$termfac==jj,finish]]
risk - predict(cph.approve,new=block,type=risk,ref=sample)
newsurv - outer(surv,risk,^)
days -
as.Date(outer(surv.approve$time[stratbins[stratbins$termfac==jj,start]:stratbins[stratbins$termfac==jj,finish]],
block$HedgeDate,+))
fund - t(t(newsurv)*block$CommitAmt)
fund - rbind(block$CommitAmt,fund)
fund - -diff(fund)
fund - as.data.frame(t(fund))
fund$acct - as.integer(rownames(fund))
ncols - ncol(fund)-1
fundlong -
reshape(fund,dir=long,varying=list(colnames(fund)[1:ncols]),idvar=acct,timevar=daycnt)
fundlong - fundlong[order(fundlong$acct,fundlong$daycnt),]
days - matrix(days,nrow=nrow(days))
days - t(days)
days - cbind(fund[,acct],days)
days - as.data.frame(days)
colnames(days)[1] - acct
daylong -
reshape(days,dir=long,varying=list(colnames(days)[2:(ncols+1)]),idvar=acct,timevar=dt)
daylong - daylong[order(daylong$acct,daylong$dt),]
daylong$V2 - as.Date(daylong$V2,origin = 1970-01-01)
fundlong$dt - as.character(daylong$V2)
sqlSave(ch1,fundlong,tmpfund)
sqlQuery(ch1, insert into RRfund (acct,dt,fund,daycnt) select
acct,dt,V1,daycnt from tmpfund)
sqlQuery(ch1,drop table tmpfund)
}
Output looks like:
acctdt funddaycnt
3623963 2012-11-16 00:00:00 472.99329489343 1
3623963 2012-11-17 00:00:00 5842.482289437712
3623963 2012-11-18 00:00:00 7807.171026727333
3623963 2012-11-19 00:00:00 7244.017697003454
3623963 2012-11-20 00:00:00 7073.831095927985
3623963 2012-11-21 00:00:00 8745.915155128846
3623963 2012-11-22 00:00:00 9473.871528069497
3623963 2012-11-23 00:00:00 12627.88904229168
3623963 2012-11-24 00:00:00 19598.56847130979
3623963 2012-11-25 00:00:00 56094.581213680210
3623963 2012-11-26 00:00:00 25690.439183149 11
3623963 2012-11-27 00:00:00 25420.991525671412
3623963 2012-11-28 00:00:00 30322.375086543413
3623963 2012-11-29 00:00:00 18651.113684675814
3623963 2012-11-30 00:00:00 20291.452806729215
Stephen
-Original Message-
From: Terry Therneau [mailto:thern...@mayo.edu]
Sent: Friday, February 01, 2013 10:11 AM
To: r-help@r-project.org; Bond, Stephen
Subject: Re: obtainl survival curves for single strata
Stephen,
I can almost but not quite get my arms around what you are asking. A bit
more detail
would help. Like
cph.approve = what kind of object, generated by function __
But, let me make a guess:
cph is the result of coxph, and the model has both covariates and a strata
You want predictioned survival curves, more than 1, of the type covariates
= a, b,c,
strata=1 covariates = d,e, f, strata=2, ... for arbitrary covariates and
strata.
Now, the predicted survival curves for different strata are different lengths.
The survfit.coxph routine gets around this by being verbose: it expects you to
give it
covariate sets, and returns
all of the strata for each covariate. This allows it to give a compact result.
You can
always do:
newpred - survfit(cox-model-fit, newdata=something)
newpred[5,17] # survival curve for the 5th strata, covariates from
[R] obtainl survival curves for single strata
Dear useRs, What is the syntax to obtain survival curves for single strata on many subjects? I have a model based on Surv(time,response) object, so there is a single row per subject and no start,stop and no switching of strata. The newdata has many subjects and each subject has a strata and the survival based on the subject risk and the subject strata is needed. If I do newpred - survfit(cph.approve,new=newapp,se=F) I get all strata for every subject. Attempting newpred - survfit(cph.approve,new=newapp,id=CertId,se=F) Error in survfit.coxph(cph.approve, new = newapp, id = CertId, se = F) : The individual option is only valid for start-stop data newpred - survfit(cph.approve,new=newapp,indi=T,se=F) Error in survfit.coxph(cph.approve, new = newapp, indi = T, se = F) : The individual option is only valid for start-stop data Please, advise if obtaining a single strata for a basic (time,response) model is possible? Due to differing lengths of the surv for different strata this will not go in a wide data.frame without padding. Thanks everybody and have a great day. Stephen B __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] coxph with smooth survival
Hello users, I would like to obtain a survival curve from a Cox model that is smooth and does not have zero differences due to no events for those particular days. I have: sum((diff(surv))==0) [1] 18 So you can see 18 days where the survival curve did not drop due to no events. Is there a way to ask survfit to fit a nice spline for the survival?? Note: I tried survreg and it did not work, but maybe I did not do it properly?? Thank you very much. Stephen B __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] coxph with smooth survival
I also tried fitting a spline to the resulting survival curve and the result was horrible. maybe spline won't work or knots need special handling. overall, I must have the final point of the smooth survival to be same as the final point of the raw Cox survival and no flat days, the drops should be spread around. Thanks again everybody. Stephen __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Rterm does not load personal library
Greetings, I am trying to run a short script from a shell: c:\projects\hellRscript --default-packages=mypack X:/4Stephen/commit/curve.R X:/4Stephen/commit/run1.out Loading required package: utils Warning message: package 'RODBC' was built under R version 2.12.2 Error: could not find function normalize Execution halted Warning message: closing unused RODBC handle 1 function 'normalize' is in mypack and I even attempt to load it from inside the script: ## curve.R dt - Sys.Date() library(RODBC) pipe - odbcConnect(commit) library(mypack) # with quotes it's the same result normalize(pipe) any ideas how to make Rterm load the library?? Running from inside ESS works fine. Also from RGui. It does not work if I start a session by R -vanilla on the command line. Thanks and happy holidays. Stephen __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Rterm does not load personal library
Uwe, Thank you for the suggestion. I checked and there are no dependencies involved. I discovered the problem is calling an old version c:\projects\hellwhich Rscript which Rscript /cygdrive/c/Program Files/R/R-2.12.1/bin/Rscript Mypack is installed under 2.15. Will change my paths and that should fix it. Would be nice if Rscript produced a message saying it could not find the lib. It returns nothing so it looks like the load worked, while it did not: c:\projects\hellR --vanilla R --vanilla R version 2.12.1 (2010-12-16) Copyright (C) 2010 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i386-pc-mingw32/i386 (32-bit) library(RODBC) library(RODBC) Warning message: package 'RODBC' was built under R version 2.12.2 library(mypack) library(mypack) normalize normalize Error: object 'normalize' not found Execution halted Stephen -Original Message- From: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de] Sent: Friday, December 21, 2012 1:01 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] Rterm does not load personal library On 21.12.2012 18:16, Bond, Stephen wrote: Greetings, I am trying to run a short script from a shell: c:\projects\hellRscript --default-packages=mypack X:/4Stephen/commit/curve.R X:/4Stephen/commit/run1.out Loading required package: utils Warning message: package 'RODBC' was built under R version 2.12.2 Error: could not find function normalize Execution halted Warning message: closing unused RODBC handle 1 function 'normalize' is in mypack and I even attempt to load it from inside the script: ## curve.R dt - Sys.Date() library(RODBC) pipe - odbcConnect(commit) library(mypack) # with quotes it's the same result normalize(pipe) any ideas how to make Rterm load the library?? Running from inside ESS works fine. Also from RGui. It does not work if I start a session by R -vanilla on the command line. I guess some dependencies on other base R packages are not declared and those packages are not loaded by your package. If you start RGui some packages are loaded that are not loaded if you start via --vanilla, but we cannot know more without information about mypack. Uwe Ligges Thanks and happy holidays. Stephen __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Rterm does not load personal library
I have to apologize to everybody. The thing was failing because my library did not have that particular function in R 2.12. The library actually loads fine, just the function is missing. Merry Christmas to those celebrating. Stephen -Original Message- From: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de] Sent: Friday, December 21, 2012 1:01 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] Rterm does not load personal library On 21.12.2012 18:16, Bond, Stephen wrote: Greetings, I am trying to run a short script from a shell: c:\projects\hellRscript --default-packages=mypack X:/4Stephen/commit/curve.R X:/4Stephen/commit/run1.out Loading required package: utils Warning message: package 'RODBC' was built under R version 2.12.2 Error: could not find function normalize Execution halted Warning message: closing unused RODBC handle 1 function 'normalize' is in mypack and I even attempt to load it from inside the script: ## curve.R dt - Sys.Date() library(RODBC) pipe - odbcConnect(commit) library(mypack) # with quotes it's the same result normalize(pipe) any ideas how to make Rterm load the library?? Running from inside ESS works fine. Also from RGui. It does not work if I start a session by R -vanilla on the command line. I guess some dependencies on other base R packages are not declared and those packages are not loaded by your package. If you start RGui some packages are loaded that are not loaded if you start via --vanilla, but we cannot know more without information about mypack. Uwe Ligges Thanks and happy holidays. Stephen __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] rms R code
Greetings, useRs. Does anybody have replication of the examples from the RMS book by Harrell coded in R? I find that most the code does not work and it takes too much time to debug. For example from p.276 age.t - w[,age] f.full - lrm(cvd~scored(rx)+rcs(dtime,5)+age.t+wt.t+pf.t+hx+sbp+ekg.t+sz.t+sg.t+ap.t+bm+hg.t,x=T,y=T) Error in model.frame.default(formula = cvd ~ scored(rx) + rcs(dtime, 5) + : invalid type (list) for variable 'age.t' thanks everybody. Stephen __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] rms R code
I have used the errata from there, but have not found where to download working R code from. Stephen -Original Message- From: Sarah Goslee [mailto:sarah.gos...@gmail.com] Sent: Monday, December 17, 2012 4:22 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] rms R code Did you try the website for the book? http://biostat.mc.vanderbilt.edu/wiki/Main/RmS On Mon, Dec 17, 2012 at 4:06 PM, Bond, Stephen stephen.b...@cibc.com wrote: Greetings, useRs. Does anybody have replication of the examples from the RMS book by Harrell coded in R? I find that most the code does not work and it takes too much time to debug. For example from p.276 age.t - w[,age] f.full - lrm(cvd~scored(rx)+rcs(dtime,5)+age.t+wt.t+pf.t+hx+sbp+ekg.t+sz.t+sg.t+ap.t+bm+hg.t,x=T,y=T) Error in model.frame.default(formula = cvd ~ scored(rx) + rcs(dtime, 5) + : invalid type (list) for variable 'age.t' thanks everybody. Stephen -- Sarah Goslee http://www.functionaldiversity.org __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] RODBC to MS SQL Server update error
Is this a bug: Trying to update when the where condition gives zero rows throws an error on MS SQL server sqlQuery(pipe,select * from ComDetailCurrent where RateTypeId is null;) [1] ProcessDate SourceSystemIdAccountNumber Xref1 0 rows (or 0-length row.names) sqlQuery(pipe,update ComDetailCurrent set RateTypeId=1 where RateTypeId is null;) [1] [RODBC] ERROR: Could not SQLExecDirect 'update ComDetailCurrent set RateTypeId=1 where RateTypeId is null;' On the other hand with MySQL sqlQuery(mys,update rate_hist set val=5.0 where dtend'1986-05-01') [1] No Data Please, try to update a table on SQL server and give it a where condition such that zero rows will be updated. And post the result on this forum :-) Thanks everybody and have a great day. Stephen B __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Hmisc describe error
Describe fails for me with a message similar to what was an issue in 2008 and
got fixed according to posts.
R version 2.15.0 (2012-03-30)
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i386-pc-mingw32/i386 (32-bit)
# output truncated
options(chmhelp = FALSE, help_type = text)
.help.ESS - help
options(STERM='iESS', editor='gnuclient.exe')
load('prostate.sav')
library(rms)
Loading required package: Hmisc
Loading required package: survival
Loading required package: splines
Hmisc library by Frank E Harrell Jr
Type library(help='Hmisc'), ?Overview, or ?Hmisc.Overview')
to see overall documentation.
NOTE:Hmisc no longer redefines [.factor to drop unused levels when
subsetting. To get the old behavior of Hmisc type dropUnusedLevels().
Attaching package: 'Hmisc'
The following object(s) are masked from 'package:survival':
untangle.specials
The following object(s) are masked from 'package:base':
format.pval, round.POSIXt, trunc.POSIXt, units
Attaching package: 'rms'
The following object(s) are masked from 'package:survival':
Surv
describe(prostate)
Error in format(dates(x)) : could not find function dates
Thanks everybody.
Stephen
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Hmisc describe error
Just as a clarification: I downloaded 'prostate.sav' from F. Harrell's website.
For some reason
data(prostate)
Warning message:
In data(prostate) : data set 'prostate' not found
I don't have the prostate data set as is.
-Original Message-
From: David Winsemius [mailto:dwinsem...@comcast.net]
Sent: Monday, October 01, 2012 12:53 PM
To: Bond, Stephen
Cc: r-help@r-project.org
Subject: Re: [R] Hmisc describe error
On Oct 1, 2012, at 9:33 AM, Bond, Stephen wrote:
Describe fails for me with a message similar to what was an issue in 2008 and
got fixed according to posts.
R version 2.15.0 (2012-03-30)
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i386-pc-mingw32/i386 (32-bit)
# output truncated
options(chmhelp = FALSE, help_type = text)
.help.ESS - help
options(STERM='iESS', editor='gnuclient.exe')
load('prostate.sav')
library(rms)
Loading required package: Hmisc
Loading required package: survival
Loading required package: splines
Hmisc library by Frank E Harrell Jr
Type library(help='Hmisc'), ?Overview, or ?Hmisc.Overview')
to see overall documentation.
NOTE:Hmisc no longer redefines [.factor to drop unused levels when
subsetting. To get the old behavior of Hmisc type dropUnusedLevels().
Attaching package: 'Hmisc'
The following object(s) are masked from 'package:survival':
untangle.specials
The following object(s) are masked from 'package:base':
format.pval, round.POSIXt, trunc.POSIXt, units
Attaching package: 'rms'
The following object(s) are masked from 'package:survival':
Surv
describe(prostate)
Error in format(dates(x)) : could not find function dates
Iget anerror when I use that code.
Error in describe(prostate) : object 'prostate' not found
And when I attempt loading it get:
data(prostate)
So you must have a different set of objects or packages loaded on your
installation than I do.
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] odfWeave fails to load
Uwe, Just a user's perspective: there are too many packages that work only on the maintainer's box and it would benefit the community if there were stricter standards for allowing people to post a package. Open systems like Ubuntu have a ratings sytem that allows users to review packages, so the few bad apples are properly labelled and can be avoided by the community. Kind regards Stephen B -Original Message- From: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de] Sent: Wednesday, May 30, 2012 4:23 AM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] odfWeave fails to load See http://cran.r-project.org/web/checks/check_results_odfWeave.html which indicates the package has some problems. Hence CRAN does not make binaries available. Please contact the maintainer. Best, Uwe Ligges On 29.05.2012 16:23, stephenb wrote: R version 2.15.0 (2012-03-30) Copyright (C) 2012 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i386-pc-mingw32/i386 (32-bit) package 'survey' successfully unpacked and MD5 sums checked package 'odfWeave.survey' successfully unpacked and MD5 sums checked library(odfWeave.survey) Loading required package: odfWeave Error: package 'odfWeave' could not be loaded any ideas, anybody?? I had odfSweave on 2.12, but no such thing on 2.15 just odfWeave.survey and it won't load. -- View this message in context: http://r.789695.n4.nabble.com/odfWeave-fails-to-load-tp4631700.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] ff problems
this works: bigcl - read.table(file=bigCL1.csv,sep=',',header=T , colClasses=c(rep(factor,3),numeric,NULL,integer, numeric,integer,NULL,numeric,NULL),nrow=1000) this doesn't: bigcl - read.table.ffdf(file=bigCL1.csv,sep=',',header=T , colClasses=c(TERMGROUP=factor,SN=factor, INS=factor,INCENTIVE=numeric,SRC=character,FLAG=integer, RELAGE=numeric,BURNSUM=integer,prod=character,BAL=numeric,origbal=NULL) ) Error in ff(initdata = initdata, length = length, levels = levels, ordered = ordered, : vmode 'character' not implemented neither this: bigcl - read.table.ffdf(file=bigCL1.csv,sep=',',header=T , + colClasses=c(TERMGROUP=factor,SN=factor, + INS=factor,INCENTIVE=numeric,SRC=NULL,FLAG=integer, + RELAGE=numeric,BURNSUM=integer,prod=NULL,BAL=numeric,origbal=NULL) + ) Error in repnam(colClasses, colnames(x), default = NA) : the following argument names do not match'SRC','prod','origbal' if I load as is without classes then: bigcl - read.table.ffdf(file=bigCL1.csv,sep=',',header=T) bigcl$term - factor(bigcl$term) Error in sort.list(y) : 'x' must be atomic for 'sort.list' Have you called 'sort' on a list? attempting to save it: ffsave(bigcl,file=fcl) Error in system(cmd, input = filelist, intern = TRUE) : 'zip' not found attempting to get rid of the unwanted columns: bigcl - bigcl[,c(-5,-9,-11)] Error: cannot allocate vector of size 129.8 Mb In addition: Warning messages: 1: In class(df) - data.frame : Reached total allocation of 1535Mb: see help(memory.size) 2: In class(df) - data.frame : Reached total allocation of 1535Mb: see help(memory.size) 3: In class(df) - data.frame : Reached total allocation of 1535Mb: see help(memory.size) 4: In class(df) - data.frame : Reached total allocation of 1535Mb: see help(memory.size) ff was supposed to be gentle on RAM usage?? I am attaching the first 1000 lines from the csv file. thanks a lot. Stephen __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] ff usage for glm
Thomas, I tried biglm and it does not work see http://r.789695.n4.nabble.com/unable-to-get-bigglm-working-ATTN-Thomas-Lumley-td2276524.html#a2278381 . There are other posts from people who cannot get biglm working and others who get strange results. Please, advise if you can help. I have row based native code which works, but it is inconvenient as it does not produce an R object, which can be fed to anova etc. offered it to the developer forum, but message is still waiting for mod approval. regards Stephen B -Original Message- From: Thomas Lumley [mailto:tlum...@uw.edu] Sent: Friday, March 30, 2012 7:32 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] ff usage for glm On Sat, Mar 31, 2012 at 9:05 AM, Bond, Stephen stephen.b...@cibc.com wrote: Greetings useRs, Can anyone provide an example how to use ff to feed a very large data frame to glm? The data.frame cannot be loaded in R using conventional read.csv as it is too big. glm(...,data=ff.file) ?? I shouldn't think glm() will work on data that are too big to read into R. However, bigglm() from the biglm package should work. You just need to write a function that supplies chunks of data from ff.file as requested (see the example on ?bigglm). I haven't used ff, but it looks from the documentation as though chunk() will do all the difficult parts. -thomas -- Thomas Lumley Professor of Biostatistics University of Auckland __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] ff usage for glm
Greetings useRs, Can anyone provide an example how to use ff to feed a very large data frame to glm? The data.frame cannot be loaded in R using conventional read.csv as it is too big. glm(...,data=ff.file) ?? Thank you Stephen B __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Weighted mad
Greetings UseRs, Pls advise if there is a way to write a func that can be supplied to aggregate to compute weighted MeanAbsolute Dev (MAD). I am having trouble passing the correct weights from each group level and cannot see the code behind aggregate. But maybe 'aggregate' is not the best way to do that. m1 - aggregate(pool[,c(SMM)],by=list(time=pool$ym),weighted.mean,w=pool$wght) Error in weighted.mean.default(X[[1L]], ...) : 'x' and 'w' must have the same length Apparently the grouping does not work on the additional argument. I am using weighted mean here just to be explicit and avoid supplying a custom function gor weighted MAD, which is not difficult to write by itself. It's making it work with aggreagte that is the problem. aggregate function (x, ...) UseMethod(aggregate) environment: namespace:stats Does not show anything... Stephen B __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] glmD error
I am attemting to use glmD in order to use latex(obj) afterwards.
Glm works fine, glmD throws an error:
conv.go -
glmD(cbind(go.cnt,tot.cnt-go.cnt)~sn+rcs(relAge,4)+termfac+rate:termfac+
+
I(relAge*term(term-1.25))+I((prevbal/tot.cnt)*1e-4)+prod+newmort+
+
I(dnorm(relAge*term-round(relAge*term/12)*12,sd=1.2)*as.integer((relAge*term)3))+
+ I(term %in% c(36,60,72,120))
+ ,data=conv,
+ family=quasibinomial)
Error in if (!length(fname) || !any(fname == zname)) { :
missing value where TRUE/FALSE needed
conv.go -
glm(cbind(go.cnt,tot.cnt-go.cnt)~sn+rcs(relAge,4)+termfac+rate:termfac+
+
I(relAge*term(term-1.25))+I((prevbal/tot.cnt)*1e-4)+prod+newmort+
+
I(dnorm(relAge*term-round(relAge*term/12)*12,sd=1.2)*as.integer((relAge*term)3))+
+ I(term %in% c(36,60,72,120))
+ ,data=conv,
+ family=quasibinomial)
latex(conv.go)
Error in format == :
comparison (1) is possible only for atomic and list types
str(conv)
'data.frame': 25537 obs. of 18 variables:
$ prod : Factor w/ 4 levels CD,CL,CO,..: 2 2 2 1 2 2 2 4 2 2 ...
$ strt :Class 'Date' num [1:25537] 10623 11354 11382 12084 12084 ...
$ matym :Class 'Date' num [1:25537] 13180 13180 13180 13180 13180 ...
$ ins: Factor w/ 2 levels 0,1: 1 1 1 1 1 1 1 1 1 1 ...
$ ym :Class 'Date' num [1:25537] 13149 13149 13149 13149 13149 ...
$ tot.cnt: int 24 117 25 17 25 26 102 95 52 52 ...
$ prevbal: num 1435380 8991665 1652058 1056949 2240357 ...
$ RAM: num 139 160 150 107 160 ...
$ age: num 83.4 59.4 58.9 35.3 35.3 ...
$ rate : num 6.46 7.4 7.6 6 5.44 ...
$ go : num 0 0 0 0 0 0 0 0 0 0 ...
$ go.cnt : num 0 0 0 0 0 0 0 0 0 0 ...
$ SMM.go : num 0 0 0 0 0 0 0 0 0 0 ...
$ term : num 84 60 59 36 36 24 12 6 84 84 ...
$ relAge : num 0.993 0.99 0.999 0.981 0.981 ...
$ sn : Factor w/ 12 levels Jan,Feb,Mar,..: 1 1 1 1 1 1 1 1 1 2 ...
$ termfac: Factor w/ 7 levels 1,2,3,4,..: 6 5 5 3 3 2 1 1 6 6 ...
$ newmort: int 0 0 0 0 0 0 0 0 0 0 ...
Thanks everybody.
Stephen B
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Create design matrix
Greetings useRs, What is the easiest way to create a design matrix of several factor variables? Function gendata in Design seems to do that for a fitted model, but how to do that only on several factor vectors?? The result should be a df with one row for each distinct combination of levels of factors eg for (M,F) (Y,O) We get M Y M O F Y F O In reality I will have more than 1000 rows so doing by hand not good. Maybe there is a way with outer, but I couldn't see it. All the best to everybody. Stephen __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 install a locally built package
Is there a way to put all R code in a single file?
I have too many small files now, and it is inconvenient to edit (I still have
to put everything in one buffer) and when I add just one new func I have to go
through the process of manually editing all help files one by one. When I put
all the code in one file only the first func is loaded.
Thanks everybody
Stephen B
-Original Message-
From: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de]
Sent: Friday, September 16, 2011 3:43 AM
To: Bond, Stephen
Cc: r-help@r-project.org
Subject: Re: [R] how to install a locally built package
On 15.09.2011 20:38, Bond, Stephen wrote:
I got it working by typing a string in getdfv.Rd after \title{
\title{ getdf
%% ~~function to do ... ~~
}
Strange why the skeleton would not do that given it did \name{getdfv}
\alias{getdfv}
One reason is that we try to force people to deal with the documentation which
we feel is very important, even if the stuff is just used by yourself.
Best,
Uwe Ligges
Anyway, I'm happy now.
Stephen Bond
-Original Message-
From: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de]
Sent: Thursday, September 15, 2011 12:15 PM
To: Bond, Stephen
Cc: r-help@r-project.org
Subject: Re: [R] how to install a locally built package
On 15.09.2011 17:47, Bond, Stephen wrote:
Uwe,
That gave me the same error like CMD install
install.packages(C:/Temp/mypack_1.0.tar.gz, repos=NULL,
type=source)
install.packages(C:/Temp/mypack_1.0.tar.gz, repos=NULL,
type=source) Loading required package: stats Loading required
package: utils Loading required package: graphics Loading required
package: splines
* installing *source* package 'mypack' ...
** R
** preparing package for lazy loading
** help
Warning:
C:/DOCUME~1/BondStep/LOCALS~1/Temp/Rtmpw5N7dm/R.INSTALL6ccc14e6/mypac
k/man/mypack-package.Rd:33: All text must be in a section
*** installing help indices
Error in Rd_info(db[[i]]) : Rd files must have a non-empty \title.
See chapter 'Writing R documentation' in manual 'Writing R Extensions'.
* removing 'c:/PROGRA~1/R/R-212~1.1/library/mypack'
* restoring previous 'c:/PROGRA~1/R/R-212~1.1/library/mypack'
Warning message:
In install.packages(C:/Temp/mypack_1.0.tar.gz, repos = NULL, type =
source) :
installation of package 'C:/Temp/mypack_1.0.tar.gz' had non-zero
exit status
is there a way to skip the Rd part? This is for private use only and there
is no help or data files.
Then you have to delete the Rd file that has been generated by
package.skeleton.
Please read the manuals Writing R Extensions and R Installation and
Administration.
Best,
Uwe Ligges
Thank you.
Stephen Bond
-Original Message-
From: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de]
Sent: Thursday, September 15, 2011 11:43 AM
To: Bond, Stephen
Cc: r-help@r-project.org
Subject: Re: [R] how to install a locally built package
On 15.09.2011 17:32, Bond, Stephen wrote:
Hello useRs,
I am trying to put my funcs in a package to avoid clutter in my workspace
as the funcs are now in .Rprofile.
All plain R code no compilations. Using win XP with a full cygwin
install and R2.12 I first did
package.skeleton(mypack,list=getdfv, namespace=T) # just a
single func
which created a folder with the required stuff in it.
Calling R CMD build creates a tar.gz file with warnings about DOS paths.
Trying to install from the R GUI complains
install.packages(mypack,repos=C:/Temp/mypack_1.0.tar.gz)
Warning: unable to access index for repository
C:/Temp/mypack_1.0.tar.gz/bin/windows/contrib/2.12
Warning message:
In getDependencies(pkgs, dependencies, available, lib) :
package 'mypack' is not available
install.packages(C:/Temp/mypack_1.0.tar.gz, repos=NULL,
type=source)
seems to be what you want.
Uwe Ligges
putting the mypack folder in zip and then trying to install from the
zip produces no error message but then
library(mypack)
Error in library(mypack) : 'mypack' is not a valid installed
package
Just putting the folder in the library folder of R did not work
either
Finally
R CMD Install complains about missing title in the Rd file. I do not have a
Rd file as the skeleton did not create it.
Please, advise how to make simple R code available to be called
without showing in ls()
Thank you.
Stephen Bond
[[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] how to install a locally built package
Uwe,
Are u saying
1) I can add the new func in one of the existing .R files in the R dir??
Or
2) add a new .R file in the same dir and ignore the lack of a matching .Rd
file?
Thank you.
Stephen Bond | Senior Analyst | Treasury Analytics | 416-956-3092
-Original Message-
From: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de]
Sent: Thursday, October 06, 2011 9:12 AM
To: Bond, Stephen
Cc: r-help@r-project.org
Subject: Re: [R] how to install a locally built package
On 06.10.2011 15:10, Bond, Stephen wrote:
Is there a way to put all R code in a single file?
I have too many small files now, and it is inconvenient to edit (I still have
to put everything in one buffer) and when I add just one new func I have to
go through the process of manually editing all help files one by one. When I
put all the code in one file only the first func is loaded.
He? You really do not need to edit anything when you add a function.
Just don't use package skeleton but just add the function.
Uwe Ligges
Thanks everybody
Stephen B
-Original Message-
From: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de]
Sent: Friday, September 16, 2011 3:43 AM
To: Bond, Stephen
Cc: r-help@r-project.org
Subject: Re: [R] how to install a locally built package
On 15.09.2011 20:38, Bond, Stephen wrote:
I got it working by typing a string in getdfv.Rd after \title{
\title{ getdf
%% ~~function to do ... ~~
}
Strange why the skeleton would not do that given it did \name{getdfv}
\alias{getdfv}
One reason is that we try to force people to deal with the documentation
which we feel is very important, even if the stuff is just used by yourself.
Best,
Uwe Ligges
Anyway, I'm happy now.
Stephen Bond
-Original Message-
From: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de]
Sent: Thursday, September 15, 2011 12:15 PM
To: Bond, Stephen
Cc: r-help@r-project.org
Subject: Re: [R] how to install a locally built package
On 15.09.2011 17:47, Bond, Stephen wrote:
Uwe,
That gave me the same error like CMD install
install.packages(C:/Temp/mypack_1.0.tar.gz, repos=NULL,
type=source)
install.packages(C:/Temp/mypack_1.0.tar.gz, repos=NULL,
type=source) Loading required package: stats Loading required
package: utils Loading required package: graphics Loading required
package: splines
* installing *source* package 'mypack' ...
** R
** preparing package for lazy loading
** help
Warning:
C:/DOCUME~1/BondStep/LOCALS~1/Temp/Rtmpw5N7dm/R.INSTALL6ccc14e6/mypa
c
k/man/mypack-package.Rd:33: All text must be in a section
*** installing help indices
Error in Rd_info(db[[i]]) : Rd files must have a non-empty \title.
See chapter 'Writing R documentation' in manual 'Writing R Extensions'.
* removing 'c:/PROGRA~1/R/R-212~1.1/library/mypack'
* restoring previous 'c:/PROGRA~1/R/R-212~1.1/library/mypack'
Warning message:
In install.packages(C:/Temp/mypack_1.0.tar.gz, repos = NULL, type =
source) :
installation of package 'C:/Temp/mypack_1.0.tar.gz' had
non-zero exit status
is there a way to skip the Rd part? This is for private use only and there
is no help or data files.
Then you have to delete the Rd file that has been generated by
package.skeleton.
Please read the manuals Writing R Extensions and R Installation
and Administration.
Best,
Uwe Ligges
Thank you.
Stephen Bond
-Original Message-
From: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de]
Sent: Thursday, September 15, 2011 11:43 AM
To: Bond, Stephen
Cc: r-help@r-project.org
Subject: Re: [R] how to install a locally built package
On 15.09.2011 17:32, Bond, Stephen wrote:
Hello useRs,
I am trying to put my funcs in a package to avoid clutter in my workspace
as the funcs are now in .Rprofile.
All plain R code no compilations. Using win XP with a full cygwin
install and R2.12 I first did
package.skeleton(mypack,list=getdfv, namespace=T) # just a
single func
which created a folder with the required stuff in it.
Calling R CMD build creates a tar.gz file with warnings about DOS paths.
Trying to install from the R GUI complains
install.packages(mypack,repos=C:/Temp/mypack_1.0.tar.gz)
Warning: unable to access index for repository
C:/Temp/mypack_1.0.tar.gz/bin/windows/contrib/2.12
Warning message:
In getDependencies(pkgs, dependencies, available, lib) :
package 'mypack' is not available
install.packages(C:/Temp/mypack_1.0.tar.gz, repos=NULL,
type=source)
seems to be what you want.
Uwe Ligges
putting the mypack folder in zip and then trying to install from
the zip produces no error message but then
library(mypack)
Error in library(mypack) : 'mypack' is not a valid installed
package
Just putting the folder in the library folder of R did not work
either
Finally
R CMD Install complains about missing title in the Rd file. I do not have
a Rd file as the skeleton did not create it.
Please, advise how to make simple R
[R] how to install a locally built package
Hello useRs, I am trying to put my funcs in a package to avoid clutter in my workspace as the funcs are now in .Rprofile. All plain R code no compilations. Using win XP with a full cygwin install and R2.12 I first did package.skeleton(mypack,list=getdfv, namespace=T) # just a single func which created a folder with the required stuff in it. Calling R CMD build creates a tar.gz file with warnings about DOS paths. Trying to install from the R GUI complains install.packages(mypack,repos=C:/Temp/mypack_1.0.tar.gz) Warning: unable to access index for repository C:/Temp/mypack_1.0.tar.gz/bin/windows/contrib/2.12 Warning message: In getDependencies(pkgs, dependencies, available, lib) : package 'mypack' is not available putting the mypack folder in zip and then trying to install from the zip produces no error message but then library(mypack) Error in library(mypack) : 'mypack' is not a valid installed package Just putting the folder in the library folder of R did not work either Finally R CMD Install complains about missing title in the Rd file. I do not have a Rd file as the skeleton did not create it. Please, advise how to make simple R code available to be called without showing in ls() Thank you. Stephen Bond [[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] how to install a locally built package
Uwe, That gave me the same error like CMD install install.packages(C:/Temp/mypack_1.0.tar.gz, repos=NULL, type=source) install.packages(C:/Temp/mypack_1.0.tar.gz, repos=NULL, type=source) Loading required package: stats Loading required package: utils Loading required package: graphics Loading required package: splines * installing *source* package 'mypack' ... ** R ** preparing package for lazy loading ** help Warning: C:/DOCUME~1/BondStep/LOCALS~1/Temp/Rtmpw5N7dm/R.INSTALL6ccc14e6/mypack/man/mypack-package.Rd:33: All text must be in a section *** installing help indices Error in Rd_info(db[[i]]) : Rd files must have a non-empty \title. See chapter 'Writing R documentation' in manual 'Writing R Extensions'. * removing 'c:/PROGRA~1/R/R-212~1.1/library/mypack' * restoring previous 'c:/PROGRA~1/R/R-212~1.1/library/mypack' Warning message: In install.packages(C:/Temp/mypack_1.0.tar.gz, repos = NULL, type = source) : installation of package 'C:/Temp/mypack_1.0.tar.gz' had non-zero exit status is there a way to skip the Rd part? This is for private use only and there is no help or data files. Thank you. Stephen Bond -Original Message- From: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de] Sent: Thursday, September 15, 2011 11:43 AM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] how to install a locally built package On 15.09.2011 17:32, Bond, Stephen wrote: Hello useRs, I am trying to put my funcs in a package to avoid clutter in my workspace as the funcs are now in .Rprofile. All plain R code no compilations. Using win XP with a full cygwin install and R2.12 I first did package.skeleton(mypack,list=getdfv, namespace=T) # just a single func which created a folder with the required stuff in it. Calling R CMD build creates a tar.gz file with warnings about DOS paths. Trying to install from the R GUI complains install.packages(mypack,repos=C:/Temp/mypack_1.0.tar.gz) Warning: unable to access index for repository C:/Temp/mypack_1.0.tar.gz/bin/windows/contrib/2.12 Warning message: In getDependencies(pkgs, dependencies, available, lib) : package 'mypack' is not available install.packages(C:/Temp/mypack_1.0.tar.gz, repos=NULL, type=source) seems to be what you want. Uwe Ligges putting the mypack folder in zip and then trying to install from the zip produces no error message but then library(mypack) Error in library(mypack) : 'mypack' is not a valid installed package Just putting the folder in the library folder of R did not work either Finally R CMD Install complains about missing title in the Rd file. I do not have a Rd file as the skeleton did not create it. Please, advise how to make simple R code available to be called without showing in ls() Thank you. Stephen Bond [[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] how to install a locally built package
I got it working by typing a string in getdfv.Rd after \title{
\title{ getdf
%% ~~function to do ... ~~
}
Strange why the skeleton would not do that given it did
\name{getdfv}
\alias{getdfv}
Anyway, I'm happy now.
Stephen Bond
-Original Message-
From: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de]
Sent: Thursday, September 15, 2011 12:15 PM
To: Bond, Stephen
Cc: r-help@r-project.org
Subject: Re: [R] how to install a locally built package
On 15.09.2011 17:47, Bond, Stephen wrote:
Uwe,
That gave me the same error like CMD install
install.packages(C:/Temp/mypack_1.0.tar.gz, repos=NULL, type=source)
install.packages(C:/Temp/mypack_1.0.tar.gz, repos=NULL, type=source)
Loading required package: stats
Loading required package: utils
Loading required package: graphics
Loading required package: splines
* installing *source* package 'mypack' ...
** R
** preparing package for lazy loading
** help
Warning:
C:/DOCUME~1/BondStep/LOCALS~1/Temp/Rtmpw5N7dm/R.INSTALL6ccc14e6/mypack/man/mypack-package.Rd:33:
All text must be in a section
*** installing help indices
Error in Rd_info(db[[i]]) : Rd files must have a non-empty \title.
See chapter 'Writing R documentation' in manual 'Writing R Extensions'.
* removing 'c:/PROGRA~1/R/R-212~1.1/library/mypack'
* restoring previous 'c:/PROGRA~1/R/R-212~1.1/library/mypack'
Warning message:
In install.packages(C:/Temp/mypack_1.0.tar.gz, repos = NULL, type =
source) :
installation of package 'C:/Temp/mypack_1.0.tar.gz' had non-zero exit
status
is there a way to skip the Rd part? This is for private use only and there is
no help or data files.
Then you have to delete the Rd file that has been generated by
package.skeleton.
Please read the manuals Writing R Extensions and R Installation and
Administration.
Best,
Uwe Ligges
Thank you.
Stephen Bond
-Original Message-
From: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de]
Sent: Thursday, September 15, 2011 11:43 AM
To: Bond, Stephen
Cc: r-help@r-project.org
Subject: Re: [R] how to install a locally built package
On 15.09.2011 17:32, Bond, Stephen wrote:
Hello useRs,
I am trying to put my funcs in a package to avoid clutter in my workspace as
the funcs are now in .Rprofile.
All plain R code no compilations. Using win XP with a full cygwin install
and R2.12
I first did
package.skeleton(mypack,list=getdfv, namespace=T) # just a single func
which created a folder with the required stuff in it.
Calling R CMD build creates a tar.gz file with warnings about DOS paths.
Trying to install from the R GUI complains
install.packages(mypack,repos=C:/Temp/mypack_1.0.tar.gz)
Warning: unable to access index for repository
C:/Temp/mypack_1.0.tar.gz/bin/windows/contrib/2.12
Warning message:
In getDependencies(pkgs, dependencies, available, lib) :
package 'mypack' is not available
install.packages(C:/Temp/mypack_1.0.tar.gz, repos=NULL, type=source)
seems to be what you want.
Uwe Ligges
putting the mypack folder in zip and then trying to install from the zip
produces no error message but then
library(mypack)
Error in library(mypack) : 'mypack' is not a valid installed package
Just putting the folder in the library folder of R did not work either
Finally
R CMD Install complains about missing title in the Rd file. I do not have a
Rd file as the skeleton did not create it.
Please, advise how to make simple R code available to be called without
showing in ls()
Thank you.
Stephen Bond
[[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] model formula
Hello useRs, Pls help with removing a single interaction term from a formula: summary( glm.turn.2 - glm(cbind(turn.cnt,tot.cnt-turn.cnt)~sn+poly(relAge,2,raw=T)+termfac+rate:termfac,data=fix, family=quasibinomial) ) Gives Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) -7.028467 0.106002 -66.305 2e-16 *** snFeb 0.156963 0.023660 6.634 3.27e-11 *** snMar 0.317540 0.022883 13.876 2e-16 *** snApr 0.526084 0.022004 23.908 2e-16 *** snMay 1.026710 0.020347 50.460 2e-16 *** snJun 0.841044 0.021318 39.452 2e-16 *** snJul 0.668790 0.022530 29.685 2e-16 *** snAug 0.544267 0.022580 24.104 2e-16 *** snSep 0.389667 0.023363 16.679 2e-16 *** snOct 0.351294 0.023586 14.894 2e-16 *** snNov 0.391464 0.024057 16.272 2e-16 *** snDec -0.373369 0.028755 -12.985 2e-16 *** poly(relAge, 2, raw = T)1 2.952887 0.067455 43.776 2e-16 *** poly(relAge, 2, raw = T)2 -1.783907 0.064074 -27.841 2e-16 *** termfac2 -0.681719 0.128571 -5.302 1.14e-07 *** termfac3 -1.032416 0.146396 -7.052 1.77e-12 *** termfac4 -1.267011 0.108940 -11.630 2e-16 *** termfac5 -1.009922 0.213129 -4.739 2.15e-06 *** termfac6 3.300301 0.203465 16.221 2e-16 *** termfac1:rate 0.012274 0.019895 0.6170.537 termfac2:rate 0.121724 0.013472 9.036 2e-16 *** termfac3:rate 0.175232 0.018987 9.229 2e-16 *** termfac4:rate 0.197726 0.005787 34.164 2e-16 *** termfac5:rate 0.145622 0.027295 5.335 9.56e-08 *** termfac6:rate -0.362379 0.025261 -14.345 2e-16 *** and I would like to remove the termfac1:rate interaction. Is there a way to do that? Thank you Stephen Bond [[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] model formula
Using glm.fit I can remove the unwanted column from model.matrix(glm.turn.2) manually and then call glm.fit directly. Q: how can I stuff the output from glm.fit into a regular glm object, so I can use prediction functions, etc.?? the help does not describe how to use output from glm.fit in the way a glm model can be used. Thanks everybody. Stephen -Original Message- From: Uwe Ligges [mailto:lig...@statistik.tu-dortmund.de] Sent: Thursday, August 11, 2011 11:40 AM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] model formula On 11.08.2011 17:27, Bond, Stephen wrote: Hello useRs, Pls help with removing a single interaction term from a formula: summary( glm.turn.2- glm(cbind(turn.cnt,tot.cnt-turn.cnt)~sn+poly(relAge,2,raw=T)+termfac+rate:termfac,data=fix, family=quasibinomial) ) Gives Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) -7.028467 0.106002 -66.305 2e-16 *** snFeb 0.156963 0.023660 6.634 3.27e-11 *** snMar 0.317540 0.022883 13.876 2e-16 *** snApr 0.526084 0.022004 23.908 2e-16 *** snMay 1.026710 0.020347 50.460 2e-16 *** snJun 0.841044 0.021318 39.452 2e-16 *** snJul 0.668790 0.022530 29.685 2e-16 *** snAug 0.544267 0.022580 24.104 2e-16 *** snSep 0.389667 0.023363 16.679 2e-16 *** snOct 0.351294 0.023586 14.894 2e-16 *** snNov 0.391464 0.024057 16.272 2e-16 *** snDec -0.373369 0.028755 -12.985 2e-16 *** poly(relAge, 2, raw = T)1 2.952887 0.067455 43.776 2e-16 *** poly(relAge, 2, raw = T)2 -1.783907 0.064074 -27.841 2e-16 *** termfac2 -0.681719 0.128571 -5.302 1.14e-07 *** termfac3 -1.032416 0.146396 -7.052 1.77e-12 *** termfac4 -1.267011 0.108940 -11.630 2e-16 *** termfac5 -1.009922 0.213129 -4.739 2.15e-06 *** termfac6 3.300301 0.203465 16.221 2e-16 *** termfac1:rate 0.012274 0.019895 0.6170.537 termfac2:rate 0.121724 0.013472 9.036 2e-16 *** termfac3:rate 0.175232 0.018987 9.229 2e-16 *** termfac4:rate 0.197726 0.005787 34.164 2e-16 *** termfac5:rate 0.145622 0.027295 5.335 9.56e-08 *** termfac6:rate -0.362379 0.025261 -14.345 2e-16 *** and I would like to remove the termfac1:rate interaction. Is there a way to do that? Thank you And what do you want to do with the observations with termfac=1 then? Ignore? Also move to the Intercept? Uwe Ligges Stephen Bond [[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] coxph fails to survfit
Responding to T. Therneau: sessionInfo() R version 2.10.1 (2009-12-14) i386-pc-mingw32 locale: [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 LC_MONETARY=English_United States.1252 [4] LC_NUMERIC=C LC_TIME=English_United States.1252 attached base packages: [1] splines stats graphics grDevices utils datasets methods base other attached packages: [1] survival_2.35-8 Also: s1 - predict(mod1,newdata=inc[50050:50100,], + type=expected) Error in predict.coxph(mod1, newdata = inc[50050:50100, ], type = expected) : Method not yet finished Stephen -Original Message- From: Terry Therneau [mailto:thern...@mayo.edu] Sent: Friday, February 04, 2011 8:54 AM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: coxph fails to survfit I have a model with quant vars only and the error message does not make sense: Could you tell us what version of S and of the survival package you are using? You can get this with sessionInfo(), see the posting guide for details. This would help me identify the issue. I was planning to post an update to CRAN this weekend that has changes to survfit.coxph, addressing a different issue than yours -- models with a strata by covariate interaction --- but this would be a very good time to address your issue too if it is a code bug. Terry Therneau (author of survival). __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] coxph fails to survfit
I have a model with quant vars only and the error message does not make sense: (mod1 - coxph(Surv(time=strt,time2=stp,event=(resp==1))~ +incpost+I(amt/1e5)+rate+strata(termfac), subset=dt2010-08-30, data=inc,method=efron)) Call: coxph(formula = Surv(time = strt, time2 = stp, event = (resp == 1)) ~ +incpost + I(amt/1e+05) + rate + strata(termfac), data = inc, subset = dt 2010-08-30, method = efron) coef exp(coef) se(coef) z p incpost 0.2563 1.292 0.02479 10.34 0.0e+00 I(amt/1e+05) -0.0532 0.948 0.00487 -10.92 0.0e+00 rate -0.0507 0.951 0.00945 -5.36 8.2e-08 Likelihood ratio test=295 on 3 df, p=0 n= 1192634 length(mod1$xlevels) [1] 0 # now calling survfit with just a few rows from the data s1 - survfit(mod1,newdata=inc[50050:50100,],se.fit=F,individual=T,type=aa) Error in `contrasts-`(`*tmp*`, value = contr.treatment) : contrasts can be applied only to factors with 2 or more levels inc[50050:50100,] dt inc incpost strt stp resp rateamt p1 matfac termfac vol 50050 2009-02-05 0.00.000 10 5.35 266833 C 5 3 14229571 50051 2009-02-06 -0.095750.001 20 5.35 266833 C 5 3 14229571 50052 2009-02-09 -0.038750.004 50 5.35 266833 C 5 3 14229571 50053 2009-02-10 0.008500.005 60 5.35 266833 C 5 3 14229571 50054 2009-02-11 0.002500.006 70 5.35 266833 C 5 3 14229571 50055 2009-02-12 -0.047250.007 80 5.35 266833 C 5 3 14229571 50056 2009-02-13 -0.000750.008 90 5.35 266833 C 5 3 14229571 50057 2009-02-16 -0.104250.00 11 120 5.35 266833 C 5 3 14229571 50058 2009-02-17 -0.104250.00 12 130 5.35 266833 C 5 3 14229571 50059 2009-02-18 -0.029250.00 13 140 5.35 266833 C 5 3 14229571 50060 2009-02-19 -0.013000.00 14 150 5.35 266833 C 5 3 14229571 50061 2009-02-20 -0.040750.00 15 160 5.35 266833 C 5 3 14229571 50062 2009-02-23 -0.031000.00 18 190 5.35 266833 C 5 3 14229571 50063 2009-02-24 0.019750.00 19 200 5.35 266833 C 5 3 14229571 50064 2009-02-25 0.130500.00 20 210 5.35 266833 C 5 3 14229571 50065 2009-02-26 0.127250.00 21 220 5.35 266833 C 5 3 14229571 ... further output truncated Please, comment if you see what's the issue. Thank you. Stephen [[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] coxph fails to survfit
Responding to the suggestion by D Winsemius : s1 - survfit(mod1,newdata=inc[50050:50100,c(strt,stp,incpost, amt,rate, termfac)], + se.fit=F,individual=T,type=aa) Error in Surv(time = strt, time2 = stp, event = (resp == 1)) : object 'resp' not found it appears it wants to fit a survival curve instead of predicting a survival curve for the subject using the explanatory vars. I think I probably need predict.coxph(...,type=expected) or write code for matching time dependent risk to the proper time index of the survival curve. Thanks David. Stephen -Original Message- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Thursday, February 03, 2011 3:10 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] coxph fails to survfit On Feb 3, 2011, at 2:14 PM, Bond, Stephen wrote: I have a model with quant vars only and the error message does not make sense: (mod1 - coxph(Surv(time=strt,time2=stp,event=(resp==1))~ +incpost +I(amt/1e5)+rate+strata(termfac), subset=dt2010-08-30, data=inc,method=efron)) Call: coxph(formula = Surv(time = strt, time2 = stp, event = (resp == 1)) ~ +incpost + I(amt/1e+05) + rate + strata(termfac), data = inc, subset = dt 2010-08-30, method = efron) coef exp(coef) se(coef) z p incpost 0.2563 1.292 0.02479 10.34 0.0e+00 I(amt/1e+05) -0.0532 0.948 0.00487 -10.92 0.0e+00 rate -0.0507 0.951 0.00945 -5.36 8.2e-08 Likelihood ratio test=295 on 3 df, p=0 n= 1192634 length(mod1$xlevels) [1] 0 # now calling survfit with just a few rows from the data s1 - survfit (mod1,newdata=inc[50050:50100,],se.fit=F,individual=T,type=aa) Error in `contrasts-`(`*tmp*`, value = contr.treatment) : contrasts can be applied only to factors with 2 or more levels I obviously cannot test my theory, but would have thought the call would be: s1 - survfit(mod1, newdata=inc[50050:50100,], c(incpost, amt, rate, termfac) ], se.fit=F,individual=T,type=aa) I.e., a data.frame for newdata that matched the variables used in the RHS of the formula. inc[50050:50100,] dt inc incpost strt stp resp rateamt p1 matfac termfac vol 50050 2009-02-05 0.00.000 10 5.35 266833 C 5 3 14229571 50051 2009-02-06 -0.095750.001 20 5.35 266833 C 5 3 14229571 50052 2009-02-09 -0.038750.004 50 5.35 266833 C 5 3 14229571 50053 2009-02-10 0.008500.005 60 5.35 266833 C 5 3 14229571 50054 2009-02-11 0.002500.006 70 5.35 266833 C 5 3 14229571 50055 2009-02-12 -0.047250.007 80 5.35 266833 C 5 3 14229571 50056 2009-02-13 -0.000750.008 90 5.35 266833 C 5 3 14229571 50057 2009-02-16 -0.104250.00 11 120 5.35 266833 C 5 3 14229571 50058 2009-02-17 -0.104250.00 12 130 5.35 266833 C 5 3 14229571 50059 2009-02-18 -0.029250.00 13 140 5.35 266833 C 5 3 14229571 50060 2009-02-19 -0.013000.00 14 150 5.35 266833 C 5 3 14229571 50061 2009-02-20 -0.040750.00 15 160 5.35 266833 C 5 3 14229571 50062 2009-02-23 -0.031000.00 18 190 5.35 266833 C 5 3 14229571 50063 2009-02-24 0.019750.00 19 200 5.35 266833 C 5 3 14229571 50064 2009-02-25 0.130500.00 20 210 5.35 266833 C 5 3 14229571 50065 2009-02-26 0.127250.00 21 220 5.35 266833 C 5 3 14229571 ... further output truncated Please, comment if you see what's the issue. Thank you. Stephen [[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 West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] coxph strange result
The following fit does not make sense to me, please, correct me if I have a logical error. moddowsn Call: coxph(formula = Surv(start, stop, resp) ~ sn + matfac2, data = coxsn1, method = efron) coef exp(coef) se(coef) z p sn2 0.0497 1.051 0.02030 2.450 1.4e-02 sn3 -0.0532 0.948 0.02038 -2.610 9.0e-03 sn4 -0.0410 0.960 0.01979 -2.073 3.8e-02 sn5 -0.0776 0.925 0.01954 -3.973 7.1e-05 sn6 -0.1133 0.893 0.01839 -6.161 7.2e-10 sn7 -0.1252 0.882 0.01846 -6.781 1.2e-11 sn8 -0.1222 0.885 0.01994 -6.130 8.8e-10 sn9 -0.0507 0.951 0.02047 -2.478 1.3e-02 sn10 -0.0444 0.957 0.02056 -2.159 3.1e-02 sn11 -0.0433 0.958 0.02157 -2.008 4.5e-02 sn12 -0.0114 0.989 0.02037 -0.557 5.8e-01 matfac22 -0.2599 0.771 0.01727 -15.048 0.0e+00 matfac25 -0.1804 0.835 0.00924 -19.512 0.0e+00 Likelihood ratio test=651 on 13 df, p=0 n= 253802 This would indicate that in sn6 to sn8 there is less of a chance of an event. ?? do the relative frequencies implied by the following table make any sense?? table(coxsn1$matfac2,coxsn1$sn,coxsn1$resp) , , = 0 1 2 3 4 5 6 7 8 9101112 1 3407 3177 3425 3348 3975 3564 3181 3077 2894 2610 3441 3443 2 920 1005 1142 1327 1645 1530 1330 1184 964 864 888 860 5 9036 9507 10258 11888 16826 15575 13394 12346 9938 9001 8970 8599 , , = 1 1 2 3 4 5 6 7 8 9101112 1 1453 1459 1186 1496 1295 1754 1429 1153 1106 1234 965 1532 2 312 290 330 390 454 539 479 367 295 276 256 267 5 2994 3207 3371 3629 4095 5581 5837 3844 3400 3199 2705 3084 Apparently the frequency of an event is higher in the summer months. I apologize for not being able to disclose the dataset, but think that the table provides enough to address the question. Thanks everybody. Stephen B [[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] coxph linear.predictors
Re: 1. X*beta != linear.predictor.
The equality is stated in three different help docs, which is misleading,
especially in light of the way glm is set up. I felt like was wrestling with
SAS :-)
The relative risk was the original idea behind cox regression, but it can be
used for many non-relative purposes. If we want to calculate death probability
in each period, then lp is no longer shift invariant.
Re: 2. Survfit is too slow.
It seems that the implementation follows the procedure in the original Cox
paper, which calls iterative optimization for each death time.
My subjects are mortgages and both the estimation and the prediction samples
are several hundred thousand. The call appears to recalculate/optimize
everything even though only the $surv changes. Since each subject belongs to a
single strata, most of the calculations are redundant.
I am not much of a programmer and could never figure out how to use the R
profiler, so cannot be exact here, but the simple exponentiation takes no time
and survfit takes several secs for each subject.
So I did:
survlong - survfit(modlong) # a single call suffices
bl1 - c(1,cumsum(survlong$strata)+1)
bl2 - cumsum(survlong$strata) # get the start and end of each strata
for (jj in 1:nrow(newapp)){
strat=as.integer(newapp[jj,termfac])
surv - survlong$surv[(b1[strat]):(b2[strat])] # extract the strata
risk - predict(modlong,new=newapp[jj,],type=risk)# it seems there is no
# optimization here
newsurv - surv^risk # we done
... rest of code
}
As a package maintainer, you have to decide whether including any of the above
and below is useful or users can figure out things on their own. Or maybe
survfit can be made smart and subsequent calls on the same model will use the
first call to survfit?? It's your call :-)
Kind regards
Stephen B
-Original Message-
From: Terry Therneau [mailto:thern...@mayo.edu]
Sent: Thursday, October 28, 2010 6:39 PM
To: Bond, Stephen; David Winsemius
Cc: r-help@r-project.org
Subject: Re: [R] coxph linear.predictors
Gentlemen,
I read R-news in batch mode so I'm often a day behind. Let me try to
answer some of the questions.
1. X*beta != linear.predictor.
I'm sorry if the documentation isn't all it could be. Between the book,
tech report, and help I've written about 400 pages, but this particular
topic isn't yet in it. The final snipe about being opaque like SAS
was really unfair.
The Cox model is a relative risk model, if lp is a linear predictor then
so is lp +c for any constant; they are equally good and equally valid.
The linear.predictor component in a coxph fit is (X-means) * beta. The
computation exp(lp) occurs multiple times downstream and this keep the
exp function from overflowing when there is something like a Date object
as a predictor. Adding this constant changes not a single downstream
calcuation.
2. Survfit is too slow.
I'd like to hear more about this. My work mostly involves modest data
sets so perhaps I haven't seen it. Accuracy and maintainability have
been my first worries.
3. Baseline survival.
Let xbase be a particular set of values for the x covariates (one for
each). The survival curve for a given xbase is obtained from survfit
fit - coxph(
sfit - survfit(fit, newdata=xbase)
chaz - -log(sfit$surv) #cumulative hazard
(The xbase vector will need to have variable names for the function to
know which value goes to which of course).
The cumulative hazard for any other subject will be
newhaz - chaz * exp(fit$coef%*% (x-xbase))
There is not a simple transformation of the standard error from one fit
to another, however. You will need to call survfit with a data frame
for newdata, which will return one curve per row with the proper values.
In my view there is no such thing as A baseline survival curve. Any
xbase you chose is a baseline. However, it is wise to choose something
near the center of the data space in order to avoid numeric problems
with the exp function above. I would never ever chose a vector of
zeros, although some text books do -- it saves them about 8 characters
of typing in the newhaz formula above.
Terry Therneau
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] coxph linear.predictors
To close the issue: S.ave - survfit(fit) # this is average survival S.1 - (S.ave$surv)^exp(fit$linear.predictors[1]) # get subject 1 S.1 [1] 0.848450223861993 0.696973203043377 0.530604282995790 0.391136296063732 0.228370373409774 0.132991279576665 0.071742127376216 [8] 0.030988671457923 0.012340713882587 0.004553665545169 0.000865566322530 0.000135746409544 survfit(fit,new=ovarian[1,])$surv # subject 1 from survfit [1] 0.848450223861993 0.696973203043377 0.530604282995791 0.391136296063732 0.228370373409774 0.132991279576665 0.071742127376216 [8] 0.030988671457923 0.012340713882587 0.004553665545169 0.000865566322530 0.000135746409544 A single call to survfit will suffice for thousands of subjects, avoiding the iterative loop calculations inside survfit. Stephen Bond -Original Message- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Wednesday, October 27, 2010 1:15 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] coxph linear.predictors On Oct 27, 2010, at 12:12 PM, Bond, Stephen wrote: I would like to be able to construct hazard rates (or unconditional death prob) Hazards are not probabilities (since probabilities are constrained to the range [0,1] and hazards are unbounded upward.) for many subjects from a given survfit. This will involve adjusting the ( n.event/n.risk) with (coxph object )$linear.predictors I must be having another silly day as I cannot reproduce the linear predictor: fit - coxph(Surv(futime, fustat) ~ age, data = ovarian) fit$linear.predictors[1] [1] 2.612756 That's the linear predictor (the beta*X) and that particular number only applies to the first case. coef(fit)*model.matrix(fit)[1,1] age 11.69021 I don't know what that might be and you are not telling us what you think it is. The above is based on the help listing for coxph.object coefficients: the coefficients of the linear predictor, which multiply the columns of the model matrix. If the model is over-determined there will be missing values in the vector corresponding to the redundant columns in the model matrix. Also, please comment whether n.event/n.risk The Nelson-Aalen estimator of the cumulative hazard as a function of intervals prior to t is sum( n-event(t)/ n.risk(t)) gives the baseline hazard exp(alpha) ? No. The baseline hazard, as you are calling this, would be an estimate for persons with all covariates = 0, so in this case is for women of age=0. (Not a particularly interpretable result in many situations. The baseline hazard following treated ovarian cancer for neonates is not medically sensible.) What is the purpose of this request? Is someone telling you you need to provide estimates for instantaneous hazards? -- David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] coxph linear.predictors
coef(fit)*model.matrix(fit)[1,1] age 11.69021 I don't know what that might be and you are not telling us what you think it is. I think this the calculation of the linear predictor, multiplying (the beta*X) I expected that coef(fit)*model.matrix(fit)[1,1]= fit$linear.predictors[1] And since it does not, I am lost. Ultimately I need unconditional death probability from day 1 to day t in the end. Survfit takes too much time as it does a lot of things I do not need. I hope I should be able to get the death prob for every subject from a single call to survfit and then rescale with exp(linear.predictor) for each subject. Just a comment: as there is not deltat argument, I think that coxph assumes deltat=1 so the hazard is the conditional probability over the next time interval ( t(i-1), t(i) ). With dense data sampled daily this is a day interval. Your comments are appreciated. Stephen Bond -Original Message- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Wednesday, October 27, 2010 1:15 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] coxph linear.predictors On Oct 27, 2010, at 12:12 PM, Bond, Stephen wrote: I would like to be able to construct hazard rates (or unconditional death prob) Hazards are not probabilities (since probabilities are constrained to the range [0,1] and hazards are unbounded upward.) for many subjects from a given survfit. This will involve adjusting the ( n.event/n.risk) with (coxph object )$linear.predictors I must be having another silly day as I cannot reproduce the linear predictor: fit - coxph(Surv(futime, fustat) ~ age, data = ovarian) fit$linear.predictors[1] [1] 2.612756 That's the linear predictor (the beta*X) and that particular number only applies to the first case. coef(fit)*model.matrix(fit)[1,1] age 11.69021 I don't know what that might be and you are not telling us what you think it is. The above is based on the help listing for coxph.object coefficients: the coefficients of the linear predictor, which multiply the columns of the model matrix. If the model is over-determined there will be missing values in the vector corresponding to the redundant columns in the model matrix. Also, please comment whether n.event/n.risk The Nelson-Aalen estimator of the cumulative hazard as a function of intervals prior to t is sum( n-event(t)/ n.risk(t)) gives the baseline hazard exp(alpha) ? No. The baseline hazard, as you are calling this, would be an estimate for persons with all covariates = 0, so in this case is for women of age=0. (Not a particularly interpretable result in many situations. The baseline hazard following treated ovarian cancer for neonates is not medically sensible.) What is the purpose of this request? Is someone telling you you need to provide estimates for instantaneous hazards? -- David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] coxph linear.predictors
What you suggest is exactly what I intended. The big problem is that X*beta != linear.predictor # it works for glm That's a big loose end. The survival documentation on p41 implies that the above should be an equality, but it is not. And I could not produce the surv vector using a variety of formulas, so I guess this is the spirit of SAS disguising itself as an R package. Stephen Bond -Original Message- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Wednesday, October 27, 2010 2:17 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] coxph linear.predictors On Oct 27, 2010, at 1:45 PM, Bond, Stephen wrote: coef(fit)*model.matrix(fit)[1,1] age 11.69021 I don't know what that might be and you are not telling us what you think it is. I think this the calculation of the linear predictor, multiplying (the beta*X) I expected that coef(fit)*model.matrix(fit)[1,1]= fit $linear.predictors[1] And since it does not, I am lost. Ultimately I need unconditional death probability from day 1 to day t in the end. Now it sounds like you do want cap-lambda, the cumulative hazard. Let's call it Lambda Since S(t) = exp(-Lambda(t)), you can invert to get Lambda(t) = - log(S(t)) ?survfit Possibly: plot(-log( survfit(fit)$surv )) Survfit takes too much time as it does a lot of things I do not need. I hope I should be able to get the death prob for every subject from a single call to survfit and then rescale with exp(linear.predictor) for each subject. I have run into slowness with that approach myself. The rms package provides survest and Predict methods, but if you wanted to stay in survival, then why not create a survfit()$surv object based on a vector of means and then calculate an exp(newdata) by which you would multiply the S(mean|t)? There is also a predict.coxph function with its type=risk parameter that you might want to check. Just a comment: as there is not deltat argument, I think that coxph assumes deltat=1 so the hazard is the conditional probability over the next time interval ( t(i-1), t(i) ). With dense data sampled daily this is a day interval. That's not the way I understand it. The Lambda() and S() functions only change at the particular event times on a continuous time so they are not assumed to have integer values. Pretty sure that Therneau avoids presenting extraction methods for instantaneous hazards from Cox models. In his book Modeling Survival Data he mostly sticks to discussing cumulative hazard functions rather than estimating an interval or instantaneous hazard. And fortunately he visits rhelp so it's quite possible that my errors in this area will get expert review and correction. Your comments are appreciated. Stephen Bond -Original Message- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Wednesday, October 27, 2010 1:15 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] coxph linear.predictors On Oct 27, 2010, at 12:12 PM, Bond, Stephen wrote: I would like to be able to construct hazard rates (or unconditional death prob) Hazards are not probabilities (since probabilities are constrained to the range [0,1] and hazards are unbounded upward.) for many subjects from a given survfit. This will involve adjusting the ( n.event/n.risk) with (coxph object )$linear.predictors I must be having another silly day as I cannot reproduce the linear predictor: fit - coxph(Surv(futime, fustat) ~ age, data = ovarian) fit$linear.predictors[1] [1] 2.612756 That's the linear predictor (the beta*X) and that particular number only applies to the first case. coef(fit)*model.matrix(fit)[1,1] age 11.69021 I don't know what that might be and you are not telling us what you think it is. The above is based on the help listing for coxph.object coefficients: the coefficients of the linear predictor, which multiply the columns of the model matrix. If the model is over-determined there will be missing values in the vector corresponding to the redundant columns in the model matrix. Also, please comment whether n.event/n.risk The Nelson-Aalen estimator of the cumulative hazard as a function of intervals prior to t is sum( n-event(t)/ n.risk(t)) gives the baseline hazard exp(alpha) ? No. The baseline hazard, as you are calling this, would be an estimate for persons with all covariates = 0, so in this case is for women of age=0. (Not a particularly interpretable result in many situations. The baseline hazard following treated ovarian cancer for neonates is not medically sensible.) What is the purpose of this request? Is someone telling you you need to provide estimates for instantaneous hazards? -- David Winsemius, MD West Hartford, CT David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https
Re: [R] coxph linear.predictors
I think I found the first piece of it. Typical SAS: fit$linear.predictors[1]-coef(fit)*(ovarian$age[1]-mean(ovarian$age)) age -4.4408920985e-16 Well misrepresented in several places... Stephen Bond -Original Message- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Wednesday, October 27, 2010 1:15 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] coxph linear.predictors On Oct 27, 2010, at 12:12 PM, Bond, Stephen wrote: I would like to be able to construct hazard rates (or unconditional death prob) Hazards are not probabilities (since probabilities are constrained to the range [0,1] and hazards are unbounded upward.) for many subjects from a given survfit. This will involve adjusting the ( n.event/n.risk) with (coxph object )$linear.predictors I must be having another silly day as I cannot reproduce the linear predictor: fit - coxph(Surv(futime, fustat) ~ age, data = ovarian) fit$linear.predictors[1] [1] 2.612756 That's the linear predictor (the beta*X) and that particular number only applies to the first case. coef(fit)*model.matrix(fit)[1,1] age 11.69021 I don't know what that might be and you are not telling us what you think it is. The above is based on the help listing for coxph.object coefficients: the coefficients of the linear predictor, which multiply the columns of the model matrix. If the model is over-determined there will be missing values in the vector corresponding to the redundant columns in the model matrix. Also, please comment whether n.event/n.risk The Nelson-Aalen estimator of the cumulative hazard as a function of intervals prior to t is sum( n-event(t)/ n.risk(t)) gives the baseline hazard exp(alpha) ? No. The baseline hazard, as you are calling this, would be an estimate for persons with all covariates = 0, so in this case is for women of age=0. (Not a particularly interpretable result in many situations. The baseline hazard following treated ovarian cancer for neonates is not medically sensible.) What is the purpose of this request? Is someone telling you you need to provide estimates for instantaneous hazards? -- David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] specify strata in survfit
Is it possible to get survfit to produce the survival line for a single strata like preddow - survfit(modall,newdata=newdat,se.fit=F,strata=2) # the strata argument is being ignored in the call above Or even get a more economical/faster calculation of the hazard directly from the coxph object?? Stephen [[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] see code of plot.survfit
How can I expose the code behind plot.survfit?? Thanks a lot. Stephen B [[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] aggregate with cumsum
Gabor,
You are suggesting some very advanced usage that I do not understand, but it
seems this is not what I meant when I said loop.
I have a df with 47k rows and each of these is fed to a 'predict' which will
output about 62 rows, so the number of groups is very large and I implied that
I would go through the 47k x 62 rows with
For (jj in (set of 47k values)) # tmp.df=big.df[big.df$group==jj,] to subset
# and then sum
Which is very slow. I discovered that even creating the dataset is super slow
as I use write.table
The clogging comes from
write.table(tmp,predcom.csv,row.names=FALSE,col.names=FALSE,append=TRUE,sep=',')
Can anybody suggest a faster way of appending to a text file??
All comments are appreciated.
Stephen B
-Original Message-
From: Gabor Grothendieck [mailto:ggrothendi...@gmail.com]
Sent: Tuesday, October 12, 2010 4:16 PM
To: Bond, Stephen
Cc: r-help@r-project.org
Subject: Re: [R] aggregate with cumsum
On Tue, Oct 12, 2010 at 1:40 PM, Bond, Stephen stephen.b...@cibc.com wrote:
Hello everybody,
Data is
myd - data.frame(id1=rep(c(a,b,c),each=3),id2=rep(1:3,3),val=rnorm(9))
I want to get a cumulative sum over each of id1. trying aggregate does not
work
myd$pcum - aggregate(myd[,c(val)],list(orig=myd$id1),cumsum)
Please suggest a solution. In real the dataframe is huge so looping with for
and subsetting is not a great idea (still doable, though).
Looping can be slow but its not necessarily so. Here are three
approaches to using ave with cumsum to solve this problem. The
benchmark shows that the loop is actually the fastest:
N - 1e4
k - 10
myd - data.frame(id1=rep(letters[1:k],each=N),id2=rep(1:k,N),val=rnorm(k*N))
library(rbenchmark)
benchmark(order = relative, replications = 100,
loop = { loop - myd
for(i in 2:3) loop[, i] - ave(myd[, i], myd[, 1], FUN = cumsum)
},
nonloop1 = { nonloop1 - transform(myd,
id2 = ave(id2, id1, FUN = cumsum),
val = ave(val, id1, FUN = cumsum)
)},
nonloop2 = {
f - function(i) ave(myd[, i], myd[, 1], FUN = cumsum)
nonloop2 - replace(myd, 2:3, lapply(2:3, f))
}
)
identical(loop, nonloop1)
identical(loop, nonloop2)
The output on my laptop is:
test replications elapsed relative user.self sys.self user.child sys.child
1 loop 1008.52 1.00 8.07 0.10 NANA
3 nonloop2 1008.94 1.049296 8.29 0.17 NANA
2 nonloop1 100 11.65 1.367371 10.71 0.22 NANA
--
Statistics Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] aggregate with cumsum
Hello everybody, Data is myd - data.frame(id1=rep(c(a,b,c),each=3),id2=rep(1:3,3),val=rnorm(9)) I want to get a cumulative sum over each of id1. trying aggregate does not work myd$pcum - aggregate(myd[,c(val)],list(orig=myd$id1),cumsum) Please suggest a solution. In real the dataframe is huge so looping with for and subsetting is not a great idea (still doable, though). Thank you Stephen B [[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] Comparing COXPH models, one with age as a continuous variable, one with age as a three-level factor
Totally agreed. I made a mistake in calling the categorization a GAM. If we apply a step function to the continuous age we get a limited range ordinal variable. Categorizing is creating several binary variables from the continuous (with treatment contrasts). Stephen B -Original Message- From: Frank Harrell [mailto:harre...@gmail.com] On Behalf Of Frank Harrell Sent: Thursday, September 02, 2010 4:23 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] Comparing COXPH models, one with age as a continuous variable, one with age as a three-level factor On Thu, 2 Sep 2010, stephenb wrote: sorry to bump in late, but I am doing similar things now and was browsing. IMHO anova is not appropriate here. it applies when the richer model has p more variables than the simpler model. this is not the case here. the competing models use different variables. A simple approach is to have the factor variable in the model and to formally test for added information given by the continuous variable (linear, quadratic, spline, etc). AIC could also be used. you are left with IC. by transforming a continuous variable into categorical you are smoothing, which is the idea of GAM. if you look at what is offered in GAMs you may find better approximations f(age) as well as tools for testing among different f(age) transformations. I don't follow that comment. Smoothing uses the full continuous variable to begin with. A restricted cubic spline function in age is a simple approach. E.g.: require(rms) dd - datadist(mydata); options(datadist='dd') f - cph(Surv(dtime,death) ~ rcs(age,4) + sex, data=mydata) plot(Predict(f, age)) Note that you can always expect the categorized version of age not to fit the data except sometimes when behavior is dictated by law (driving, drinking, military service, medicare). Frank regards. S. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] sqldf syntax
I had checked those references before posting, actually. SQLite has a very limited implementation of the standard. To do a single table update I would not go to sql. It's easy enough to do in R. The problem is when I need to do an update from a left outer join, which I had to do with sqlSave (to a mySQL table), then sqlQuery, then sqlFetch. sqlSave is amazingly slow, takes half an hour. (Would never do that at home :-) just too lazy to write a formal table def and use load data infile from a csv dump. Also not happy with Dates becoming years in the transition. Will check the other suggestion about data.table and report. Cheers everybody. Stephen B -Original Message- From: Gabor Grothendieck [mailto:ggrothendi...@gmail.com] Sent: Thursday, August 26, 2010 4:26 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] sqldf syntax On Thu, Aug 26, 2010 at 2:31 PM, Bond, Stephen stephen.b...@cibc.com wrote: Please correct the following sqldf(update esc left join forwagg on esc.ym=forwagg.Date set esc.ri2=forwagg.N1 where esc.age=12,select * from main.esc) Error in sqliteExecStatement(con, statement, bind.data) : RS-DBI driver: (error in statement: near left: syntax error) 1. sqldf takes one sql argument whereas the above has two sql arguments; however, the one argument may be a vector of sql commands. See ?sqldf and the examples on the sqldf home page http://sqldf.googlecode.com 2. there is an error in the syntax of your update statement. For correct syntax see the sqlite site: http://sqlite.org/lang_update.html -- Statistics Software Consulting GKX Group, GKX Associates Inc. tel: 1-877-GKX-GROUP email: ggrothendieck at gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] sqldf syntax
Please correct the following sqldf(update esc left join forwagg on esc.ym=forwagg.Date set esc.ri2=forwagg.N1 where esc.age=12,select * from main.esc) Error in sqliteExecStatement(con, statement, bind.data) : RS-DBI driver: (error in statement: near left: syntax error) Thanks. Stephen [[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] transfer glm model from SAS
Please, tell what is the best way to create an R glm object with parameters etc estimated in SAS? I have a large dataset and bigglm fails to converge, so estimation is done is SAS. However, there are a lot of predictions that are much more easily done in R and I would like to use glm.predict and the flexibility of R. bigglm ends with an error, so no object is created. Thanks everybody. Stephen [[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] forcing a zero level in contr.sum
Peter, I am asked to show seasonal effects vs a general averaged level without any seasonal effects. That's why contr.sum and not treatment. The estimates show that some coefs are not stat different from 0, which is also the enforced total average. So, my understanding of model selection leads me to force certain seasons to be zero. For example: Y= a0 + b1 + Xb | factor=1 Y= a0 +Xb| factor=2 Y= a0 + b2 + Xb | factor=3 b1+b2=0 | which in this case means b1=-b2 but with more levels not necessarily Stephen Bond | -Original Message- From: Peter Dalgaard [mailto:pda...@gmail.com] Sent: Thursday, July 08, 2010 4:09 AM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] forcing a zero level in contr.sum Bond, Stephen wrote: I need to use contr.sum and observe that some levels are not statistically different from the overall mean of zero. What is the proper way of forcing the zero estimate? It seems the column corresponding to that level should become a column of zeros. Is there a way to achieve that without me constructing the design matrix? As Berwin (indirectly) points out, you probably overlooked the how.many argument to C(). However, are you _sure_ that you want to do this? This is like testing that one treatment is exactly equal to the average of all other treatments, which is a rather strange hypothesis. -- Peter Dalgaard Center for Statistics, Copenhagen Business School 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] forcing a zero level in contr.sum
I need to use contr.sum and observe that some levels are not statistically different from the overall mean of zero. What is the proper way of forcing the zero estimate? It seems the column corresponding to that level should become a column of zeros. Is there a way to achieve that without me constructing the design matrix? Thank you. Stephen Bond [[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] forcing a zero level in contr.sum
Clarifying my question: options(contrasts = c(contr.sum, contr.poly)) contrasts(fixw$snconv) [,1] [,2] [,3] [,4] [,5] 0110000 0301000 0500100 0600010 0700001 09 -1 -1 -1 -1 -1 I need to force the coefficient on level 03 to be zero. Thank you. Stephen Bond -Original Message- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Wednesday, July 07, 2010 3:44 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] forcing a zero level in contr.sum On Jul 7, 2010, at 3:13 PM, Bond, Stephen wrote: I need to use contr.sum and observe that some levels are not statistically different from the overall mean of zero. What is the proper way of forcing the zero estimate? It seems the column corresponding to that level should become a column of zeros. Is there a way to achieve that without me constructing the design matrix? Thank you. lm( formula = z ~ x + y + 0, ...) _might_ do something close to what you want. David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] forcing a zero level in contr.sum
Please, do not post if you do not know the answer. People will see this has answers and skip. I tried with mat1=contrasts(fixw$snconv) mat1=mat1[,-2] summary(frm2sum - glm(resp.frm ~ C(snconv,contr=mat1)+mprime+mshape,data=fixw,family=quasibinomial)) the unwanted level is still there. Unbelievable. Stephen Bond -Original Message- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Wednesday, July 07, 2010 4:15 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] forcing a zero level in contr.sum On Jul 7, 2010, at 4:04 PM, Bond, Stephen wrote: Clarifying my question: options(contrasts = c(contr.sum, contr.poly)) contrasts() [,1] [,2] [,3] [,4] [,5] 0110000 0301000 0500100 0600010 0700001 09 -1 -1 -1 -1 -1 I need to force the coefficient on level 03 to be zero. ?factor ?relevel Perhaps worth a try: fixw$snconv - relevel(fixw$snconv, ref=03) (But I wonder if using contr.sum will ever generally satisfy that goal, since contr.sum calculates the difference from a grand mean and this will only work if a) the GM=0 and b) there is only one term on the RHS of the model, and c) probably a bunch of other restrictions.) -- David. Thank you. Stephen Bond -Original Message- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Wednesday, July 07, 2010 3:44 PM To: Bond, Stephen Cc: r-help@r-project.org Subject: Re: [R] forcing a zero level in contr.sum On Jul 7, 2010, at 3:13 PM, Bond, Stephen wrote: I need to use contr.sum and observe that some levels are not statistically different from the overall mean of zero. What is the proper way of forcing the zero estimate? It seems the column corresponding to that level should become a column of zeros. Is there a way to achieve that without me constructing the design matrix? Thank you. lm( formula = z ~ x + y + 0, ...) _might_ do something close to what you want. David Winsemius, MD West Hartford, CT David Winsemius, MD West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] correlation in glm residuals
Is there a library dealing with correlation in the residuals of a glm? I have bin3alt -glm(respalt~ t+sn+c5.vrm,data=dfalt,family=quasibinomial) bin3alt Call: glm(formula = respalt ~ t + sn + c5.vrm, family = quasibinomial, data = dfalt) Coefficients: (Intercept) t2 t3 t4 t5 t6 sn2 sn3 c5.vrm -3.35957 1.81455 0.96161 -0.37701 -2.32657 -3.75074 0.24266 0.39056 0.06673 Degrees of Freedom: 230 Total (i.e. Null); 222 Residual Null Deviance: 107000 Residual Deviance: 2290AIC: NA dfalt$pears - residuals(bin3alt,type=pearson) arima(dfalt$pears[dfalt$t==4],order=c(1,0,0))) Call: arima(x = dfalt$pears[dfalt$t == 4], order = c(1, 0, 0)) Coefficients: ar1 intercept 0.6333-0.5091 s.e. 0.1257 1.0490 Not all levels of the t factor, show correlation, but some do. The factor is not a random effect it is month of ageing. Also, if I use the Cochrane Orcutt manually, should I use response or pearson residuals? I know of lme, but think it requires a random effect. Stephen Bond [[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] how to remove interactions of factor with continuous var
I need to remove certain interactions and keep only the one between the second level of the factor and the continuous var t2 bin4 - glm(resp2~ t*t2+c5.vrm,data=dfa,family=quasibinomial) summary(bin4) Call: glm(formula = resp2 ~ t * t2 + c5.vrm, family = quasibinomial, data = dfa) Deviance Residuals: Min 1Q Median 3Q Max -6.5464 -3.0720 -1.8135 0.4896 28.9207 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) -1.492379 0.147779 -10.099 2e-16 *** t2 1.381409 0.193110 7.153 1.19e-11 *** t3 0.223979 0.275389 0.813 0.416898 t4 -1.751518 0.495391 -3.536 0.000494 *** t5 -2.455172 0.913674 -2.687 0.007747 ** t6 -2.347207 1.131134 -2.075 0.039120 * t2 -0.030819 0.013144 -2.345 0.019914 * c5.vrm 0.427088 0.140085 3.049 0.002574 ** t2:t2 -0.015909 0.007292 -2.182 0.030174 * t3:t2 -0.016170 0.010479 -1.543 0.124218 t4:t20.024352 0.016900 1.441 0.151009 t5:t2 -0.012071 0.034230 -0.353 0.724680 t6:t2 -0.041797 0.047106 -0.887 0.375880 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for quasibinomial family taken to be 41.75819) Trying with - t3:t2 etc did not work neither did putting quotes around the terms. I need to have the continuous var c5.vrm same for all levels of the factor so estimating separate equations for each level will not work. I could use offset, but seems quite an ugly solution, there should be a way to specify exactly which interactions are desired. Thank you all. Stephen Bond [[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] change axis labels in persp
I need to change one of the axis produced by persp by giving it labels for the tickmarks. The dimension has months (a factor) so the default decimals don't look good. The graph will finally become a pdf or emf file to be embedded in a document, so persp3d will not work, I think. Thank you all. Stephen [[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] how to specify between group variance in lme
Hello r-help, I am using lme with two specs for the variance func varComb(varFixed(~1/n)),varPower(~Age)) this produces worse forecasts than the lm model with simple weights=n I think due to the fact that the lme spec works on variance inside the group. I need to show it that 1/n scales the variance between groups. Is that possible? I cannot disclose my dataset, but could post plots if that is possible somehow, let me know. Anova shows that everything in random=~poly(age,2) is significant. Thank you all very much. Stephen [[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] how to specify between group variance in lme
Clarification: Lm is much better than the base forecast from lme level=0, Level=1 produces a much tighter fit than lm. I was expecting that level=0 would produce something very close to lm, but it does not. [[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] definition of AIC and BIC in gls
Hello everybody, Please help with connecting the AIC and BIC numbers printed by summary.gls to the logLik number. 1. is the logLik number the true ML or density scaling constants have been omitted? 2. what is the formula for calculating the AIC and BIC from logLik (and how can I see it)? I tried printing summary.gls but it says object not found. Thank you very much. Stephen [[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] should months be an ordered factor?
This is not a purely stats question as it depends on the implementation of mixed models in R. The help says that ordered factors are treated differently than unordered, which is not very informative. More explicitly: does ordering simply imply that February follows January or that February has got more than January of sth that they both have? Or pls recommend another way of treating monthly effects. Thank you [[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] how to specify two variance effects in gls
Hello everybody, I have a dataset where each row has number of subjects and that gives me natural weights for the variance function. Additionally I see that variance increases with Age, which is a regressor. So using gls I have weights=varFixed(~1/n) but don't know how to include the extra effect of the regressor. Fitted values show a quadratic curve vs. age, not sure if that helps. Thanks everybody. Stephen [[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.
