Bug in coef<-.varIdent method (nlme package) (PR#9831)

Tue, Aug 7, 2007 12:44 PM

This is a multi-part message in MIME format.
--------------040502080208060001050807
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit

Hello,

1. It appears that "coef<-.varIdent" method does not work properly in 
some instances.
Execution error:

Error in `coef<-.varIdent`(`*tmp*`, value = c(11, 12)) :
        Cannot change the length of the varIdent parameter after 
initialization


occurs when  "coef<-.varIdent" is applied to an initialized object of 
class varIdent  with some of the coefficients  being _fixed_.
Attached files: 'varIdentOrthoEmail.txt'  and 'varIdentOrthoEmail.Rout' 
illustrate the problem.

2. The code for  "coef<-.varIdentX" method in  
'varIdentmethodsRevised.txt' file  illustrates how to fix this problem.

3. Specifically, to fix the problem the line

       if (length(value) != nGroups  - 1)  

   in  the "coef<-.varIdent" method should be replaced
   with the following two lines :

        nFixed  <- sum(as.numeric(attr(object,"whichFix")))  # inserted 
new line
        if (length(value) != nGroups - nFixed - 1) {               # 
modified original line

4. Note: Although I am using lme "3.1-80", the  related problem PR#9765 
is fixed manually by over-writing varIdent function.

Thank you

Andrzej Galecki



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



ls() 
require(nlme)
sessionInfo()


## New class varIdentX and methods defined.
### coef<-.varIdentX illustrates how to fix a bug in coef<-.varIdent
# Note: PR#9765 fixed in varIdent() and varIdentX() 
source("C:/temp/varIdentmethodsRevised.txt")
ls()



#### Chunk 1: Everything is OK here
# value= and fixed= arguments
# variance component at age=12 is fixed

val <- c("10"=1.10,"14"=1.14)
vf <- varIdent(value=val, form=~1|age, fixed=c("12"=1.12))
vfi <- Initialize(vf,Orthodont)
str(vfi)
coef(vfi)
coef(vfi, unconstrained = FALSE, allCoef = TRUE)  
vfiCopy <- vfi        # copy of an initialized object

#### Chunk 2: Bug in coef<-.varIdent illustrated

length(vfiCopy)             # length is 2
coef(vfiCopy) <- c(11,12)   # Execution error
# Error in `coef<-.varIdent`(`*tmp*`, value = c(11, 12)) : 
#        Cannot change the length of the varIdent parameter after initialization

#### Chunk 3: Consequently the same execution error in gls

gls.error  <- gls(distance ~ age, 
                weights  = vfi,
                data=Orthodont)

#Error in `coef<-.varIdent`(`*tmp*`, value = c(0.095310179804325, 0.131028262406404 : 
#         Cannot change the length of the varIdent parameter after initialization


### Chunk 1A: 
  
val <- c("10"=1.10,"14"=1.14)
vf <- varIdentX(value=val, form=~1|age, fixed=c("12"=1.12))
vfi <- Initialize(vf,Orthodont)
str(vfi)              # Note: class varIdentX
coef(vfi)
coef(vfi, unconstrained = FALSE, allCoef = TRUE)  
vfiCopy <- vfi        # copy of an initialized object

#### Chunk 2A: Bug in coef<-.varIdent  *** corrected ***

length(vfiCopy)             # length is 2
coef(vfiCopy) <- c(11,12)   # NO Execution error

#### Chunk 3A: No execution error in gls

gls.noerror  <- gls(distance ~ age, 
                weights  = vfi,
                data=Orthodont)



--------------040502080208060001050807
Content-Type: text/plain;
 name="varIdentXOrthoEmail.Rout"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="varIdentXOrthoEmail.Rout"

character(0)

Loading required package: nlme
[1] TRUE

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"

[1] "coef.varIdentX"       "coef<-.varIdentX"     "Initialize.varIdentX"
[4] "varIdent"             "varIdentX"

Classes 'varIdent', 'varFunc'  atomic [1:2] 0.0953 0.1310
  ..- attr(*, "groupNames")= chr [1:4] "8" "10" "14" "12"
  ..- attr(*, "fixed")= Named num 0.113
  .. ..- attr(*, "names")= chr "12"
  ..- attr(*, "formula")=Class 'formula' length 2 ~1 | age
  .. .. ..- attr(*, ".Environment")=<R_GlobalEnv> 
  ..- attr(*, "groups")= chr [1:108] "8" "10" "12" "14" ...
  ..- attr(*, "whichFix")= logi [1:3] FALSE FALSE  TRUE
  ..- attr(*, "weights")= Named num [1:108] 1.000 0.909 0.893 0.877 1.000 ...
  .. ..- attr(*, "names")= chr [1:108] "8" "10" "12" "14" ...
  ..- attr(*, "logLik")= num -9.17

[1] 0.09531018 0.13102826

8   10   14   12 
1.00 1.10 1.14 1.12

[1] 2

Error in `coef<-.varIdent`(`*tmp*`, value = c(11, 12)) : 
        Cannot change the length of the varIdent parameter after initialization

+                 weights  = vfi,
+                 data=Orthodont)
Error in `coef<-.varIdent`(`*tmp*`, value = c(0.095310179804325, 0.131028262406404 : 
        Cannot change the length of the varIdent parameter after initialization

Classes 'varIdentX', 'varFunc'  atomic [1:2] 0.0953 0.1310
  ..- attr(*, "groupNames")= chr [1:4] "8" "10" "14" "12"
  ..- attr(*, "fixed")= Named num 0.113
  .. ..- attr(*, "names")= chr "12"
  ..- attr(*, "formula")=Class 'formula' length 2 ~1 | age
  .. .. ..- attr(*, ".Environment")=<R_GlobalEnv> 
  ..- attr(*, "groups")= chr [1:108] "8" "10" "12" "14" ...
  ..- attr(*, "whichFix")= logi [1:3] FALSE FALSE  TRUE
  ..- attr(*, "weights")= Named num [1:108] 1.000 0.909 0.893 0.877 1.000 ...
  .. ..- attr(*, "names")= chr [1:108] "8" "10" "12" "14" ...
  ..- attr(*, "logLik")= num -9.17

[1] 0.09531018 0.13102826

8   10   14   12 
1.00 1.10 1.14 1.12

[1] 2

+                 weights  = vfi,
+                 data=Orthodont)

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

# source("C:/RBookX/QTools/RcodeFix/varIdentmethodsRevised.txt")
# After testing replace varIdentX -> varIdent
# corrections in varIdent


# Note: PR#9765 fixed
# Revised varIdent Function  with PR#9765 fixed


# === List of functions: ===
# 1. varIdent PR#9765 fixed
# 2. varIdentX PR#9765 fixed
# 3. coef<-.varIdentX  with corrections
# 4. Initialize.varIdentX  same as  Initialize.varIdent  
# 5. coef.varIdentX same as coef.varIdentX 

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
}
# End of revised varIdent function


varIdentX <- 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("varIdentX", "varFunc")
    value
}


