[Rd] Bug in varIdent (nlme package) (PR#9765)

From: <agalecki_at_umich.edu>
Date: Fri, 29 Jun 2007 01:18:25 +0200 (CEST)


This is a multi-part message in MIME format.

--------------000706080709050507050106
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit

Hello,

  1. The following code (below) and output <varIdentOutput.txt> illustrate that argument fixed= in varIdent function is ignored during the execution.
  2. File named <RevisedvarIdent.txt> shows how to fix the problem. The changed line in this file is labeled using comment " # !!! This line changed".

Thank you

Andrzej Galecki

### Code used for illustration starts here #### library(nlme)
sessionInfo()

Orth <- subset(Orthodont, Subject %in% c("M01","F01")) # Argument fixed in varIdent is ignored
vf1.incorrect <- varIdent(form=~1|Sex,fixed=c(Female=0.5)) vf1I <- Initialize(vf1.incorrect, data=Orth) varWeights(vf1I)

# Tentative varIdent function (revised)
source("C:/temp/RevisedvarIdent.txt")

vf2.correct <- varIdent(form=~1|Sex,fixed=c(Female=0.5)) vf2I <- Initialize(vf2.correct, data=Orth) coef(vf2I)
varWeights(vf2I)
# rm(varIdent) # Remove tentative varIdent function

### Code ends here ####

--------------000706080709050507050106
Content-Type: text/plain;
 name="RevisedvarIdent.txt"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="RevisedvarIdent.txt"

varIdent <-
function (value = numeric(0), form = ~1, fixed = NULL) {

    if (is.null(getGroupsFormula(form))) {

        value <- numeric(0)
        attr(value, "fixed") <- NULL

    }
    else {
        if ((lv <- length(value)) > 0) {
            if (is.null(grpNames <- names(value)) && (lv > 1)) {
                stop("Initial values must have group names in varIdent")

}
value <- unlist(value) if (any(value <= 0)) { stop("Initial values for \"varIdent\" must be > 0.")
}
value <- log(value) } else grpNames <- NULL attr(value, "groupNames") <- grpNames if (!is.null(fixed)) { # !!! This line changed fix <- attr(value, "fixed") <- log(unlist(fixed)) if (is.null(fixNames <- names(fix))) { stop("Fixed parameters must have names in varIdent")
}
if (!is.null(attr(value, "groupNames"))) { attr(value, "groupNames") <- c(attr(value, "groupNames"), fixNames)
}
}

    }
    attr(value, "formula") <- asOneSidedFormula(form)     class(value) <- c("varIdent", "varFunc")     value
}

--------------000706080709050507050106
Content-Type: text/plain;
 name="varIdentOutput.txt"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="varIdentOutput.txt"

> library(nlme)
> sessionInfo()
R version 2.5.0 (2007-04-23)
i386-pc-mingw32

locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

attached base packages:

[1] "stats"     "graphics"  "grDevices" "utils"     "datasets"  "methods"  
[7] "base"     

other attached packages:

    nlme
"3.1-80"
>

> # Argument fixed in varIdent is ignored
> Orth <- subset(Orthodont, Subject %in% c("M01","F01"))
> # Argument fixed in varIdent is ignored
> vf1.incorrect  <- varIdent(form=~1|Sex,fixed=c(Female=0.5))
> vf1I <- Initialize(vf1.incorrect, data=Orth)
> varWeights(vf1I)
  Male   Male   Male   Male Female Female Female Female 
     1      1      1      1      1      1      1      1 
>
> # Tentative varIdent function (revised)
> source("C:/temp/RevisedvarIdent.txt")
> 
> vf2.correct  <- varIdent(form=~1|Sex,fixed=c(Female=0.5))
> vf2I <- Initialize(vf2.correct, data=Orth)
> varWeights(vf2I)
  Male   Male   Male   Male Female Female Female Female 
     1      1      1      1      2      2      2      2 
> rm(varIdent) # Remove tentative varIdent function

--------------000706080709050507050106--



R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Fri 29 Jun 2007 - 15:35:24 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sat 30 Jun 2007 - 00:35:37 GMT.

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