0

I have been given an Rdata file containing a large number of inputs and outputs from a regression model. I have been able to extract the data analyzed by the model and reproduce the parameter estimates. However, when I attempt to use the original predict statement I receive an error even though the predict statement does not return an error when applied to the model stored in the Rdata file.

I am hoping there is enough information presented below that someone may be able to tell me how to correct my predict statement, my.probs, even though I am not providing a functional reproducible example. This, I think, is the first time I have ever posted a question here without providing such an example. The data set contains > 100,000 observations, is somewhat sensitive and I am unsure how I would reproduce the Rdata file.

library(msm)
library(MASS)
library(pscl)

# model output returned when extracting the model name from the `Rdata` file
original.model
# Call:
# zeroinfl(formula = AA ~ log(BB) + CC + DD + CC:DD | log(BB) + DD, data = original.data, 
#     offset = log(EE), dist = "negbin")
# 
# Count model coefficients (negbin with log link):
#           (Intercept)  log(BB)      CC3      CC4      CC5 DDPrivate CC3:DDPrivate CC4:DDPrivate CC5:DDPrivate  
#              -2.05317  0.31178 -0.41402 -0.71208 -0.92290   0.17878      -0.18476      -0.18674       0.07307  
# Theta = 0.8551 
# 
# Zero-inflation model coefficients (binomial with logit link):
#           (Intercept)  log(BB)        DDPrivate  
#                1.6724 -0.5022            0.9742  
#  
# Warning message:
# In deparse(x$call, width.cutoff = floor(getOption("width") * 0.85)) :
#   invalid 'cutoff' value for 'deparse', using default

# data for new observation for use in the predict statement
new.data
#         DD      EE   CC               BB
# 1  Private       1    4         1118.948

str(new.data)
#'data.frame': 1 obs. of  4 variables:
# $ DD       : Factor w/ 2 levels "Public","Private": 2
# $ EE       : num 1
# $ CC       : Factor w/ 4 levels "2","3","4","5": 3
# $ BB: num 1119

original.probs <- predict(original.model, new.data, type='prob')
original.probs
# truncated probabilities returned by the predict statement.  These sum to one if vector not truncated
c(0.7534319,    0.1552296,    0.05681916,   0.02133936,   0.008116065,  0.003110019,  0.001197667)

# reproduce the original model
my.version <- zeroinfl(formula = AA ~ log(BB) + CC + DD + CC:DD | log(BB) + DD, offset = log(EE), dist = "negbin")

# Error returned by the predict statement
my.probs <- predict(my.version, new.data, type='prob')
my.probs
# Error in exp(X %*% object$coefficients$count + offsetx)[, 1] : 
#   incorrect number of dimensions
# In addition: Warning message:
# In X %*% object$coefficients$count + offsetx :
#   Recycling array of length 1 in array-vector arithmetic is deprecated.
#   Use c() or as.vector() instead.
Mark Miller
  • 12,483
  • 23
  • 78
  • 132
  • It's easier to help you if you include a simple [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) with sample input and desired output that can be used to test and verify possible solutions. – MrFlick Apr 27 '20 at 19:34
  • 1
    Yes, I know that. I have been here 10 years. I just do not know how to reproduce the Rdata file that would be the first step in doing what you suggest. – Mark Miller Apr 27 '20 at 19:35
  • Well, it's very unclear how we can help you then. Maybe try turning on debugging for the `preoduct()` function and step though it in each of the cases to see when something different happens. – MrFlick Apr 27 '20 at 19:40

1 Answers1

0

The predict function worked after I grouped the input variables into a data.frame and included the data option in the zeroinfl model statement:

my.data <- data.frame(AA = AA,
                      BB = BB,
                      CC = CC,
                      DD = DD,
                      EE = EE)

my.version <- zeroinfl(formula = AA ~ log(BB) + CC + DD + CC:DD | log(BB) + DD, 
                               offset = log(EE), dist = "negbin", data = my.data)

summary(my.version)

my.probs <- predict(my.version, new.data, type='prob')
my.probs
Mark Miller
  • 12,483
  • 23
  • 78
  • 132