alexpkeil1 / qgcomp Goto Github PK

The model confidence and pointwise bounds shown by plot.qgcomp are too wide due to bugs in calculation. These calculation errors are fixed in version 2.1.2, which will be shortly available on CRAN, but can be installed now via

devtools::install_github("alexpkeil1/qgcomp")

This install also includes the functions modelbound.boot and pointwisebound.boot to estimate standard errors/confidence intervals for the qgcomp regression fits.

The current versions of pointwisebound.(no)boot functions appear to have reverted to an old version that is for linear models only (affects qgcomp v 2.5 to 2.5.7). I am working on bringing back the old version. This affects plots for ".boot" objects as well as pointwise effect measure intervals, but does not otherwise impact results.

Hello, I am wondering if there's a way to remove the intercept when using the qgcomp.boot function? I've tried adding a zero or -1 to the model formula, which works when using the R function 'glm,' but I still see an intercept term present in the model results.

You may see an error like this:

Error in `$<-.data.frame`(`*tmp*`, "weights", value = c(1, 1, 1, 1, 1, : replacement has 332 rows, data has 1664

This appears to be related to having missing values in the data.

If using .fixbreaks = TRUE (default) in qgcomp.partials, a bug presents where the quantile definitions for one variable are imposed on other variables in an arbitrary fashion. If using v2.9.1 or v2.9.0, it is recommended to set .fixbreaks=FALSE (the old default behavior) until you can install v2.10.0, which is currently in testing.

Would be nice to have a longitudinal qgcomp with some of the following possibilities:

• Effects of increasing all prior exposures by one quantile
• Quantiles for time-varying exposures that are based on single, common set of exposures OR quantiles that vary over time
• Prompt effects with mixed effects models
• Incorporation of arbitrary longitudinal models

From an email:

In mice.impute.leftcenslognorm (), the examples given can not run normally, always an Error is as follows: Error: No tidy method recognized for this list.

qgcomp v2.8.0

This appears to be due to some change in another package. Will investigate fixes.

Hi Alex,
Thanks for working on this fabulous package enabling the assessment of the effect of the exposure mixture. I am wondering if there is a way to extract the exposure index from the fit object? Such index should be useful for downstream analyses, such as mutually adjustment for multiple groups of exposure mixture, just in the similar way of group WQS.

Thanks.

Kai

There is a bug in ggplot2 code that prevents qgcomp plots (bar plots with weights) from printing correctly under some pretty specific circumstances (e.g. bars on one side of the plot are much wider than the bars on the other side of the plot).

I can't fix this on my end, but will update this issue once I have a chance to put in a bug report for ggplot2.

There are currently no fixes other than plotting the weights manually, but I can't provide a general solution to that.

Below is a reproducible example that uses some of the underlying qgcomp plotting code, which is just a wrapper around ggplot2 functions:

library(ggplot2)

pw = c(x1=.9, x3=0.1)
nw = c(x2=0.5, z=0.5)
nms = names(sort(c(pw,nw)))

ggplot() + 
  stat_identity(aes(x=v, y=w), position = "identity", geom="bar", 
                data=data.frame(w=pw, v=names(pw))) + 
  scale_x_discrete(name="Variable", limits=nms, breaks=nms, labels=nms, drop=FALSE) 
ggplot() + 
  stat_identity(aes(x=v, y=w), position = "identity", geom="bar", 
                data=data.frame(w=nw, v=names(nw))) + 
  scale_x_discrete(name="Variable", limits=nms, breaks=nms, labels=nms, drop=FALSE)

How should the NHANES weight be handled in the fitting of the qgcomp model? Is there an appropriate parameter that can be used?

My outcome is binary and I fitting a model using qgcomp.boot, rr = TRUE, B=5000, q=4

The summary estimates from qgcomp.boot shows

                            Estimate           Std.Error
         Intercept         -4.109              0.3428
        psi1               **0.06719**        0.1635

The summary from modelbound.boot(...) shows

                       quantile      midpoint         linepred         r
               0      0                 0.125             -4.109        **0.0456**
               1      1                 0.375                               
               2      2                 0.625
               3      3                 0.875

the vignette https://cran.r-project.org/web/packages/qgcomp/qgcomp.pdf, pg 14 says the model estimate of the regression line E(Y|X) are expected to include the true value of E(Y|X) in 95% of studies. My y is 1 or 0, how does the estimate r= 0.0456 translate to E(Y|X) for my binary outcome ? Please help. Thanks.

In version 1.3.0, a change in the print function lead to improper reporting of confidence intervals. The correct confidence intervals can be obtained by object$ci.coef until version 2.0 shows up on CRAN, or the master branch can be installed via github: via

devtools::install_github("alexpkeil1/qgcomp")

This isn't a "bug" just a question of how estimating a competing risks model might be accomplished.

It seems that the Surv package in R doesn't do competing risks. For that, you have to use the cmprsk package. BUT then qgcomp can’t estimate a model created from cmprsk...

Is there a way around this?

Thank you for your help.

Version 2.2.2 (and possibly earlier versions) when displaying the results of a qgcomp fit with a gaussian family. The p value is displayed under a column "t value" and the corresponding t-statistics are not in the coefficient table. These values can be obtained directly from the qgcomp fit via the following example:

gf = qgcomp.noboot( #omitted)
# t statistic
gf$tstat
# pvalue
gf$pval

v2.2.4 fixes this issue and can be installed via:
devtools::install_github("alexpkeil1/qgcomp")

Excuse me. I am studying the association between particulate matter components (Organic matter, black carbon, sulfate, nitrate, ammonium) and death, although all five components increase the risk of death in single pollutant model, however, when I use joint exposure with qgcomp, two pollutants (sulfate, nitrate) have negative weights. I am very difficult to interpret it.

I have been attempting to adjust the yaxis using ylim as in a ggplot2 object, but I have been unsuccessful. Any suggestions?

ggplot2 instituted "visual breaking changes" that unfortunately make the weight plots look bad and miscommunicate the results. This is a reminder to shift to base R plotting functions in the future to reduce this somewhat needless dependency. Versions 2.3.5+ institute a modified plot function to deal with ggplot2 3.3.0+ (these updated versions of qgcomp also force an update of ggplot2, so user beware!)

Install version 2.2.5+:

devtools::install_github("alexpkeil1/qgcomp")

Is this easy?

Whe we use qgcomp.boot with bootstrapping method. How to estimate the weigths of joint exposure effect. Is it similar to 'qgcomp.noboot' to allow further to calculate the exposure weights?

This is just a placeholder issue because so many issues come down to "this is fixed in the newest version, which is not yet on CRAN"

Install the most recent developmental version here (this may not always be working, so use judiciously!) by running this in R (being sure you have the 'devtools' package installed):

 devtools::install_github("alexpkeil1/qgcomp")

See example below: if one of the non-exposure variables is coded as a factor, then

library(qgcomp)
N = 100
dat = data.frame(qgcomp:::.rmvnorm(N, rep(0, 3), diag(rep(1, 3))))
names(dat) <- c("y",'x','z')

a1<-qgcomp.noboot(y~., expnms='x', data = dat, family=gaussian(), q=4)
b1<-pointwisebound.noboot(a1, alpha=0.05, pointwiseref=1)
b1


# now redefine "z" to be a factor variable
dat$z = as.factor(sample(1:3, N, replace=TRUE))
a2<-qgcomp.noboot(y~., expnms='x', data = dat, family=gaussian(), q=4)
b2<-pointwisebound.noboot(a2, alpha=0.05, pointwiseref=1)
b2
# throws error:
#  Error in `[.data.frame`(x$fit$data, names(zvars)) : 
# undefined columns selected 

Hi！ I'm looking for advice on how to interpret the results of qgcomp.boot. I am using multivariate logistic regression and x1,x2,x3 are all negatively correlated with Y. In qgcomp.boot analysis (adjusted for the same variables) I noticed that the weight value of x1 is negative and the effect sizes are in a different order than in the logistic regression analysis, how should I interpret the anomalous results? Is using WQS a better option? How can I understand the weight values?

Hi,

I see that you use for instance:

Please use the following instead:

if (parplan) {
  oplan <- future::plan(strategy = future::multisession)
  on.exit(plan(oplan), add = TRUE)
}

This is the proper way to do it, cf. help("plan", package = "future"). It also avoids using transparent, which should not be used in package code and will eventually be deprecated.

When I use factor(gender) in my R model with survival::Surv() and then use plot() to plot my bootstrapped survival curve I obtained using qgcomp.cox.boot I get the error: object 'gender' not found.

However, when I convert gender to a numeric variable it runs fine.

This is more complicated when I have a multi-class variable like race. I can't convert that to an ordinal variable. How do I get around this to plot the survival curves?

Thank you!!

qc.survfit1 <- qgcomp.cox.boot(Surv(disease_time, disease_state) ~ arsenic+barium+cadmium,expnms=Xnm,
deg=3,data=metals[,c(Xnm, 'disease_time', 'disease_state')],
q=100,B=500,MCsize=500)
df<-pointwisebound.boot(qc.survfit1 , alpha=0.05, pointwiseref=1)
df<-as.data.frame(df)
names(df)
ggplot(df,aes(quantile,rr))+geom_line()+
geom_ribbon(aes(ymin=ll.rr,ymax=ul.rr),alpha=0.3)

I want to use this to plot exposure-response curve, but it does run.

Error in pointwisebound.boot(qc.survfit1, alpha = 0.05, pointwiseref = 1) :
This function does not work with this type of qgcomp fit

Please consider adding URL and BugReports fields to this package DESCRIPTION file so that it is easy to find this repo from https://cran.r-project.org/package=qgcomp.

m4<-qgcomp.noboot(death~as.factor(city)+as.factor(year)+
tmean.win+tmean.sum+ tsd.sum+tsd.win+no2+co+so2+03+offset(Logpop),
expnms=comp, data=data, q=20, family=quasipoisson())

when I add", family=quasipoisson()",one error exist:

Error in print.qgcompfitCxD : object "testtype" not found
In addition: Warning message:In print.agcompfit(x) :
The quasipoisson distribution has not been tested with agcomp! Please use with extreme caution and check results thoroughly with simulated data to ensure it works.

When breaks are specified the plot function for qgcomp.cox.boot results in an error:

Error in 0:(x$q - 1) : argument of length 0

Is this because the quartiles are set to NULL?

I apologize, but I'm a little confused following the introduction posted for qgcomp.

When estimating the hazard function using qgcomp.cox.boot you explain you can obtain the hazard ratios as a function of overall exposure:

examining the overall hazard ratio as a function of overall exposure
hrs_q = exp(matrix(c(0,0,1,1,2,4,3,9), ncol=2, byrow=TRUE)%*%qc.survfit4$msmfit$coefficients)
colnames(hrs_q) = "Hazard ratio"
print("Hazard ratios by quartiles (min-25%,25-50%, 50-75%, 75%-max)")

However, it doesn't seem that you can easily obtain the confidence intervals of these hazards ratios? I'm probably missing something.

Hi there! The broom dev team just ran reverse dependency checks on the upcoming broom 0.7.0 release and found new NOTEs for the CRAN version of this package. I've pasted the results below, which seem to result from our decision to no longer export the finish_glance() function in favor of adding more specific calls to stats::BIC() et al. (see this issue).

• checking dependencies in R code ... NOTE

  Missing or unexported object: ‘broom::finish_glance’
  

We hope to submit this new version of the package to CRAN in the coming weeks. If you encounter any problems fixing these issues, please feel free to reach out!🙂

using "predict" on a qgcomp object gives the following error for me:

>fitted_qgcomp = qgcomp.noboot(defect~., expnms=names(X), data = andata[,c(names(X), names(W), "defect")])
>predict.qgcompfit(fitted_qgcomp)
Error in UseMethod("predict") : 
  no applicable method for 'predict' applied to an object of class "data.frame"

Until this is fixed, a workaround is to explicitly call the qgcomp function that predicts:
>qgcomp:::predict.qgcompfit(fitted_qgcomp)

This bug is in v2.1.2 and was introduced when fixing other plotting bugs.

v2.1.3 fixes this issue and can be installed via:
devtools::install_github("alexpkeil1/qgcomp")

Example here:

 set.seed(50)
 n=100
 dat <- data.frame(y=rbinom(n, 1, 0.5), x1=runif(n), x2=runif(n), z=runif(n))


 ft2 = qgcomp.boot(f=y ~ x1 + x2, expnms = c('x1', 'x2'), data=dat, q=4, family=binomial(), 
              B = 200, parallel=TRUE, bayes = FALSE)
 plot(ft2)

How the plot looks in v2.1.3

How the plot looks in v2.1.2 (model bounds are incorrect)

The bounds can also just be turned off via

plot(ft2,modelband = FALSE)

There is a error when trying to plot an object from qgcomp.cox.boot if non-exposure covariates are included, when using the CRAN version of the package (1.2.0).

This issue is fixed in v1.2.6, which (until 1.3.0 is released on CRAN) can be installed via devtools::install_github("alexpkeil1/qgcomp")

This error is raised when the data given to qgcomp.boot is a tibble. Here is a minimal working example:

library(dplyr)
library(qgcomp)

N = 100
dat <- tibble(i = 1:N) %>%
  mutate(
    u = runif(N)*0.75 + .125,
    x1 = rnorm(N, u),
    x2 = rnorm(N, u),
    y = rnorm(N) + x1 - x2*x1
  )

# error related

qgcomp.boot(y ~ x1 + x2, expnms = c("x1", "x2"), data=dat)
#Error: Must use a vector in `[`, not an object of class matrix.
#Call `rlang::last_error()` to see a backtrace

Until this is fixed, a simple workaround is to convert the data to a data.frame, as in here:

qgcomp.boot(y ~ x1 + x2, expnms = c("x1", "x2"), data=data.frame(dat))

When we conduct an non-linear terms using qgcomp.boot. How to interpret the results of non-linearity ? What’s the means of the p vaule of psi1 and psi2, respectively.

Nonsense results when using the code supplied in the help files of the qgcomp package can occur when using a survival model with the "pool" function from mice.

One can still pool the results by hand using Rubin's rules, but this issue should be fixed internally in the next version of qgcomp

Just trying to understand the default qgcomp is using when there is a lonely PSU. From what I understand the options can be:

1. Assigning that participant to another PSU (done by the user)
2. Centering the lonely PSU around the sample grand mean
3. Make no contribution to the variance

I know that I have a lonely PSU, but qgcomp.cox.noboot and qgcomp.cox.boot run without a warning despite the lonely PSU. So there must be a default way of handling it. Maybe there is something I'm missing? Just want to know what the default is so I can clarify for reviewers.

Thanks!

Other distributions are currently being tested. Currently, the "poisson" distribution will not produce errors when used with ggcomp.boot and qgcomp.noboot functions. There is some further testing that needs to be done to ensure that there are no remaining issues to resolve. Use with caution and, if possible, test that results make sense using simulated data!