georgheinze / coxphf

2 stars 0 forks source link

Unable to reconstruct coxphftest example #13

Closed karl-an closed 1 year ago

karl-an commented 1 year ago

Unable to reconstruct coxphftest example as I get the following error:

Error in terms.formula(formula, data = data) : invalid term in model formula

I am a little lost here. I am trying to extract the model test results for the summary table: Likelihood ratio test=0.9842619 on 3 df, p=0.80506, n=10 Wald test = 0.758843 on 3 df, p = 0.8592837

> library(survival)
> testdata <- data.frame(list(start=c(1, 2, 5, 2, 1, 7, 3, 4, 8, 8),
+                             stop =c(2, 3, 6, 7, 8, 9, 9, 9,14,17),
+                             event=c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0),
+                             x1 =c(1, 0, 0, 1, 0, 1, 1, 1, 0, 0),
+                             x2 =c(0, 1, 1, 1, 0, 0, 1, 0, 1, 0),
+                             x3 =c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0)))
> summary( coxphf( formula=Surv(start, stop, event) ~ x1+x2+x3, data=testdata))
coxphf(formula = Surv(start, stop, event) ~ x1 + x2 + x3, data = testdata)

Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood 

         coef se(coef) exp(coef) lower 0.95 upper 0.95     Chisq         p
x1 -0.3119630 0.901733 0.7320086  0.1331772   3.898528 0.1437621 0.7045693
x2  0.8176286 1.007481 2.2651220  0.4284707  18.866707 0.8833246 0.3472927
x3 -0.1496022 1.026306 0.8610504  0.1150115   5.119178 0.0257560 0.8724976

Likelihood ratio test=0.9842619 on 3 df, p=0.80506, n=10
Wald test = 0.758843 on 3 df, p = 0.8592837

Covariance-Matrix:
           x1         x2         x3
x1  0.8131224 -0.2902941  0.3959233
x2 -0.2902941  1.0150183 -0.4745365
x3  0.3959233 -0.4745365  1.0533046
> # testing H0: x1=0, x2=0
> coxphftest( formula=Surv(start, stop, event) ~ x1+x2+x3, test=~x1+x2, data=testdata)
Error in terms.formula(formula, data = data) : 
  invalid term in model formula
> fitt<- coxphf( formula=Surv(start, stop, event) ~ x1+x2+x3*stop, data=testdata, pl=FALSE)
> test <- coxphf(formula=Surv(start, stop, event) ~ x1+x2+x3*stop, data=testdata, adapt=c(1, 1, 0, 0))
> # PLR p-value for x3 + x3:stop
> pchisq((fitt$loglik[2]-test$loglik[2])*2, 2, lower.tail=FALSE)
[1] 0.9970177
gregorsteiner commented 1 year ago

Hi, I could not reproduce this. Your code runs without errors on my machine, see below:

library(coxphf)
library(survival)
testdata <- data.frame(list(start=c(1, 2, 5, 2, 1, 7, 3, 4, 8, 8),
                                           stop =c(2, 3, 6, 7, 8, 9, 9, 9,14,17),
                                           event=c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0),
                                           x1 =c(1, 0, 0, 1, 0, 1, 1, 1, 0, 0),
                                           x2 =c(0, 1, 1, 1, 0, 0, 1, 0, 1, 0),
                                           x3 =c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0)))
summary( coxphf( formula=Surv(start, stop, event) ~ x1+x2+x3, data=testdata))
#> coxphf(formula = Surv(start, stop, event) ~ x1 + x2 + x3, data = testdata)
#> 
#> Model fitted by Penalized ML
#> Confidence intervals and p-values by Profile Likelihood 
#> 
#>          coef se(coef) exp(coef) lower 0.95 upper 0.95     Chisq         p
#> x1 -0.3119630 0.901733 0.7320086  0.1331772   3.898528 0.1437621 0.7045693
#> x2  0.8176286 1.007481 2.2651220  0.4284707  18.866707 0.8833246 0.3472927
#> x3 -0.1496022 1.026306 0.8610504  0.1150115   5.119178 0.0257560 0.8724976
#> 
#> Likelihood ratio test=0.9842619 on 3 df, p=0.80506, n=10
#> Wald test = 0.758843 on 3 df, p = 0.8592837
#> 
#> Covariance-Matrix:
#>            x1         x2         x3
#> x1  0.8131224 -0.2902941  0.3959233
#> x2 -0.2902941  1.0150183 -0.4745365
#> x3  0.3959233 -0.4745365  1.0533046
coxphftest( formula=Surv(start, stop, event) ~ x1+x2+x3, test=~x1+x2, data=testdata)
#> coxphftest(formula = Surv(start, stop, event) ~ x1 + x2 + x3, 
#>     data = testdata, test = ~x1 + x2)
#> Model fitted by Penalized ML 
#> 
#> Factors fixed as follows:
#> x1 x2 x3 
#>  0  0 NA 
#> 
#> Likelihoods:
#> Restricted model       Full model       difference 
#>       -7.0492717       -6.6053140        0.4439576 
#> 
#> Likelihood ratio test=0.8879153 on 2 df, p=0.6414926

