search

harrelfe / rmsb

Regression Modeling Strategies Bayesian
16 stars 2 forks source link

Problem with several rmsb errors #2

Closed albertocarm closed 3 months ago

albertocarm commented 4 years ago

I found several possible bugs while using the functions of the new rmsb package. I report the first one, summary(fit) or plot(summary(fit) does not seem to work in this example. Please check and see what you can do. Thank you.


library(rms)
#> Loading required package: Hmisc
#> Loading required package: lattice
#> Loading required package: survival
#> Loading required package: Formula
#> Loading required package: ggplot2
#> 
#> Attaching package: 'Hmisc'
#> The following objects are masked from 'package:base':
#> 
#>     format.pval, units
#> Loading required package: SparseM
#> 
#> Attaching package: 'SparseM'
#> The following object is masked from 'package:base':
#> 
#>     backsolve
library(Hmisc)
library(rmsb)
#> Warning: package 'rmsb' was built under R version 4.1.0

DAT<-readRDS("C:/Users/Alberto/Desktop/caterina bayes/datos", refhook = NULL)

DAT$RR_PAC5<-DAT$RR_PAC4
DAT$RR_PAC5[DAT$RR_PAC5==1]<-0
DAT$RR_PAC5[DAT$RR_PAC5==2]<-1
DAT$RR_PAC5[DAT$RR_PAC5==3]<-2
DAT$RR_PAC5<-factor(DAT$RR_PAC5)

dd<-datadist(DAT)
options(datadist='dd')
bsix2 <- blrm(HE6  ~ cir * rcs(Age,3)+ linf+ RR_PAC5 + Esquema2+Estadio3  + pol(EORTC), 
             ~  pol(EORTC), cppo=function(y) y, data=DAT, file='C:/Users/Alberto/Desktop/caterina bayes/mod_finalint.rds') 
bsix2
#> Frequencies of Missing Values Due to Each Variable
#>      HE6      cir      Age     linf  RR_PAC5 Esquema2 Estadio3    EORTC 
#>        0        0        0        0        0        0        0        2 
#> 
#> Constrained Partial Proportional Odds Ordinal Logistic Model
#>  
#>  blrm(formula = HE6 ~ cir * rcs(Age, 3) + linf + RR_PAC5 + Esquema2 + 
#>      Estadio3 + pol(EORTC), ppo = ~pol(EORTC), cppo = function(y) y, 
#>      data = DAT, file = "C:/Users/Alberto/Desktop/caterina bayes/mod_finalint.rds")
#>  
#>  
#>  Frequencies of Responses
#>  
#>     0  8.3 16.7   25 33.3 41.7   50 58.3 66.7   75 83.3 91.7  100 
#>     5    1    4    1    9    2   43   18   35   20   48    8   23 
#>  
#>  
#>                    Mixed Calibration/             Discrimination               Rank Discrim.    
#>                Discrimination Indexes                    Indexes                     Indexes    
#>     Obs217    LOO log L-458.6+/-12.05    g   1.15 [0.896, 1.414]    C   0.651 [0.629, 0.671]    
#>  Draws4000     LOO IC    917.2+/-24.1    gp 0.228 [0.188, 0.269]    Dxy 0.301 [0.258, 0.342]    
#>    Chains4    Effective p27.62+/-2.08    EV  0.177 [0.119, 0.23]                                
#>     p   13     B 0.214 [0.203, 0.225]    v   1.18 [0.634, 1.712]                                
#>                                          vp 0.042 [0.029, 0.055]                                
#>  
#>                 Mode Beta Mean Beta Median Beta S.E.   Lower   Upper  
#>  cir=1          -0.3048   -0.2117   -0.2311     2.5772 -5.4525  4.6297
#>  Age             0.0781    0.0814    0.0810     0.0422 -0.0042  0.1628
#>  Age'           -0.1145   -0.1195   -0.1186     0.0529 -0.2267 -0.0194
#>  linf            0.2816    0.2868    0.2916     0.3336 -0.3699  0.9325
#>  RR_PAC5=1      -0.5909   -0.5951   -0.5896     0.3487 -1.2717  0.0925
#>  RR_PAC5=2      -0.7137   -0.7289   -0.7364     0.3644 -1.4355 -0.0187
#>  Esquema2       -0.0156   -0.0222   -0.0305     0.3128 -0.6002  0.6123
#>  Estadio3=1      0.0299    0.0286    0.0284     0.3145 -0.5739  0.6349
#>  Estadio3=2     -0.3872   -0.4075   -0.4141     0.5537 -1.5124  0.6412
#>  EORTC           0.2385    0.2528    0.2524     0.0980  0.0573  0.4420
#>  EORTC^2        -0.0013   -0.0014   -0.0014     0.0007 -0.0028 -0.0001
#>  cir=1 * Age    -0.0077   -0.0100   -0.0097     0.0550 -0.1200  0.0943
#>  cir=1 * Age'    0.0103    0.0139    0.0117     0.0723 -0.1251  0.1555
#>  EORTC x f(y)    0.0686    0.0800    0.0776     0.0558 -0.0266  0.1880
#>  EORTC^2 x f(y) -0.0003   -0.0004   -0.0004     0.0004 -0.0012  0.0003
#>                 Pr(Beta>0) Symmetry
#>  cir=1          0.4655     1.02    
#>  Age            0.9755     1.08    
#>  Age'           0.0092     0.96    
#>  linf           0.8090     0.98    
#>  RR_PAC5=1      0.0460     1.01    
#>  RR_PAC5=2      0.0217     1.00    
#>  Esquema2       0.4622     0.99    
#>  Estadio3=1     0.5370     1.04    
#>  Estadio3=2     0.2260     1.01    
#>  EORTC          0.9955     1.00    
#>  EORTC^2        0.0190     1.02    
#>  cir=1 * Age    0.4298     0.97    
#>  cir=1 * Age'   0.5700     1.03    
#>  EORTC x f(y)   0.9290     1.12    
#>  EORTC^2 x f(y) 0.1668     0.88    

summary(bsix2) # not working
#> Error in xd %*% beta: argumentos no compatibles

plot(summary(bsix2)) # not working
#> Error in xd %*% beta: argumentos no compatibles

´´´
harrelfe commented 4 years ago

Since you do not wish to follow guidelines for posting minimal working examples I don't know how to help.

albertocarm commented 4 years ago

I'm so sorry I did it wrong. I have reduced the code to the essentials, added a simulated database, and after the model I include the functions that I can't make work. Thank you very much for your time and valuable support, Professor Harrell.

library(rms)
library(Hmisc)
library(rmsb)
data<- data.frame(HE6=as.integer(cut2(sample(1:100, 200, replace = TRUE),g=10)), Age = sample(1:85, 200, replace = TRUE), EORTC = sample(1:100, 200, replace = TRUE), cir=rbinom(200, 1,.5))
head(data)

dd<-datadist(data)
options(datadist='dd')

bsx <- blrm(HE6  ~ cir*rcs(Age,3)+ pol(EORTC), 
            ~ cir*rcs(Age,3)+ pol(EORTC), cppo=function(y) y, data=data) 
summary (bsx)
plot(summary(bsx)) 

newdata <- data.frame(cir=c(0,1), Age=53, EORTC= 75)
predict(bsx, newdata, type='fitted') # 

ggplot (Predict (bsx), abbrev =TRUE , ylab=NULL) # ci out of range

exprob <- ExProb (bsx)

fun5<-function(x) exprob (x, y=5)
fun6<-function(x) exprob (x, y=6)
fun9<-function(x) exprob (x, y=9)

M <- Mean (bsx)
qu <-Quantile (bsx)
med <- function (lp) qu(.5 , lp)

nom.ord<- nomogram(bsx, fun=list(Mean=M, "Prob y>=5"=fun5,"Prob y >=6"=fun6, "Prob y>=9"=fun9), lp=F)

plot(nom.ord)
albertocarm commented 4 years ago

I downloaded the new version of rms/rmsb. The plot for partial effects now works fine. The function 'summary' works for the previous example, but not if more variables are added. For example:

data<- data.frame(HE6=as.integer(cut2(sample(1:100, 200, replace = TRUE),g=10)), 
                  Age = sample(1:85, 200, replace = TRUE), EORTC = sample(1:100, 200, replace = TRUE), 
                  linf=rbinom(200, 1,.5),
                  cir=rbinom(200, 1,.5))
bsx <- blrm(HE6  ~ cir*rcs(Age,3)+ pol(EORTC)+linf, 
            ~ cir*rcs(Age,3)+ pol(EORTC), cppo=function(y) y, data=data) 
summary (bsx)
plot(summary(bsx)) 
plot(summary (bsx, ycut=5))

The nomogram doesn't work either.

exprob <- ExProb (bsx)

fun5<-function(x) exprob (x, y=5)
fun6<-function(x) exprob (x, y=6)
fun9<-function(x) exprob (x, y=9)

M <- Mean (bsx)
qu <-Quantile (bsx)
med <- function (lp) qu(.5 , lp)

nom.ord<- nomogram(bsx, fun=list(Mean=M, "Prob y>=5"=fun5,"Prob y >=6"=fun6, "Prob y>=9"=fun9), lp=F)

plot(nom.ord)