Confidence Intervals Errors

[aneumann-science]

Dear Andrey,

First, thanks for developing lme4qtl, I have been looking for this kind of package for a while.

Second, I attempted to compute confidence intervals with the example data as well as with my own data and receive lots of warnings and errors. I used the code from the paper supplement and applied it to the quick start example:

library(lme4) # needed for `VarCorr` function
library(lme4qtl)

# load synthetic data set `dat40` distributed within `lme4qtl`
# - table of phenotypes `dat40`
# - the double kinship matrix `kin2`
data(dat40)

# (1) model continiuous trait `trait1`
mod <- relmatLmer(trait1 ~ AGE + SEX + (1|FAMID) + (1|ID), dat40, relmat = list(ID = kin2))

# `?lme4::profile`
prof <- profile(mod, which = "theta_", prof.scale = "varcov")

# `?lme4qtl::varpropProf`
prof_prop <- varpropProf(prof)

ci <- confint(prof_prop, level = 0.95)
ci

I get following warnings from profile():

Warning messages:
1: In sqrt(m) : NaNs produced
2: In sqrt(m) : NaNs produced
3: In sqrt(m) : NaNs produced
4: In sqrt(m) : NaNs produced
5: In sqrt(m) : NaNs produced
6: In sqrt(m) : NaNs produced
7: In zetafun(np, ns) : NAs detected in profiling
8: In nextpar(mat, cc, i, delta, lowcut, upcut) :
  Last two rows have identical or NA .zeta values: using minstep

and then when using varpropProf():

Error in na.fail.default(data.frame(x = as.numeric(obj1), y = as.numeric(obj2))) : 
  missing values in object

Is there any error in the script or is this bug? Are there alternative ways to obtain confidence intervals? I would appreciate it, if you could help with obtaining confidence intervals.

[variani]

Hi,

I am not sure what is happening. I run some scripts on the toy dataset (dat40) and the GAIT2 data, but I can't draw any conclusion. I observe the following:

• The GAIT2 example is reproducible (good news).
• I can't perform profiling on dat40 when fitting data with lme4qtl. Though, I can profile and further estimate CI with lme4 (when dropping the (1|ID) term).

At the same time, I find this lme4 issue that shows this kind of errors occurs for models fitted by lme4.

In conclusion, it is hard to say whether the errors come due to lme4qtl functionality on the top of lme4; or "bad" data examples and numerical problems that might happen for lme4 (consequently, for lme4qtl).

I will keep exploring this issue, but I actually don't have much time. Feel free to do some experiments and I will be happy to follow them.

Best,
Andrey

[variani]

Another quick comment:

library(lme4)
?confint.merMod

There you will see that the Wald method gives you CI estimates, but only for fixed effects:

fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
fm1W <- confint(fm1, method = "Wald") # very fast, but ....

[variani]

I added a relevant example in REDME.Rmd + examples that mostly fail in demo/profile.R.

m2 <- relmatLmer(trait2 ~ (1|ID), dat40, relmat = list(ID = kin2)) 
VarProp(m2)
prof <- profile(m2)
prof_prop <- varpropProf(prof)
confint(prof_prop)

Going back to the first post here:

• warnings with NA come from singular fits when profiling (I guess fitting problems are not expected at this stage);
• error comes from fitting splines (calling interpSpline); see, for example, a small function lme4:::logProf

For me it is still an open question on how to perform profiling. Any suggestions and tips are welcome.

[variani]

I tried the delta method in demo/se.R and added some links at the top of script.

data(dat40, package = "lme4qtl")
m <- relmatLmer(trait1 ~ AGE + SEX + (1|FAMID) + (1|ID), dat40, relmat = list(ID = kin2))

# hessian
f <- update(m, devFunOnly = TRUE)

th <- getME(m, "theta")

library(numDeriv)
h <- numDeriv::hessian(f, th)

library(MASS)
cov <- MASS::ginv(h)

# delta method
nth <- length(th)
terms <- paste0("x", seq(nth), "*", "x", seq(nth)) # "x1*x1" "x2*x2"
terms <- c(terms, "1") # # "x1*x1" "x2*x2", "1"

total <- paste(terms, collapse = " + ") # "x1*x1 + x2*x2 + 1"

forms <- lapply(seq(nth), function(i) 
 paste("~", terms[i], "/ (", total, ")"))

library(msm)  
ses <- sapply(forms, function(f)
  msm::deltamethod(as.formula(f), th, cov, ses = TRUE))

I get the following results.

    grp       est         se
1    ID 0.8954191 0.05862602
2 FAMID 0.0000000 0.00000000