sessionInfo()
#> R version 4.2.2 (2022-10-31 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 17763)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=English_Austria.1252  LC_CTYPE=English_Austria.1252   
#> [3] LC_MONETARY=English_Austria.1252 LC_NUMERIC=C                    
#> [5] LC_TIME=English_Austria.1252    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] survival_3.4-0 coxphf_1.13.3 
#> 
#> loaded via a namespace (and not attached):
#>  [1] rstudioapi_0.14   knitr_1.42        magrittr_2.0.3    splines_4.2.2    
#>  [5] lattice_0.20-45   R.cache_0.16.0    rlang_1.0.6       fastmap_1.1.0    
#>  [9] styler_1.9.0      tools_4.2.2       grid_4.2.2        xfun_0.37        
#> [13] R.oo_1.25.0       cli_3.4.1         withr_2.5.0       htmltools_0.5.4  
#> [17] yaml_2.3.5        digest_0.6.29     lifecycle_1.0.3   Matrix_1.5-1     
#> [21] purrr_1.0.1       vctrs_0.5.2       R.utils_2.12.2    fs_1.5.2         
#> [25] glue_1.6.2        evaluate_0.20     rmarkdown_2.20    reprex_2.0.2     
#> [29] compiler_4.2.2    generics_0.1.3    R.methodsS3_1.8.2

Are you sure you are using version 1.13.3? We recently fixed an issue with terms(), so this might be the problem.

karl-an commented 1 year ago

There seems to be a problem with one of my packages. But its hard to know which one it is from the error message.


> library(coxphf)
> library(survival)
> testdata <- data.frame(list(start=c(1, 2, 5, 2, 1, 7, 3, 4, 8, 8),
+                             stop =c(2, 3, 6, 7, 8, 9, 9, 9,14,17),
+                             event=c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0),
+                             x1 =c(1, 0, 0, 1, 0, 1, 1, 1, 0, 0),
+                             x2 =c(0, 1, 1, 1, 0, 0, 1, 0, 1, 0),
+                             x3 =c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0)))
> summary( coxphf( formula=Surv(start, stop, event) ~ x1+x2+x3, data=testdata))
coxphf(formula = Surv(start, stop, event) ~ x1 + x2 + x3, data = testdata)

Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood 

         coef se(coef) exp(coef) lower 0.95 upper 0.95     Chisq         p
x1 -0.3119630 0.901733 0.7320086  0.1331772   3.898528 0.1437621 0.7045693
x2  0.8176286 1.007481 2.2651220  0.4284707  18.866707 0.8833246 0.3472927
x3 -0.1496022 1.026306 0.8610504  0.1150115   5.119178 0.0257560 0.8724976

Likelihood ratio test=0.9842619 on 3 df, p=0.80506, n=10
Wald test = 0.758843 on 3 df, p = 0.8592837

Covariance-Matrix:
           x1         x2         x3
x1  0.8131224 -0.2902941  0.3959233
x2 -0.2902941  1.0150183 -0.4745365
x3  0.3959233 -0.4745365  1.0533046
> 
> coxphftest( formula=Surv(start, stop, event) ~ x1+x2+x3, test=~x1+x2, data=testdata)
Error in terms.formula(formula, data = data) : 
  invalid term in model formula
> sessionInfo()
R version 4.3.0 (2023-04-21 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22621)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8    LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                           LC_TIME=English_United States.utf8    

time zone: Asia/Taipei
tzcode source: internal

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] modelsummary_1.4.1         multiColl_2.0              logistf_1.25.0             interplot_0.2.3           
 [5] arm_1.13-1                 lme4_1.1-33                Matrix_1.5-4               abind_1.4-5               
 [9] gtsummary_1.7.2            igraph_1.4.2               fuzzyjoin_0.1.6            stringdist_0.9.10         
