[R] error using user-defined link function with mixed models (LMER)
Jessi Brown
jessilbrown at gmail.com
Sun Feb 11 23:38:52 CET 2007
Hmmm.. I had actually at one point tried the newer (>R 2.4.0) method
of specifying only a logexp link function and using it with the
binomial family, and once again, it seems to work as expected for
glm's, and not so much with LMER. The error message when using
summary() was the same. Here's output from a non-nest exposure
logistic regression:
> apfa.lmer.b.b.1<-lmer(Success~MeanAge+I(MeanAge^2)+(1|Territory), data=apfa4, family=binomial(link="logit"), method="Laplace", control=list(msVerbose=1))
relative tolerance set to 3.18713398786316e-05
0 313.762: 3.69488 -0.100835 0.00165992 0.145015
1 313.758: 3.69488 -0.100834 0.00166623 0.145015
2 313.758: 3.69488 -0.100834 0.00166623 0.145015
3 313.758: 3.69488 -0.100834 0.00166623 0.145015
> str(apfa.lmer.b.b.1)
Formal class 'glmer' [package "lme4"] with 33 slots
..@ family :List of 11
.. ..$ family : chr "binomial"
.. ..$ link : chr "logit"
.. ..$ linkfun :function (mu)
.. ..$ linkinv :function (eta)
.. ..$ variance :function (mu)
.. ..$ dev.resids:function (y, mu, wt)
.. ..$ aic :function (y, n, mu, wt, dev)
.. ..$ mu.eta :function (eta)
.. ..$ initialize: expression({ if (NCOL(y) == 1) { if
(is.factor(y)) y <- y != levels(y)[1] n <-
rep.int(1, nobs) if (any(y < 0 | y > 1)) stop("y
values must be 0 <= y <= 1") mustart <- (weights * y +
0.5)/(weights + 1) m <- weights * y if (any(abs(m -
round(m)) > 0.001)) warning("non-integer #successes in a
binomial glm!") } else if (NCOL(y) == 2) { if
(any(abs(y - round(y)) > 0.001)) warning("non-integer
counts in a binomial glm!") n <- y[, 1] + y[, 2] y <-
ifelse(n == 0, 0, y[, 1]/n) weights <- weights * n
mustart <- (n * y + 0.5)/(n + 1) } else stop("for the binomial
family, y must be a vector of 0 and 1's\n", "or a 2 column
matrix where col 1 is no. successes and col 2 is no. failures") })
.. ..$ validmu :function (mu)
.. ..$ valideta :function (eta)
.. ..- attr(*, "class")= chr "family"
..@ weights : Named num [1:662] 1 1 1 1 1 1 1 1 1 1 ...
.. ..- attr(*, "names")= chr [1:662] "1" "2" "3" "4" ...
..@ frame :'data.frame': 662 obs. of 4 variables:
.. ..$ Success : int [1:662] 1 1 1 1 1 1 1 1 1 1 ...
.. ..$ MeanAge : num [1:662] 5.5 12.5 19.5 26.5 33.5 40.5 47.5
54.5 61.5 66.5 ...
.. ..$ I(MeanAge^2):Class 'AsIs' num [1:662] 30.2 156.2 380.2
702.2 1122.2 ...
.. ..$ Territory : Factor w/ 36 levels "9 MILE","B-1",..: 1 1 1 1
1 1 1 1 1 1 ...
.. ..- attr(*, "terms")=Classes 'terms', 'formula' length 3 Success
~ MeanAge + I(MeanAge^2) + (1 + Territory)
.. .. .. ..- attr(*, "variables")= language list(Success, MeanAge,
I(MeanAge^2), Territory)
.. .. .. ..- attr(*, "factors")= int [1:4, 1:3] 0 1 0 0 0 0 1 0 0 0 ...
.. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. ..$ : chr [1:4] "Success" "MeanAge" "I(MeanAge^2)" "Territory"
.. .. .. .. .. ..$ : chr [1:3] "MeanAge" "I(MeanAge^2)" "Territory"
.. .. .. ..- attr(*, "term.labels")= chr [1:3] "MeanAge"
"I(MeanAge^2)" "Territory"
.. .. .. ..- attr(*, "order")= int [1:3] 1 1 1
.. .. .. ..- attr(*, "intercept")= int 1
.. .. .. ..- attr(*, "response")= int 1
.. .. .. ..- attr(*, ".Environment")=length 26 <environment>
.. .. .. ..- attr(*, "predvars")= language list(Success, MeanAge,
I(MeanAge^2), Territory)
.. .. .. ..- attr(*, "dataClasses")= Named chr [1:4] "numeric"
"numeric" "other" "factor"
.. .. .. .. ..- attr(*, "names")= chr [1:4] "Success" "MeanAge"
"I(MeanAge^2)" "Territory"
..@ call : language lmer(formula = Success ~ MeanAge +
I(MeanAge^2) + (1 | Territory), data = apfa4, family =
binomial(link = "logit"), method = "Laplace", ...
..@ terms :Classes 'terms', 'formula' length 3 Success ~ MeanAge +
I(MeanAge^2)
.. .. ..- attr(*, "variables")= language list(Success, MeanAge, I(MeanAge^2))
.. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$ : chr [1:3] "Success" "MeanAge" "I(MeanAge^2)"
.. .. .. .. ..$ : chr [1:2] "MeanAge" "I(MeanAge^2)"
.. .. ..- attr(*, "term.labels")= chr [1:2] "MeanAge" "I(MeanAge^2)"
.. .. ..- attr(*, "order")= int [1:2] 1 1
.. .. ..- attr(*, "intercept")= int 1
.. .. ..- attr(*, "response")= int 1
.. .. ..- attr(*, ".Environment")=length 26 <environment>
.. .. ..- attr(*, "predvars")= language list(Success, MeanAge, I(MeanAge^2))
.. .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "numeric" "other"
.. .. .. ..- attr(*, "names")= chr [1:3] "Success" "MeanAge" "I(MeanAge^2)"
..@ flist :List of 1
.. ..$ Territory: Factor w/ 36 levels "9 MILE","B-1",..: 1 1 1 1 1 1
1 1 1 1 ...
..@ Zt :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. ..@ i : int [1:662] 0 0 0 0 0 0 0 0 0 0 ...
.. .. ..@ p : int [1:663] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ Dim : int [1:2] 36 662
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : NULL
.. .. .. ..$ : NULL
.. .. ..@ x : num [1:662] 1 1 1 1 1 1 1 1 1 1 ...
.. .. ..@ factors : list()
..@ X : num [1:662, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:662] "1" "2" "3" "4" ...
.. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. ..- attr(*, "assign")= int [1:3] 0 1 2
..@ y : num [1:662] 1 1 1 1 1 1 1 1 1 1 ...
..@ wts : num [1:662] 0.176 0.221 0.256 0.275 0.275 ...
..@ wrkres : num [1:662] 4.44 3.97 3.66 3.50 3.50 ...
..@ cnames :List of 2
.. ..$ Territory: chr "(Intercept)"
.. ..$ .fixed : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
..@ nc : Named int 1
.. ..- attr(*, "names")= chr "Territory"
..@ Gp : int [1:2] 0 36
..@ XtX :Formal class 'dpoMatrix' [package "Matrix"] with 5 slots
.. .. ..@ x : num [1:9] 39.7 0.0 0.0 1179.4 43661.6 ...
.. .. ..@ Dim : int [1:2] 3 3
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. .. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. .. ..@ uplo : chr "U"
.. .. ..@ factors : list()
..@ ZtZ :Formal class 'dsCMatrix' [package "Matrix"] with 7 slots
.. .. ..@ i : int [1:36] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ p : int [1:37] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ Dim : int [1:2] 36 36
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : NULL
.. .. .. ..$ : NULL
.. .. ..@ x : num [1:36] 1.416 0.534 1.069 0.529 1.376 ...
.. .. ..@ uplo : chr "U"
.. .. ..@ factors : list()
..@ ZtX :Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
.. .. ..@ x : num [1:108] 1.416 0.534 1.069 0.529 1.376 ...
.. .. ..@ Dim : int [1:2] 36 3
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : NULL
.. .. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. .. ..@ factors : list()
..@ Zty : num [1:36] 5.47 2.00 1.77 1.98 4.08 ...
..@ Xty : num [1:3] 100 3026 117899
..@ Omega :List of 1
.. ..$ Territory:Formal class 'dpoMatrix' [package "Matrix"] with 5 slots
.. .. .. ..@ x : num 6.9
.. .. .. ..@ Dim : int [1:2] 1 1
.. .. .. ..@ Dimnames:List of 2
.. .. .. .. ..$ : chr "(Intercept)"
.. .. .. .. ..$ : chr "(Intercept)"
.. .. .. ..@ uplo : chr "U"
.. .. .. ..@ factors :List of 1
.. .. .. .. ..$ Cholesky:Formal class 'Cholesky' [package "Matrix"]
with 5 slots
.. .. .. .. .. .. ..@ x : num 2.63
.. .. .. .. .. .. ..@ Dim : int [1:2] 1 1
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ uplo : chr "U"
.. .. .. .. .. .. ..@ diag : chr "N"
..@ L :Formal class 'dCHMsuper' [package "Matrix"] with 9 slots
.. .. ..@ x : num [1:36] 2.88 2.73 2.82 2.72 2.88 ...
.. .. ..@ super : int [1:37] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ pi : int [1:37] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ px : int [1:37] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ s : int [1:36] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ colcount: int [1:36] 1 1 1 1 1 1 1 1 1 1 ...
.. .. ..@ perm : int [1:36] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ type : int [1:6] 0 1 1 1 1 1
.. .. ..@ Dim : int [1:2] 36 36
..@ RZX :Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
.. .. ..@ x : num [1:108] 0.491 0.196 0.379 0.194 0.479 ...
.. .. ..@ Dim : int [1:2] 36 3
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : NULL
.. .. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. .. ..@ factors : list()
..@ RXX :Formal class 'dtrMatrix' [package "Matrix"] with 5 slots
.. .. ..@ x : num [1:9] 5.8 0.0 0.0 172.2 92.7 ...
.. .. ..@ Dim : int [1:2] 3 3
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. .. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. .. ..@ uplo : chr "U"
.. .. ..@ diag : chr "N"
..@ rZy : num [1:36] 1.896 0.733 0.628 0.726 1.419 ...
..@ rXy : num [1:3] 14.54 0.49 2.66
..@ devComp : Named num [1:8] 662.00 3.00 872.61 6.41 74.81 ...
.. ..- attr(*, "names")= chr [1:8] "n" "p" "yty" "logryy2" ...
..@ deviance: Named num [1:2] 314 NA
.. ..- attr(*, "names")= chr [1:2] "ML" "REML"
..@ fixef : num [1:3] 3.69488 -0.10083 0.00167
..@ ranef : num [1:36] 0.2191 0.0834 -0.1204 0.0827 0.0703 ...
..@ RZXinv :Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
.. .. ..@ x : num [1:108] 2.12e-313 -1.57e-154 0.00e+00
6.88e-310 9.60e-206 ...
.. .. ..@ Dim : int [1:2] 36 3
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : NULL
.. .. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. .. ..@ factors : list()
..@ bVar :List of 1
.. ..$ Territory: num [1, 1, 1:36] 1.84e-299 3.83e-301 1.60e-213
1.34e-300 2.47e-300 ...
..@ gradComp:List of 1
.. ..$ Territory: num [1, 1, 1:4] 2.74e-208 2.74e-208 8.39e-204 2.74e-208
..@ status : Named int [1:3] 2 0 2
.. ..- attr(*, "names")= chr [1:3] "stage" "REML" "glmm"
Now, here's the equivalent from using the logexp link function as
copied from the "family" help file example:
> apfa.lmer.b.1<-lmer(Success~MeanAge+I(MeanAge^2)+(1|Territory), data=apfa4, family=binomial(link="logexp"(apfa4$Days)), method="Laplace", control=list(msVerbose=1))
relative tolerance set to 3.49161751721116e-05
0 286.400: 5.55948 -0.0908251 0.00149257 0.145015
1 286.400: 5.55948 -0.0908251 0.00149257 0.145015
> str(apfa.lmer.b.1)
Formal class 'glmer' [package "lme4"] with 33 slots
..@ family :List of 11
.. ..$ family : chr "binomial"
.. ..$ link : chr [1:662] "logexp(7)" "logexp(7)" "logexp(7)"
"logexp(7)" ...
.. ..$ linkfun :function (mu)
.. .. ..- attr(*, "source")= chr "function(mu) qlogis(mu^(1/days))"
.. ..$ linkinv :function (eta)
.. .. ..- attr(*, "source")= chr "function(eta) plogis(eta)^days"
.. ..$ variance :function (mu)
.. ..$ dev.resids:function (y, mu, wt)
.. ..$ aic :function (y, n, mu, wt, dev)
.. ..$ mu.eta :function (eta)
.. .. ..- attr(*, "source")= chr [1:2] "function(eta) days *
plogis(eta)^(days-1) *" ...
.. ..$ initialize: expression({ if (NCOL(y) == 1) { if
(is.factor(y)) y <- y != levels(y)[1] n <-
rep.int(1, nobs) if (any(y < 0 | y > 1)) stop("y
values must be 0 <= y <= 1") mustart <- (weights * y +
0.5)/(weights + 1) m <- weights * y if (any(abs(m -
round(m)) > 0.001)) warning("non-integer #successes in a
binomial glm!") } else if (NCOL(y) == 2) { if
(any(abs(y - round(y)) > 0.001)) warning("non-integer
counts in a binomial glm!") n <- y[, 1] + y[, 2] y <-
ifelse(n == 0, 0, y[, 1]/n) weights <- weights * n
mustart <- (n * y + 0.5)/(n + 1) } else stop("for the binomial
family, y must be a vector of 0 and 1's\n", "or a 2 column
matrix where col 1 is no. successes and col 2 is no. failures") })
.. ..$ validmu :function (mu)
.. ..$ valideta :function (eta)
.. .. ..- attr(*, "source")= chr "function(eta) TRUE"
.. ..- attr(*, "class")= chr "family"
..@ weights : Named num [1:662] 1 1 1 1 1 1 1 1 1 1 ...
.. ..- attr(*, "names")= chr [1:662] "1" "2" "3" "4" ...
..@ frame :'data.frame': 662 obs. of 4 variables:
.. ..$ Success : int [1:662] 1 1 1 1 1 1 1 1 1 1 ...
.. ..$ MeanAge : num [1:662] 5.5 12.5 19.5 26.5 33.5 40.5 47.5
54.5 61.5 66.5 ...
.. ..$ I(MeanAge^2):Class 'AsIs' num [1:662] 30.2 156.2 380.2
702.2 1122.2 ...
.. ..$ Territory : Factor w/ 36 levels "9 MILE","B-1",..: 1 1 1 1
1 1 1 1 1 1 ...
.. ..- attr(*, "terms")=Classes 'terms', 'formula' length 3 Success
~ MeanAge + I(MeanAge^2) + (1 + Territory)
.. .. .. ..- attr(*, "variables")= language list(Success, MeanAge,
I(MeanAge^2), Territory)
.. .. .. ..- attr(*, "factors")= int [1:4, 1:3] 0 1 0 0 0 0 1 0 0 0 ...
.. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. ..$ : chr [1:4] "Success" "MeanAge" "I(MeanAge^2)" "Territory"
.. .. .. .. .. ..$ : chr [1:3] "MeanAge" "I(MeanAge^2)" "Territory"
.. .. .. ..- attr(*, "term.labels")= chr [1:3] "MeanAge"
"I(MeanAge^2)" "Territory"
.. .. .. ..- attr(*, "order")= int [1:3] 1 1 1
.. .. .. ..- attr(*, "intercept")= int 1
.. .. .. ..- attr(*, "response")= int 1
.. .. .. ..- attr(*, ".Environment")=length 25 <environment>
.. .. .. ..- attr(*, "predvars")= language list(Success, MeanAge,
I(MeanAge^2), Territory)
.. .. .. ..- attr(*, "dataClasses")= Named chr [1:4] "numeric"
"numeric" "other" "factor"
.. .. .. .. ..- attr(*, "names")= chr [1:4] "Success" "MeanAge"
"I(MeanAge^2)" "Territory"
..@ call : language lmer(formula = Success ~ MeanAge +
I(MeanAge^2) + (1 | Territory), data = apfa4, family =
binomial(link = logexp(apfa4$Days)), ...
..@ terms :Classes 'terms', 'formula' length 3 Success ~ MeanAge +
I(MeanAge^2)
.. .. ..- attr(*, "variables")= language list(Success, MeanAge, I(MeanAge^2))
.. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$ : chr [1:3] "Success" "MeanAge" "I(MeanAge^2)"
.. .. .. .. ..$ : chr [1:2] "MeanAge" "I(MeanAge^2)"
.. .. ..- attr(*, "term.labels")= chr [1:2] "MeanAge" "I(MeanAge^2)"
.. .. ..- attr(*, "order")= int [1:2] 1 1
.. .. ..- attr(*, "intercept")= int 1
.. .. ..- attr(*, "response")= int 1
.. .. ..- attr(*, ".Environment")=length 25 <environment>
.. .. ..- attr(*, "predvars")= language list(Success, MeanAge, I(MeanAge^2))
.. .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "numeric" "other"
.. .. .. ..- attr(*, "names")= chr [1:3] "Success" "MeanAge" "I(MeanAge^2)"
..@ flist :List of 1
.. ..$ Territory: Factor w/ 36 levels "9 MILE","B-1",..: 1 1 1 1 1 1
1 1 1 1 ...
..@ Zt :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. ..@ i : int [1:662] 0 0 0 0 0 0 0 0 0 0 ...
.. .. ..@ p : int [1:663] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ Dim : int [1:2] 36 662
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : NULL
.. .. .. ..$ : NULL
.. .. ..@ x : num [1:662] 1 1 1 1 1 1 1 1 1 1 ...
.. .. ..@ factors : list()
..@ X : num [1:662, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:662] "1" "2" "3" "4" ...
.. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. ..- attr(*, "assign")= int [1:3] 0 1 2
..@ y : num [1:662] 1 1 1 1 1 1 1 1 1 1 ...
..@ wts : num [1:662] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
..@ wrkres : num [1:662] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
..@ cnames :List of 2
.. ..$ Territory: chr "(Intercept)"
.. ..$ .fixed : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
..@ nc : Named int 1
.. ..- attr(*, "names")= chr "Territory"
..@ Gp : int [1:2] 0 36
..@ XtX :Formal class 'dpoMatrix' [package "Matrix"] with 5 slots
.. .. ..@ x : num [1:9] NaN 0 0 NaN NaN 0 NaN NaN NaN
.. .. ..@ Dim : int [1:2] 3 3
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. .. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. .. ..@ uplo : chr "U"
.. .. ..@ factors : list()
..@ ZtZ :Formal class 'dsCMatrix' [package "Matrix"] with 7 slots
.. .. ..@ i : int [1:36] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ p : int [1:37] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ Dim : int [1:2] 36 36
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : NULL
.. .. .. ..$ : NULL
.. .. ..@ x : num [1:36] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
.. .. ..@ uplo : chr "U"
.. .. ..@ factors : list()
..@ ZtX :Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
.. .. ..@ x : num [1:108] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
.. .. ..@ Dim : int [1:2] 36 3
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : NULL
.. .. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. .. ..@ factors : list()
..@ Zty : num [1:36] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
..@ Xty : num [1:3] NaN NaN NaN
..@ Omega :List of 1
.. ..$ Territory:Formal class 'dpoMatrix' [package "Matrix"] with 5 slots
.. .. .. ..@ x : num 6.9
.. .. .. ..@ Dim : int [1:2] 1 1
.. .. .. ..@ Dimnames:List of 2
.. .. .. .. ..$ : chr "(Intercept)"
.. .. .. .. ..$ : chr "(Intercept)"
.. .. .. ..@ uplo : chr "U"
.. .. .. ..@ factors :List of 1
.. .. .. .. ..$ Cholesky:Formal class 'Cholesky' [package "Matrix"]
with 5 slots
.. .. .. .. .. .. ..@ x : num 2.63
.. .. .. .. .. .. ..@ Dim : int [1:2] 1 1
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ uplo : chr "U"
.. .. .. .. .. .. ..@ diag : chr "N"
..@ L :Formal class 'dCHMsuper' [package "Matrix"] with 9 slots
.. .. ..@ x : num [1:36] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
.. .. ..@ super : int [1:37] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ pi : int [1:37] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ px : int [1:37] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ s : int [1:36] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ colcount: int [1:36] 1 1 1 1 1 1 1 1 1 1 ...
.. .. ..@ perm : int [1:36] 0 1 2 3 4 5 6 7 8 9 ...
.. .. ..@ type : int [1:6] 0 1 1 1 1 1
.. .. ..@ Dim : int [1:2] 36 36
..@ RZX :Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
.. .. ..@ x : num [1:108] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
.. .. ..@ Dim : int [1:2] 36 3
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : NULL
.. .. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. .. ..@ factors : list()
..@ RXX :Formal class 'dtrMatrix' [package "Matrix"] with 5 slots
.. .. ..@ x : num [1:9] NaN 0 0 NaN NaN 0 NaN NaN NaN
.. .. ..@ Dim : int [1:2] 3 3
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. .. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. .. ..@ uplo : chr "U"
.. .. ..@ diag : chr "N"
..@ rZy : num [1:36] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
..@ rXy : num [1:3] NaN NaN NaN
..@ devComp : Named num [1:8] 662 3 NaN NaN NaN ...
.. ..- attr(*, "names")= chr [1:8] "n" "p" "yty" "logryy2" ...
..@ deviance: Named num [1:2] NaN NaN
.. ..- attr(*, "names")= chr [1:2] "ML" "REML"
..@ fixef : num [1:3] 5.55948 -0.09083 0.00149
..@ ranef : num [1:36] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
..@ RZXinv :Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
.. .. ..@ x : num [1:108] 2.12e-313 -1.57e-154 0.00e+00
7.53e-310 1.74e-206 ...
.. .. ..@ Dim : int [1:2] 36 3
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : NULL
.. .. .. ..$ : chr [1:3] "(Intercept)" "MeanAge" "I(MeanAge^2)"
.. .. ..@ factors : list()
..@ bVar :List of 1
.. ..$ Territory: num [1, 1, 1:36] 0.00e+00 0.00e+00 0.00e+00
1.76e-315 1.76e-315 ...
..@ gradComp:List of 1
.. ..$ Territory: num [1, 1, 1:4] 2.12e-314 4.24e-314 NaN NaN
..@ status : Named int [1:3] 2 0 2
.. ..- attr(*, "names")= chr [1:3] "stage" "REML" "glmm"
It's obvious that there are a lot of NaN's throughout the second
example, which I imagine are indicative of something not going as
planned, but I'm afraid I'm not familiar enough with the mathematical
reality behind these models to know what exactly is going on.
To review, here's the logexp link function suggested in the ?family file:
logexp <- function(days = 1)
{
linkfun <- function(mu) qlogis(mu^(1/days))
linkinv <- function(eta) plogis(eta)^days
mu.eta <- function(eta) days * plogis(eta)^(days-1) *
.Call("logit_mu_eta", eta, PACKAGE = "stats")
valideta <- function(eta) TRUE
link <- paste("logexp(", days, ")", sep="")
structure(list(linkfun = linkfun, linkinv = linkinv,
mu.eta = mu.eta, valideta = valideta, name = link),
class = "link-glm")
}
Any ideas? By the way, I am extremely grateful for all of the help so
far - I know this sort of thing must be tedious for you. Thank you for
putting up with a graduate student whose strength is in finding and
following birds, not programming or even thoroughly understanding the
complex underpinnings of the sorts of models I blithely attempt to
run!
cheers, Jessi Brown
More information about the R-help
mailing list