"coef<-.varIdentX" <- function (object, ..., value) 
{
    if (!(is.null(grps <- getGroups(object)) || all(attr(object, 
        "whichFix")))) {
        value <- as.numeric(value)
        nGroups <- length(attr(object, "groupNames"))
        nFixed  <- sum(as.numeric(attr(object,"whichFix")))  # inserted line
        if (length(value) != nGroups - nFixed - 1) {# subtracted nFixed
            stop(paste("Cannot change the length of the varIdent", 
                "parameter after initialization"))
        }
        object[] <- value
        natPar <- coef(object, FALSE, allCoef = TRUE)
        attr(object, "logLik") <- sum(log(attr(object, "weights") <- 1/natPar[grps]))
    }
    object
}




### the same as varIdent

Initialize.varIdentX <- function (object, data, ...) 
{
    if (!is.null(form <- formula(object)) && !is.null(grpForm <- getGroupsFormula(form))) {
        if (length(coef(object)) > 0) {
            return(object)
        }
        strat <- attr(object, "groups") <- as.character(getGroups(data, 
            form, level = length(splitFormula(grpForm, sep = "*")), 
            sep = "*"))
        if (length((uStrat <- unique(strat))) == 1) {
            return(Initialize(varIdent(), data))
        }
        if (!is.null(fix <- attr(object, "fixed"))) {
            fixNames <- names(fix)
            if (any(is.na(match(fixNames, uStrat)))) {
                stop(paste("Fixed parameters names in varIdent", 
                  "must be a subset of groups names"))
            }
            uStratVar <- uStrat[is.na(match(uStrat, fixNames))]
            uStrat <- c(uStratVar, fixNames)
        }
        else {
            uStratVar <- uStrat
        }
        if ((nStratVar <- length(uStratVar)) == 0) {
            stop("Cannot fix variances in all groups")
        }
        if (nStratVar > 1) {
            if (length(object) <= 1) {
                oldAttr <- attributes(object)
                if (length(object) > 0) {
                  object <- rep(as.vector(object), nStratVar - 
                    1)
                }
                else {
                  object <- rep(0, nStratVar - 1)
                }
                attributes(object) <- oldAttr
                attr(object, "groupNames") <- uStrat
            }
            else {
                if (length(as.vector(object)) != (len <- (nStratVar - 
                  1))) {
                  stop(paste("Initial value for \"varIdent\" should be of length", 
                    len))
                }
                if (!is.null(stN <- attr(object, "groupNames"))) {
                  missStrat <- uStrat[is.na(match(uStrat, stN))]
                  if (length(missStrat) != 1) {
                    stop(paste("Names of starting value for \"varIdent\" object", 
                      "must contain all but one of the stratum levels"))
                  }
                  stN <- c(missStrat, stN)
                  if ((length(stN) != length(uStrat)) || any(sort(stN) != 
                    sort(uStrat))) {
                    stop("Nonexistent groups names for initial values in varIdent")
                  }
                  attr(object, "groupNames") <- stN
                }
                else {
                  attr(object, "groupNames") <- uStrat
                }
            }
        }
        else {
            oldAttr <- attributes(object)
            object <- numeric(0)
            attributes(object) <- oldAttr
            attr(object, "groupNames") <- uStrat
        }
        attr(object, "whichFix") <- !is.na(match(attr(object, 
            "groupNames")[-1], names(fix)))
        if (all(attr(object, "whichFix"))) {
            if (all(attr(object, "fixed") == 0)) {
                return(Initialize(varIdent(), data))
            }
            else {
                oldAttr <- attributes(object)
                object <- numeric(0)
                attributes(object) <- oldAttr
            }
        }
        attr(object, "logLik") <- sum(log(attr(object, "weights") <- 1/coef(object, 
            FALSE, allCoef = TRUE)[strat]))
        object
    }
    else {
        attr(object, "whichFix") <- TRUE
        NextMethod()
    }
}




coef.varIdentX <- function (object, unconstrained = TRUE, allCoef = FALSE, ...) 
{
    if (!is.null(getGroupsFormula(object)) && !is.null(wPar <- attr(object, 
        "whichFix"))) {
        if (unconstrained && !allCoef) {
            return(as.vector(object))
        }
        val <- double(length(wPar))
        if (any(wPar)) {
            val[wPar] <- attr(object, "fixed")
        }
        if (any(!wPar)) {
            val[!wPar] <- as.vector(object)
        }
        if (!unconstrained) {
            val <- c(1, exp(val))
            names(val) <- attr(object, "groupNames")
            if (!allCoef) {
                val <- val[c(FALSE, !attr(object, "whichFix"))]
            }
        }
        val
    }
    else {
        numeric(0)
    }
}

--------------040502080208060001050807--