Re: [R] my error with augPred

From: Petr Pikal <petr.pikal_at_precheza.cz>
Date: Fri 08 Sep 2006 - 06:52:36 GMT

Hi Spencer

Thank you for your reply. I tried as you shad suggested and it seems to me that problem comes from this piece of code

contr <- object$contrasts
Browse[1]>
debug: for (i in names(dataMix)) {

    if (inherits(dataMix[, i], "factor") && !is.null(contr[[i]])) {

        levs <- levels(dataMix[, i])
        levsC <- dimnames(contr[[i]])[[1]]
        if (any(wch <- is.na(match(levs, levsC)))) {
            stop(paste("Levels", paste(levs[wch], collapse = ","), 
                "not allowed for", i))
        }
        attr(dataMix[, i], "contrasts") <- contr[[i]][levs, , 
            drop = FALSE]

    }
}

especially from levs and levsC comparison.

levs are letters[1:3] and levsC is NULL because object$contrasts is

Browse[1]> object$contrasts
$x1
[1] "contr.treatment"

As a statistian amateur I can not say if it si a bug or not, but it seems to me that if in case of factors object$contrasts is always "contr.treatment" there is no way how to match it with actual levels of contrasts.

Someone more experienced has to decide about nature of this feature.

Best regards
Petr

On 6 Sep 2006 at 23:10, Spencer Graves wrote:

Date sent:      	Wed, 06 Sep 2006 23:10:17 -0700
From:           	Spencer Graves <spencer.graves@pdf.com>
To:             	Petr Pikal <petr.pikal@precheza.cz>
Copies to:      	r-help@stat.math.ethz.ch, Douglas Bates <bates@stat.wisc.edu>
Subject:        	Re: [R] my error with augPred


> Thank you for providing such a complete, self contained example.
>
> I found that 'predict.nlme' does not like a factor in the 'fixed' > argument as you used it, "fixed=list(Asym~x1, R0+lrc~1)". To see > this, I added 'x1.' as a numeric version of the factor 'x1' and reran > it successfully: > > fm2.<-update(fm1, fixed=list(Asym~x1., R0+lrc~1), > start=c(103,0,-8.5,-3)) aP2. <- augPred(fm2.) plot(aP2.) >
> Unfortunately, it looks like this work-around won't help you
> with
> your original problem, because there, the counterpart to 'x1' is an > ordered factor with more than 2 levels. >
> The error message refers to 'predict.nlme'. I know no reason
> why
> 'predict.nlme' shouldn't work with a factor with more than 2 levels in > this context. If it were my problem and it was sufficiently > important, I would make a local copy of 'predict.nlme' as follows: >
> predict.nlme <- getAnywhere("predict.nlme")
>
> Then I'd use 'debug(nlme:::predict.nlme)' to walk through the
> problem example line by line until I figured out what I had to change > to make this work. >
> I hesitate to use the "B" word, but I think it might be
> appropriate to file a bug report on this; perhaps someone else will > do that. >
> I'm sorry I couldn't solve your original problem. With luck,
> someone else will convert this example into a fix to the code.
> Spencer Graves
> > Petr Pikal wrote: > > Hallo > > > > thank you for your response. I am not sure but maybe fixed effects > > cannot be set to be influenced by a factor to be able to use > > augPred. > > > > lob<-Loblolly[Loblolly$Seed!=321,] > > set.seed(1) > > lob<-data.frame(lob, x1=sample(letters[1:3], replace=T)) # add a > > #factor > > lob<-groupedData(height~age|Seed, data=lob) > > fm1 <- nlme(height ~ SSasymp(age, Asym, R0, lrc), > > data = lob, > > fixed = Asym + R0 + lrc ~ 1, > > random = Asym ~ 1, > > start = c(Asym = 103, R0 = -8.5, lrc = -3.3)) > > > > fm2<-update(fm1, fixed=list(Asym~x1, R0+lrc~1), > > start=c(103,0,-8.5,-3)) > > ^^^^^^^ > > and > > > > plot(augPred(fm2)) > > > > Throws an error. > > So it is not possible to use augPred with such constructions. > > > > Best regards. > > Petr Pikal > > > > On 2 Sep 2006 at 17:58, Spencer Graves wrote: > > > > Date sent: Sat, 02 Sep 2006 17:58:05 -0700 > > From: Spencer Graves <spencer.graves@pdf.com> > > To: Petr Pikal <petr.pikal@precheza.cz> > > Copies to: r-help@stat.math.ethz.ch > > Subject: Re: [R] my error with augPred > > > > > >> <comments in line> > >> > >> Petr Pikal wrote: > >> > >>> Dear all > >>> > >>> I try to refine my nlme models and with partial success. The model > >>> is refined and fitted (using Pinheiro/Bates book as a tutorial) > >>> but when I try to plot > >>> > >>> plot(augPred(fit4)) > >>> > >>> I obtain > >>> Error in predict.nlme(object, value[1:(nrow(value)/nL), , drop = > >>> FALSE], : > >>> Levels (0,3.5],(3.5,5],(5,7],(7,Inf] not allowed for > >>> vykon.fac > >>> > >>> > >>> Is it due to the fact that I have unbalanced design with not all > >>> levels of vykon.fac present in all levels of other explanatory > >>> factor variable? > >>> > >>> > >> I don't know, but I'm skeptical. > >> > >>> I try to repeat 8.19 fig which is OK until I try: > >>> > >>> fit4 <- update(fit2, fixed = list(A+B~1,xmid~vykon.fac, scal~1), > >>> start = c(57, 100, 700, rep(0,3), 13)) > >>> > >>> I know I should provide an example but maybe somebody will be > >>> clever enough to point me to an explanation without it. > >>> > >>> > >> I'm not. > >> > >> To answer these questions without an example from you, I'd have to > >> make up my own example and try to see if I could replicate the > >> error messages you report, and I'm not sufficiently concerned about > >> this right now to do that. > >> > >> Have you tried taking an example from the book and deleting certain > >> rows from the data to see if you can force it to reproduce your > >> error? > >> > >> > >> Alternatively, have you tried using 'debug' to trace through the > >> code line by line until you learn enough of what it's doing to > >> answer your question? > >> > >> Spencer Graves > >> > >>> nlme version 3.1-75 > >>> SSfpl model > >>> R 2.4.0dev (but is the same in 2.3.1), W2000. > >>> > >>> Thank you > >>> Best regards. > >>> > >>> Petr PikalPetr Pikal > >>> petr.pikal@precheza.cz > >>> > >>> ______________________________________________ > >>> R-help@stat.math.ethz.ch mailing list > >>> https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list > >> https://stat.ethz.ch/mailman/listinfo/r-help > >> PLEASE do read the posting guide > >> http://www.R-project.org/posting-guide.html and provide commented, > >> minimal, self-contained, reproducible code. > >> > > > > Petr Pikal > > petr.pikal@precheza.cz > > > > > > ______________________________________________ > R-help@stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide

> http://www.R-project.org/posting-guide.html and provide commented, > minimal, self-contained, reproducible code.

Petr Pikal
petr.pikal@precheza.cz



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Fri Sep 08 16:58:24 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Fri 08 Sep 2006 - 07:30:06 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.