print and coef Methods for survreg Differ

Hello,



I'm working on a survreg model where the full data are subset for 
modeling individual parts of the data separately.  When subsetting, the 
fit variable ("treatment" in the example below) has levels that are not 
in the data.
A
work-around for this is to drop the levels, but it seems inaccurate to 
have the `coef()` method provide zero as the coefficient for the level 
without data.



Why does coef(model) provide zero as the coefficient for treatment 
instead of NA?  Is this a bug?



Thanks,



Bill



``` r

library(survival)

library(emmeans)



my_data <-

 data.frame(

   value=c(rep(1, 5), 6:10),

   treatment=factor(rep(c("A", "B"), each=5), levels=c("A", "B", "C"))

 )

my_data$cens <- c(0, 1)[(my_data$value == 1) + 1]



model <- survreg(Surv(time=value, event=cens)~treatment, data=my_data)

#> Warning in survreg.fit(X, Y, weights, offset, init = init, 
controlvals =

#> control, : Ran out of iterations and did not converge

coef(model)

#> (Intercept)  treatmentB  treatmentC

#>  0.08588218  2.40341893  0.00000000

model$coef

#> (Intercept)  treatmentB  treatmentC

#>  0.08588218  2.40341893          NA

model$coefficients

#> (Intercept)  treatmentB  treatmentC

#>  0.08588218  2.40341893  0.00000000

print(model)

#> Call:

#> survreg(formula = Surv(time = value, event = cens) ~ treatment,

#>     data = my_data)

#>

#> Coefficients: (1 not defined because of singularities)

#> (Intercept)  treatmentB  treatmentC

#>  0.08588218  2.40341893          NA 

#>

#> Scale= 0.09832254

#>

#> Loglik(model)= 4.9   Loglik(intercept only)= -15

#>  Chisq= 39.92 on 2 degrees of freedom, p= 2.15e-09

#> n= 10

summary(model)

#>

#> Call:

#> survreg(formula = Surv(time = value, event = cens) ~ treatment,

#>     data = my_data)

#>               Value Std. Error     z      p

#> (Intercept)  0.0859     0.0681  1.26   0.21

#> treatmentB   2.4034     0.2198 10.93 <2e-16

#> treatmentC   0.0000     0.0000    NA     NA

#> Log(scale)  -2.3195     0.0000  -Inf <2e-16

#>

#> Scale= 0.0983

#>

#> Weibull distribution

#> Loglik(model)= 4.9   Loglik(intercept only)= -15

#>  Chisq= 39.92 on 2 degrees of freedom, p= 2.1e-09

#> Number of Newton-Raphson Iterations: 30

#> n= 10

ref_grid(model)

#> Error in ref_grid(model): Something went wrong:

#>  Non-conformable elements in reference grid.



my_data_correct_levels <- my_data

my_data_correct_levels$treatment <-
droplevels(my_data_correct_levels$treatment)



model_correct <- survreg(Surv(time=value, event=cens)~treatment,
data=my_data_correct_levels)

#> Warning in survreg.fit(X, Y, weights, offset, init = init, 
controlvals =

#> control, : Ran out of iterations and did not converge

coef(model_correct)

#> (Intercept)  treatmentB

#>  0.08588218  2.40341893

print(model_correct)

#> Call:

#> survreg(formula = Surv(time = value, event = cens) ~ treatment,

#>     data = my_data_correct_levels)

#>

#> Coefficients:

#> (Intercept)  treatmentB

#>  0.08588218  2.40341893

#>

#> Scale= 0.09832254

#>

#> Loglik(model)= 4.9   Loglik(intercept only)= -15

#>  Chisq= 39.92 on 1 degrees of freedom, p= 2.65e-10

#> n= 10

summary(model_correct)

#>

#> Call:

#> survreg(formula = Surv(time = value, event = cens) ~ treatment,

#>     data = my_data_correct_levels)

#>               Value Std. Error     z      p

#> (Intercept)  0.0859     0.0681  1.26   0.21

#> treatmentB   2.4034     0.2198 10.93 <2e-16

#> Log(scale)  -2.3195     0.0000  -Inf <2e-16

#>

#> Scale= 0.0983

#>

#> Weibull distribution

#> Loglik(model)= 4.9   Loglik(intercept only)= -15

#>  Chisq= 39.92 on 1 degrees of freedom, p= 2.6e-10

#> Number of Newton-Raphson Iterations: 30

#> n= 10

ref_grid(model_correct)

#> 'emmGrid' object with variables:

#>     treatment = A, B

#> Transformation: "log"

```



<sup>Created on 2021-02-23 by the [reprex
package](https://reprex.tidyverse.org) (v1.0.0)</sup>


[[alternative HTML version deleted]]