[Rd] Bug in coef<-.varIdent method (nlme package) (PR#9831)
agalecki at umich.edu
agalecki at umich.edu
Tue Aug 7 21:44:24 CEST 2007
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"
> ls()
character(0)
> require(nlme)
Loading required package: nlme
[1] TRUE
> 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"
>
>
> ## 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()
[1] "coef.varIdentX" "coef<-.varIdentX" "Initialize.varIdentX"
[4] "varIdent" "varIdentX"
>
>
>
> #### 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)
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
> coef(vfi)
[1] 0.09531018 0.13102826
> coef(vfi, unconstrained = FALSE, allCoef = TRUE)
8 10 14 12
1.00 1.10 1.14 1.12
> vfiCopy <- vfi # copy of an initialized object
>
> #### Chunk 2: Bug in coef<-.varIdent illustrated
>
> length(vfiCopy) # length is 2
[1] 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
> # 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
>
> #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
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
> coef(vfi)
[1] 0.09531018 0.13102826
> coef(vfi, unconstrained = FALSE, allCoef = TRUE)
8 10 14 12
1.00 1.10 1.14 1.12
> vfiCopy <- vfi # copy of an initialized object
>
> #### Chunk 2A: Bug in coef<-.varIdent *** corrected ***
>
> length(vfiCopy) # length is 2
[1] 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="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--
More information about the R-devel
mailing list