[R] my error with augPred

Petr Pikal petr.pikal at precheza.cz
Fri Sep 8 08:52:36 CEST 2006 

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 at pdf.com>
To:             	Petr Pikal <petr.pikal at precheza.cz>
Copies to:      	r-help at stat.math.ethz.ch, Douglas Bates <bates at 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 at pdf.com>
> > To:             	Petr Pikal <petr.pikal at precheza.cz>
> > Copies to:      	r-help at 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 at precheza.cz
> >>>
> >>> ______________________________________________
> >>> R-help at 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 at 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 at precheza.cz
> >
> >
> 
> ______________________________________________
> R-help at 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 at precheza.cz

More information about the R-help mailing list