[13] prodlim_2023.03.31         PerformanceAnalytics_2.0.4 xts_0.13.1                 zoo_1.8-12                
[17] ggfortify_0.4.16           ranger_0.15.1              lsr_0.5.2                  readxl_1.4.2              
[21] oglmx_3.0.0.0              maxLik_1.5-2               miscTools_0.6-28           sure_0.2.0                
[25] reshape2_1.4.4             Hmisc_5.1-0                foreign_0.8-84             lubridate_1.9.2           
[29] forcats_1.0.0              stringr_1.5.0              dplyr_1.1.2                purrr_1.0.1               
[33] readr_2.1.4                tidyr_1.3.0                tibble_3.2.1               ggplot2_3.4.2             
[37] tidyverse_2.0.0            haven_2.5.2                MASS_7.3-60                survival_3.5-5            
[41] coxphf_1.13.3             

loaded via a namespace (and not attached):
  [1] splines_4.3.0           later_1.3.1             BiasedUrn_2.0.9         cellranger_1.1.0        rpart_4.1.19           
  [6] lifecycle_1.0.3         formula.tools_1.7.1     sf_1.0-12               globals_0.16.2          lattice_0.21-8         
 [11] vroom_1.6.3             insight_0.19.1          gofcat_0.1.2            backports_1.4.1         magrittr_2.0.3         
 [16] sass_0.4.6              rmarkdown_2.21          jquerylib_0.1.4         httpuv_1.6.9            zip_2.3.0              
 [21] askpass_1.1             gld_2.6.6               DBI_1.1.3               minqa_1.2.5             multcomp_1.4-23        
 [26] rvest_1.0.3             expm_0.999-7            quadprog_1.5-8          TH.data_1.1-2           nnet_7.3-19            
 [31] sandwich_3.0-2          gdtools_0.3.3           labelled_2.11.0         lava_1.7.2.1            listenv_0.9.0          
 [36] crul_1.3                units_0.8-2             parallelly_1.35.0       svglite_2.1.1           codetools_0.2-19       
 [41] DT_0.27                 xml2_1.3.4              tidyselect_1.2.0        httpcode_0.3.0          base64enc_0.1-3        
 [46] webshot_0.5.4           broom.helpers_1.13.0    jsonlite_1.8.5          e1071_1.7-13            ellipsis_0.3.2         
 [51] Formula_1.2-5           emmeans_1.8.5           systemfonts_1.0.4       tools_4.3.0             ragg_1.2.5             
 [56] DescTools_0.99.48       Rcpp_1.0.10             glue_1.6.2              gridExtra_2.3           mgcv_1.8-42            
 [61] xfun_0.39               withr_2.5.0             fastmap_1.1.1           boot_1.3-28.1           fansi_1.0.4            
 [66] openssl_2.0.6           digest_0.6.31           estimability_1.4.1      timechange_0.2.0        R6_2.5.1               
 [71] mime_0.12               mice_3.15.0             textshaping_0.3.6       colorspace_2.1-0        utf8_1.2.3             
 [76] generics_0.1.3          fontLiberation_0.1.0    data.table_1.14.8       class_7.3-22            httr_1.4.6             
 [81] htmlwidgets_1.6.2       pkgconfig_2.0.3         gtable_0.3.3            epiR_2.0.60             Exact_3.2              
 [86] htmltools_0.5.5         fontBitstreamVera_0.1.1 carData_3.0-5           kableExtra_1.3.4        scales_1.2.1           
 [91] lmom_2.9                knitr_1.42              rstudioapi_0.14         tzdb_0.3.0              uuid_1.1-0             
 [96] coda_0.19-4             checkmate_2.2.0         nlme_3.1-162            curl_5.0.0              nloptr_2.0.3           
[101] cachem_1.0.8            proxy_0.4-27            flextable_0.9.1         operator.tools_1.6.3    rootSolve_1.8.2.3      
[106] KernSmooth_2.23-21      parallel_4.3.0          pillar_1.9.0            grid_4.3.0              reshape_0.8.9          
[111] vctrs_0.6.2             pscl_1.5.5.1            VGAM_1.1-8              promises_1.2.0.1        car_3.1-2              
[116] xtable_1.8-4            cluster_2.1.4           htmlTable_2.4.1         evaluate_0.21           mvtnorm_1.1-3          
[121] cli_3.6.1               compiler_4.3.0          rlang_1.1.1             crayon_1.5.2            future.apply_1.10.0    
[126] classInt_0.4-9          plyr_1.8.8              stringi_1.7.12          pander_0.6.5            viridisLite_0.4.2      
[131] tables_0.9.17           munsell_0.5.0           fontquiver_0.2.1        interactionTest_1.2     hms_1.1.3              
[136] bit64_4.0.5             future_1.32.0           gfonts_0.2.0            shiny_1.7.4             gt_0.9.0               
[141] broom_1.0.4             bslib_0.4.2             bit_4.0.5               officer_0.6.2 

Thanks for the help!

gregorsteiner commented 1 year ago

Hi, I think it might be formula.tools, I have had similar problems in another package. maybe try uninstalling formula.tools