R devel - Mar 2004 - termplot has problems with a single term (fix included) (PR#6679)

On Thu, 18 Mar 2004 k.hansen@biostat.ku.dk wrote:
> The bug exists on R-1.9.0-alpha compiled the 10/3.
>
> Termplot has a problem if either the model only contains a single term
> or if asked to plot a single term. In addition there are problems with
> the option se = TRUE.
I can't reproduce this in either R-devel or 1.8.1, and termplot hasn't
changed since January.

I do
example(termplot)
termplot(model, terms="z")
termplot(model, terms="z", se=TRUE)
termplot(model, terms=1)
termplot(model, terms="ns(x, 6)")
termplot(model, terms="ns(x, 6)", se=TRUE)
termplot(model, terms=2)

mm<-glm(y~x)
termplot(mm)
termplot(mm,se=TRUE)
mm<-glm(y~x-1)
termplot(mm)
termplot(mm,se=TRUE)

without any error messages.  I also note that termplot() doesn't
(anymore) start with the line you say it starts with, though it did until
January.

Could you provide an example?

	-thomas
> Analysis: termplot starts with
>
>   terms <- if (is.null(terms))
>     predict(model, type = "terms", se = se)
>   else predict(model, type = "terms", se = se, terms = terms)
>   n.tms <- ncol(tms <- as.matrix(if (se)
>                                  terms$fit
>   else terms))
>
> In this case terms is now a matrix of fitted values, with a single
> column for each term in the model or in the terms argument. Because
> the predict function returns a vector if called with a single term,
> the matrix tms has no column names. This maskes the lines
>   nmt <- colnames(tms)
>   cn <- parse(text = nmt)
> fail.
>
> Solution: manually add the column names in case the number of columns
> of tms == 1. Beware that the terms argument is overwritten in the
> first line of the function (Is that a good style? Perhaps the body of
> the function needs another name for the terms object defined in the
> very first line). This means we manually need to save the terms name
> first (code below)
>
> Fixing this bug uncovers another one: in case of a single term in the
> model formula, and with the option se = TRUE, the terms objects
> consists of a list with elements fit and se.fit. In the code above
> explicit care is taken to ensure that the fit element is a matrix by
> using
>   tms <- as.matrix(if (se) terms$fit
>                    else terms)
> The element se.fit needs to be a matrix as well. I suggest a final fix
> to be
>
> {
> *  terms.names <- terms
>    terms <- if (is.null(terms))
>      predict(model, type = "terms", se = se)
>    else predict(model, type = "terms", se = se, terms = terms)
>    n.tms <- ncol(tms <- as.matrix(if (se)
>                                   terms$fit
>    else terms))
> *  if(n.tms==1)
> *    colnames(tms) <- if(is.null(terms)) attr(model$terms,
"term.labels")
> *                        else terms.names
> *  if(se)
> *    terms$se.fit <- as.matrix(terms$se.fit)
>
> * = added line.
>
>
> The fix above is somewhat ugly. Perhaps (as mentioned above) the
> object terms needs renaming (terms.something?), in which case the
> first line is redundant.
> --
> Kasper Daniel Hansen, Research Assistent
> Department of Biostatistics, University of Copenhagen
>
> ______________________________________________
> R-devel@stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-devel
>
Thomas Lumley			Assoc. Professor, Biostatistics
tlumley@u.washington.edu	University of Washington, Seattle

tlumley@u.washington.edu writes:
> On Thu, 18 Mar 2004 k.hansen@biostat.ku.dk wrote:
> 
> > The bug exists on R-1.9.0-alpha compiled the 10/3.
> >
> > Termplot has a problem if either the model only contains a single term
> > or if asked to plot a single term. In addition there are problems with
> > the option se = TRUE.
> 
> I can't reproduce this in either R-devel or 1.8.1, and termplot
hasn't
> changed since January.
> 
> I do
> example(termplot)
...
> without any error messages.  I also note that termplot() doesn't
> (anymore) start with the line you say it starts with, though it did until
> January.
> 
> Could you provide an example?
The original story was this (which I saw but lacked the time to analyze):

library(survival)
library(MASS)
library(splines)
data(Melanoma)
men<-Melanoma[Melanoma$status<3,]
y<-coxph(Surv(time,status) ~ sex + ns(thickness,df=4),data=men)
z<-coxph(Surv(time,status) ~ ns(thickness,df=4),data=men)

termplot(y)         #OK
termplot(z)         # Prints a ? prompt (from parse(text=nmt) w/nmt==NULL)
termplot(y,terms=2) # Same

This can still be seen with the current devel version.

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk)             FAX: (+45) 35327907