jthaman / citools Goto Github PK
View Code? Open in Web Editor NEWAn R Package for Quick Uncertainty Intervals
License: GNU General Public License v3.0
An R Package for Quick Uncertainty Intervals
License: GNU General Public License v3.0
Binomial (vice Bernoulli) regression is supported by GLM and can be specified by using the "weights" option to give the number of trials per observation and the response variable to indicate the proportion of "success"s. Prediction intervals, quantile estimates, and probability estimates are sensible for such cases and could be estimated via simulation.
For certain methods like add_pi.glm and add_quantile.glm, a small parametric bootstrap is performed to generate the desired statistic. Presently if these two methods are used, a separate simulation is run for each method. If the number of replicates (nSims) and new data (tb) are not changed, it would be nice to store and reuse the simulation data.
The benefits of storing simulation data are twofold:
Though we could just handle the second point by setting a seed.
Dear John,
I am running GLMM models with family=Gamma
and I want to use add_ci()
to compute confidence intervals around predictions (many thanks for developing this by the way).
Unfortunately, add_ci()
often encounters the following error:
Error in quantile.default(newX[, i], ...) : missing values and NaN's not allowed if 'na.rm' is FALSE
often but not always; I can re-run add_ci()
and it might fail or not.
And it seems to be associated with this type of warning message:
In lme4::bootMer(fit, my_pred, nsim = nSims, type = "parametric", : some bootstrap runs failed (1/100)
Here is the code I’m running:
modfinal<-glmer(cumulative_duration ~ winter + month + s_cumul_dura_small + big_treat + (1 | id_first_feederpair), data=data_FRICOE, family = Gamma(link = "log"))
,
with cumulative_duration
being the time (in seconds) spent by the common chaffinch (Fringilla coelebs) on a feeder.
newdatpred <- expand.grid(winter=as.factor("4")
,month=as.factor("Jan")
,s_cumul_dura_small=0
,big_treat=as.factor(c("NONE","PICPIC"))
,id_first_feederpair=as.factor("0060564"))
add_ci(newdatpred,modfinal,alpha=0.01,response=T,includeRanef=F,type="boot",nSims=100,names=c("lowBCI","uppBCI"))
(Note: nSims set to 100 for speed reason, but planning on using 1000 for final results)
I hope I'm giving enough useful details on this. I'm happy to share more as needed.
I attach the data table data_FRICOE.txt if it can help reproduce this error (several runs of add_ci()
might be required given it does not always yield this error)
data_FRICOE.txt
I have not tested this but i am sure add_ci.glm doesn't work if it is supplied with an inverse-gaussian fit. I added a bootstrap procedure to mitigate this issue, but it's a definitely worth correcting. I would like to be able to calculate GLM CIs and PIs with complete coverage of the exp family.
If family=gaussian is used in add_pi, all the bounds to returned as zero vectors.
Since these prediction intervals can be solved parametrically, we should just pass to add_pi.lm if the Guassian family is used in a glm model.
At present, many GLM methods use t-distributions rather than z-distributions for GLMs where the dispersion parameter is fixed at 1. (See Section 2.6.2 of Applied generalized linear mixed models: Continuous and discrete data, 2010.) These methods should use a standard normal instead. See add_ci.glm, add_pi.glm, etc.
We lack an automated testing process that we can plug into the build process. We can implement some pretty good testing with the testthat package.
Warning message produced by lme4::lmer
fixed-effect model matrix is rank deficient so dropping 7 columns / coefficients.
Result with add_ci/add_pi
Error in X %*% vcovBetaHat : non-conformable arguments
Likely what is happening is that lmer is taking your specified model, noticing that the covariance matrix won't properly invert, and is thus dropping some columns so that it behaves better. When ciTools tries to do the same thing, it attempts to use the original, un-reduced matrix and thus produces the above error.
Was also wondering if it would be hard to also support mgcv::gam (generalized additive models) and scam (shape constrained additive models) models? The recipe to calculate confidence and prediction intervals for these should be quite similar to that for GLMs right?
Beta regression is a good alternative to binomial regression when the number of trials is unknown. I would like to support this model!
I'm working on this presently on the glmm2 branch. I have a basic parametric method for CIs, but the other methods may need to be implemented through a simulation
Trying add_ci
on a model with a formula involving I(...)
, i.e.
library(ciTools)
library(lme4)
data <- data.frame(x = c(1, 2, 3, 4), y = c(1, 3, 10, 15), gr = (1, 1, 2, 2))
mod <- lmer(y ~ I(x^2) + (1|gr), data)
add_ci(data, mod)
results in Error in x^2: non-numeric argument to binary operator
. Attempting add_ci(data, mod, type = "parametric")
results in Error in `[[<-.data.frame`(`*tmp*`, names[1], value = numeric(0)) : replacement has 0 rows, data has 4
. The code completes in other cases, but returns spurious confidence intervals and the warning longer object length is not a multiple of shorter object length
. At least in the parametric case, the problem appears to be in get_x_matrix_mermod
.
Trying add_ci
on a model with no fixed effects, i.e.
data <- data.frame(x = rep(seq(5), 2), y = c(rnorm(5), rnorm(5, mean = 0.5)), gr = c(rep(1, 5), rep(2, 5)))
mod <- lmer(y ~ (1 | gr), data)
data %>% add_ci(mod, type = "parametric")
results in Error in reformulate(attributes(terms(fit))$term.labels) : 'termlabels' must be a character vector of length at least one
. Again, the problem appears to be in get_x_matrix_mermod
.
If I run:
test <- add_ci(test, testModel, alpha = 0.1, names = c(“LCB”, “UCB”))
test$pred <- NA
test$LCB <- NA
test$UCB <- NA
THEN:
newDF <- add_ci(test, testModel, alpha = 0.1, names = c(“LCB”, “UCB”))
-or -
test <- add_ci(test, testModel, alpha = 0.1, names = c(“LCB”, “UCB”))
Both the newDF and test dataframes come up with the columns still set to NA. This occurs with values in the columns too, not only after setting the columns to NA. It does not always occur on the first attempt, but once it occurs it continues to occur.
Hello
I am new with ciTools but very interested in using it for post hoc analyses for my survreg models
I have a very naive question:
is there a way to omit NA ?
fit<- survreg(Surv(left, right, type="interval2") ~ drug, dist="weibull", data=mydata)
add_quantile(mydata, fit, p = 0.75, name = c("quant", "lwr", "upr"))
error
Error in add_quantile.survreg
Check df for missingness
thank you!
When lmer is given a rank-deficient model matrix, it may go ahead and fit the model anyways, dropping columns as necessary. It produces a message when it does this if you look at the model object, fit
or summary(fit)
. You can also look at the messages via summary(fit)$fitMsgs
.
When this takes place, add_ci.lmerMod
produces an error:
Error in X %*% vcovBetaHat : non-conformable arguments
I have a fix for this in the latest push of my master fork if you want to grab it. Alternatively, over-write get_x_matrix_mermod
from helper_functions.R
with the following:
get_x_matrix_mermod <- function(tb, fit){
##This function is necessary to avoid cases where the new data upon
##which CIs are generated does not contain all levels for one or more
##factors from the original data set. New and old data sets are
##appended, the model matrix is generated, and the function returns
##only the rows corresponding to the new data.
mm <- fit@frame
rv <- names(mm)[1]
mm[[rv]] <- as.numeric(mm[[rv]])
for(i in names(mm)){
if(is.factor(mm[[i]])) mm[[i]] <- as.character(mm[[i]])
}
X <- suppressWarnings(model.matrix(reformulate(attributes(terms(fit))$term.labels),
dplyr::bind_rows(mm, tb))[-(1:nrow(fit@frame)), ])
#remove columns from X that weren't used in the model fit
X[, colnames(X)[colnames(X) %in% colnames(model.matrix(fit))]]
}
Example:
newDF <- add_ci(test, testModel, alpha = 0.1, names = c(“LCB”, “UCB”))
If test does not have the columns pred, LCB, and UCB, no warning appears. (Expected outcome: no warning [working as expected])
If I rerun the exact same line and am overwriting newDF, no warning appears (Expected outcome: should get warning)
If I run the following:
test <- add_ci(test, testModel, alpha = 0.1, names = c(“LCB”, “UCB”))
newDF2 <- add_ci(test, testModel, alpha = 0.1, names = c(“LCB”, “UCB”))
I get a warning that I am overwriting data in newDF2, even though that df is created in that line and does not yet exist. (Expected outcome: no warning)
It seems that the warning only occurs when the dataframe defined in the function has the columns of the same name even though the warning should occur when the dataframe written to has columns of the same name.
I saw the add_quatile() function briefly. Maybe I'm wrong but the quantiles are done by simulation. Is it possible to do for the quasibinomial() family? I could make a pull request and someone analyses it.
Is this literature review mentioned in these slides publicly available anywhere?
When factor levels are missing from a model, no surprise, it cannot predict rows with that factor level. The issue comes from the use of predict()
in those models updated without all levels. The reprex below illustrates the issue.
set.seed(5)
library(ciTools)
#> ciTools version 0.6.1 (C) Institute for Defense Analyses
library(tidyverse)
my_data <-
data.frame(
counts=c(18,17,15,20,10,20,25,13,12),
outcome=gl(3,1,9),
treatment=gl(3,3)
)
my_model <- glm(counts ~ outcome + treatment, data=my_data)
# Works because it is not using `predict()`
my_new_data <-
my_data %>%
add_ci(fit=my_model)
# Fails because data subsetting is not contingent on input factor levels
my_new_data <-
my_data %>%
add_ci(fit=my_model, type="boot")
#> Error in model.frame.default(Terms, newdata, na.action = na.action, xlev = object$xlevels): factor treatment has new levels 2
Created on 2021-05-17 by the reprex package (v2.0.0)
This package is quite an accomplishment. I would use this frequently in analyses for publication, but I would need to be able to cite the sources for the interval methods you used
I've done some work in the past with CI's and PI's for non linear regression using the propagate
package. Would you be interested in receiving a PR to add nls
model support to this package?
Thank you for the wonderful package!
I am a new user, and I'm having an issue with a model that I fit in one environment and need to use ciTools in another environment. My general use case is illustrated by the reprex below:
set.seed(5)
library(ciTools)
#> ciTools version 0.6.1 (C) Institute for Defense Analyses
library(tidyverse)
my_data <-
data.frame(
counts=c(18,17,15,20,10,20,25,13,12),
outcome=rep(1:3, 3),
treatment=rep(1:3, each=3)
)
fit_my_data <- function(my_data_in_the_function) {
link <- "log"
glm.D93 <- glm(counts ~ outcome + treatment, data=my_data_in_the_function, family=gaussian(link=link))
glm.D93
}
my_model <- fit_my_data(my_data)
# Works because it is not using `update()`
my_new_data <-
my_data %>%
add_ci(fit=my_model)
# Fails because `link` is not set in the current environment (it is only set in
# the function environment)
my_new_data <-
my_data %>%
add_ci(fit=my_model, type="boot")
#> Error in gaussian(link = link): object 'link' not found
# Works because now the call doesn't rely on the link from the current
# environment; the link comes from the model family.
my_model_update <- my_model
my_model_update$call$family <- my_model$family
my_new_data <-
my_data %>%
add_ci(fit=my_model_update, type="boot")
#> Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
#> prediction from a rank-deficient fit may be misleading
Created on 2021-05-17 by the reprex package (v2.0.0)
Would it also be possible for the package to support confidence and prediction intervals for gls models?
I have a binomial GLMM of the form
fit = glmer(cbind(strain, TOTAL-strain) ~ 0+(1|REGION/obs) +
COUNTRY*SAMPLE_DATE,
family=binomial(logit), data=data)
Now I would like to produce confidence intervals, taking into account my fixed factors COUNTRY & SAMPLE_DATE and random factor intercept REGION, but ignoring the observation-level random effect (put in to take into account overdispersion).
I would also like to be able to extrapolate my model beyond the original sample dates in my dataset. Is there any way to achieve something like this? Or is the only option for me to set includeRanef = FALSE and then ignore the random intercepts per region? (which would not be ideal for me)
Very useful package. I have a minor suggestion: in add_pi
(and possibly other functions - I haven't checked), an error is thrown if tb
does not include a column for the response variable. The actual values in the column are ignored, but the column has to be present. I'd guess this is due to some internal code similar to:
X <- formula(fit) %>% model.matrix(data = tb)
used to get the design matrix for simulation-based prediction intervals. Using a few more steps should eliminate the need to include the response. Something like:
chr_formula <- formula(fit) %>% deparse() %>% strsplit(' ') %>% getElement(1)
X <- as.formula(chr_formula[-1]) %>% model.matrix(data = tb)
I noticed this specifically with a Poisson GLM. add_ci
did not require the response to be in tb
.
Is there the capability to create a confidence interval around a quantile?
Some users don't fit mixed models!
When I run the survreg with boot
option, I get the error.
> ## Define a data set.
> df <- survival::stanford2
> ## remove a covariate with missing values.
> df <- df[, 1:4]
> ## next, create the Surv object inside the survreg call:
> fit <- survival::survreg(survival::Surv(time, status) ~ age + I(age^2),
+ data = df, dist = "lognormal")
> add_pi(df, fit, alpha = 0.1, names = c("lwr", "upr"),method="boot")
Error in MASS::mvrnorm(nSims, beta.hat, vcov.hat) :
incompatible arguments
I know what is causing it, basically, the scale coefficient is not being picked up with Log(scale)
.
> vcov(fit)
(Intercept) age I(age^2) Log(scale)
(Intercept) 3.453920115 -0.1844604713 2.304116e-03 0.0052810250
age -0.184460471 0.0105711147 -1.379510e-04 0.0001468898
I(age^2) 0.002304116 -0.0001379510 1.854279e-06 -0.0000039567
Log(scale) 0.005281025 0.0001468898 -3.956700e-06 0.0049915276
> coef(fit)
(Intercept) age I(age^2)
3.344231092 0.226092785 -0.003504315
Gamma and Negative Binomial.
Not a top priority though
Hi. When I try to run add_ci with a glmerMod, I get the following error:
Error in [.data.frame
(fr, vars) : undefined columns selected
After tracing and running in debug, it seems the error is being thrown by a call to lme4::bootMer and the issue is in the newdata argument.
Would it be incorrect to change my_pred_fixed_glmer_linear
to drop the newdata argument? When I changed the my_pred
function in bootstrap_ci_glmermod
from
function (fit, lvl) { predict(fit, newdata = ciTools_data$tb_temp, re.form = NULL, family = fit@resp$family$family, type = "response") }
to
my_pred = function (fit, lvl) { predict(fit, re.form = NULL, family = fit@resp$family$family, type = "response") }
it seems to work. Since the data are contained inside the glmerMod object I don't see any need to have the newdata argument.
When viewed through the help file interface in RStudio, the equations don't render properly. They seem fine when accessed from the CRAN or other external sites, so I'm not sure what the problem is. Can you double-check this and let me know if it's only a problem on my end?
Currently, when we form prediction intervals via simulation, we take the alpha/2 and 1 - alpha/2 quantiles of the simulations to compute the interval estimate. Another method would be to provide the shortest intervals such that (1 - alpha)*M simulations lie inside. M is the total numbers of simulations set by nSims.
I'm finding that we cannot use add_ci.lmer for "big data". I tried an example from the mermod vignette with 200,000 observations and found that R couldn't put the new data frame into memory. Here's the example I tried:
## linear example
x_gen_mermod <- function(ng = 8, nw = 5){
n <- ng * nw
x2 <- runif(n)
group <- rep(as.character(1:ng), each = nw)
return(tibble::tibble(x2 = x2,
group = group))
}
mm_pipe <- function(tb, ...){
model.matrix(data = tb, ...)
}
get_validation_set <- function(tb, sigma, sigmaG, beta, includeRanef, groupIntercepts){
vm <- sample_n(tb, 5, replace = F)[rep(1:5, each = 100), ]
vf <- bind_rows(vm, tb) %>%
select(-group) %>%
mm_pipe(~.*.)
vf <- vf[1:500, ]
vGroups <- if(!includeRanef) rnorm(500, 0, sigmaG) else groupIntercepts[as.numeric(vm$group)]
vm[["y"]] <- vf %*% beta + vGroups + rnorm(500, mean = 0, sd = sigma)
vm
}
y_gen_mermod <- function(tb, sigma = 1, sigmaG = 1, delta = 1, includeRanef = FALSE, validationPoints = FALSE){
groupIntercepts <- rnorm(length(unique(tb$group)), 0, sigmaG)
tf <- tb %>%
dplyr::select(-group) %>%
mm_pipe(~.*.)
beta <- rep(delta, ncol(tf))
if(validationPoints) {
vm <- get_validation_set(tb, sigma, sigmaG, beta, includeRanef, groupIntercepts)
}
tb[["y"]] <- tf %*% beta + groupIntercepts[as.numeric(tb$group)] + rnorm(nrow(tb), mean = 0, sd = sigma)
tb[["truth"]] <- tf %*% beta + groupIntercepts[as.numeric(tb$group)] * (includeRanef)
if(validationPoints) return(list(tb = tb, vm = vm)) else return(tb)
}
tb <- x_gen_mermod(10, 20000) %>%
y_gen_mermod()
fit2 <- lmer(y ~ x2 + (1|group) , data = tb)
tb %>% add_ci(fit2, type = "parametric", includeRanef = TRUE, names = c("LCB", "UCB"))
Lmer works just fine on an example data set this large, but ciTools chokes and spits out
Error: cannot allocate vector of size 298.0 Gb
We need to re-examine how we are storing things in memory and see if we can do something more efficient. I'm not sure if this bug affects the other methods as well.
Run the example code in add_ci.lmerMod.R:
dat <- lme4::sleepstudy
# Fit a linear mixed model (random intercept model)
fit <- lme4::lmer(Reaction ~ Days + (1|Subject), data = lme4::sleepstudy)
# Get the fitted values for each observation in dat, and
# append CIs for those fitted values to dat
add_ci(dat, fit, alpha = 0.5)
I'm getting a warning message related to bind_rows:
add_ci(dat, fit, alpha = 0.5)
Warning in bind_rows_(x, .id) :
binding character and factor vector, coercing into character vector
# A tibble: 180 x 6
Reaction Days Subject pred LCB0.25 UCB0.75
<dbl> <dbl> <fctr> <dbl> <dbl> <dbl>
1 249.5600 0 308 292.1888 283.0011 301.3766
2 258.7047 1 308 302.6561 293.5979 311.7143
3 250.8006 2 308 313.1234 304.1636 322.0832
4 321.4398 3 308 323.5907 314.6971 332.4843
5 356.8519 4 308 334.0580 325.1976 342.9183
6 414.6901 5 308 344.5252 335.6649 353.3856
7 382.2038 6 308 354.9925 346.0989 363.8861
8 290.1486 7 308 365.4598 356.5000 374.4196
9 430.5853 8 308 375.9271 366.8689 384.9853
10 466.3535 9 308 386.3944 377.2066 395.5821
# ... with 170 more rows
This has to do with some code in the helper functions file. If it's harmless we can just wrap the code in suppressWarnings(), but I'm not convinced it's entirely harmless.
When we compute the SEs of random effects, we use
seRandom <- arm::se.ranef(fit)[[1]][,1] ##
which results in a scalar. This is fine if all groups are the same size. If the groups are not the same size, we should match the correct random effect SE to the group. In add_ci.glmer, we do this:
ranef_name <- names(fit@cnms)[1] ## just one random effect for now
seRandom <- arm::se.ranef(fit)[[1]][,1] ## get all random effect SEs
seRandom_vec <- rep(NA, length(tb[[ranef_name]]))
seRandom_df <- tibble( ## convert to data frame
group = names(seRandom),
seRandom = seRandom
)
names(seRandom_df)[names(seRandom_df) == 'group'] <- ranef_name ## get the group names so SEs can be matched
seRandom_vec <- dplyr::left_join(tb, seRandom_df, by = ranef_name)[["seRandom"]] ## match RE SEs to each obs in tb using ranef_name as the matching criterion
`
This code is ugly and I'm not proud of it, but I think we need to do something about this because the documentation does not state that add_ci.lmer can only be used if the group sizes are equal.
Currently we simulate from the model using a parametric bootstrap. We showed in the GLM vignette that this is valid if the model is "true", and the simulation is still relatively fast, but a parametric method would be more pleasing to me.
I'm not sure if there is any literature available on this.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.