8
$\begingroup$

I want to recreate a regression model based on what was given in a scientific paper. They gave intercept and coefficient terms.

I know how to create regression models in R, but is this possible to do without the original database?

I would use these models on my own database to perform model comparison and test their predictive capabilities.

The special case here is that I am mostly interested in logistic regression. But I guess this question is scalable to all types of regression models.

So in other words: how can we create regression model objects (e.g. glm) using only beta values.

$\endgroup$
3
  • 1
    $\begingroup$ A working example of doing this with a logistic regression model, starting with the coefficients only, appears in the first block of code at stats.stackexchange.com/a/40609/919. An example with Poisson models is at stats.stackexchange.com/a/45789/919. Note that at a minimum you will also need somehow to specify the values of the independent variables. $\endgroup$
    – whuber
    Commented Aug 17, 2016 at 18:41
  • $\begingroup$ Related SO question that didn't get any answer : stackoverflow.com/questions/56703403/…. Building the P(Y=1) function from coefficients is one thing but embedding it in a standard R object is another. I have reached the conclusion that it is not something you want to achieve. A standard R model object will allow you to use a lot of functions but most of them will give you bad results because the glm model contain a lot more informations : Std. Error, z value, Pr(>|z|) for coefficients. $\endgroup$ Commented Aug 20, 2019 at 13:38
  • 1
    $\begingroup$ This does not appear to me to be possible. The regression model includes more information than just the coefficient estimators (e.g., leverage values, standard errors, etc.). For a linear model, the output is fully determined by the correlation matrix for all the vectors in the model plus the norms of those vectors. Thus, if you don't have access to the original data, you would still need to have this information to generate the full regression output. (Of course, you can generate parts of the regression output ---e.g., regression line--- with incomplete information.) $\endgroup$
    – Ben
    Commented Dec 23, 2020 at 23:38

2 Answers 2

5
$\begingroup$

You lack information about the values of the explanatory variables. Exploring this issue reveals some useful things to know about logistic regression. For instance, you can go far beyond simple rules of thumb concerning how much data you need for estimating the parameters and you can develop a better intuition for what can and cannot be accomplished depending on simple characteristics of the dataset (sizes, locations, and ranges).


Let me explain with some examples. These scatterplots show data generated with the coefficients $\beta = (1,1)$ (and the usual logistic link) for various sets of $x$ values. That is, the log odds of the expected response $y$ is $\beta_1 + \beta_2 x = 1 + x.$ (This is fully general, because every univariate logistic regression has this form when you choose suitable units of measurement for $x.$)

enter image description here

To appreciate these plots, understand that

  • The datasets are all the same size, with $n = 100$ observations each.

  • The $x$ values are equally spaced across their ranges.

  • The responses ($y$ values) of $0$ and $1$ have been vertically jittered to avoid overplotting.

  • The gray curves are all parts of the same model with $\beta=(1,1)$ used to generate the data. They show how the chance of a response equal to $1$ varies with $x.$

  • The posted equations are the fitted logistic models to each dataset in the form "$\hat\beta_0(\operatorname{se}(\hat\beta_0))\ + \hat\beta_1(\operatorname{se}(\hat\beta_1))x:$" that is, the standard errors of estimate appear in parentheses after each estimated parameter. (The fits were computed with the glm function in R.)

Observe how different these datasets appear. They differ in the ranges and locations of the $x$ values:

  1. The ranges across the top in examples 1 through 4 are all equal to $2$ but the locations are shifted to right as you go from left to right.

  2. The locations across the bottom in examples 5 through 8 are all centered at $x=0$ but the ranges increase from $2$ to $20.$

Consequently, most of these datasets cover small and very different portions of the full logistic sigmoidal curve, as illustrated by the gray curves. That's how they manage to have such different appearances. (But, as you can check by inspection, all coefficient estimates are within a couple of standard errors of the true model coefficients $\beta = (1,1).$)


A very close examination of the standard errors will reveal a general truth: the standard errors are smallest when most of the logistic curve is included, as in examples 7 and 8, and is largest when the curve varies only a little, as in examples 1 (the left tail), 4 (the right tail), and even 5 (the middle). This implies you can get some hints about the likely $x$ values from the p-values reported in the paper, provided you know how many $x$ values are involved (which is usually the case): caeteris paribus, tiny p-values suggest you're in a situation like examples 6 through 8 while large p-values indicate otherwise.

(This gets much more complicated with multiple regression because the p-values depend on the geometric relationships among the various explanatory variables, making it unlikely you could infer much about the $x$ values used in the paper.)

Now (returning to the simple logistic regression case with a single explanatory variable), if you can estimate the range of the $x$ values and the number of $x$ values, you can simulate sample data by generating likely values of $x$ and drawing random binomial responses according to the model. This is illustrated in the first block of R code at https://stats.stackexchange.com/a/40609/919, but to be fully explicit I offer this general-purpose implementation (used to create the plots).

#
# Generate data according to a specified logistic regression model.
# `x` are the explanatory values (which can be a matrix for multiple regression).
# `beta` are the coefficients.
# When `intercept` is TRUE, beta[1] is taken to be the intercept term.
#
# Returns a matrix of `n` columns, one for each independent data set.
#
rlreg <- function(n = 1, beta, x, intercept = TRUE) {
  if(isTRUE(intercept)) x <- cbind(1, x)
  if (!is.matrix(x)) x <- matrix(x, ncol = length(beta))

  y.hat <- x %*% beta        # The inverse link values
  h <- 1 / (1 + exp(-y.hat)) # The expectations
  matrix(rbinom(n * length(h), 1, h), ncol = n)
}

As an example of the use of rlreg, here are three independent simulated datasets sharing a common set of $x$ values (leftmost column).

set.seed(17)
df <- data.frame(x = seq(-4, 4, length.out = 9))
df$y <- rlreg(3, c(1, 1), df$x)
df
   x y.1 y.2 y.3
1 -4   0   0   0
2 -3   1   0   0
3 -2   0   0   1
4 -1   1   1   1
5  0   1   0   0
6  1   1   0   1
7  2   1   1   1
8  3   1   1   1
9  4   1   1   1
$\endgroup$
5
  • 1
    $\begingroup$ Great and detailed answer. Thank you. $\endgroup$
    – WojciechF
    Commented Mar 24 at 10:59
  • $\begingroup$ thanks whuber: would you pls also share the r script graphic output part? $\endgroup$
    – Maximilian
    Commented Mar 29 at 11:38
  • 1
    $\begingroup$ @Maximilian It's really too long to post here, the methods are routine, and it would be distracting anyway. But--briefly--it exploits geom_jitter in the ggplot2 package to plot the points at random heights, geom_text to post the formulas, and facet_wrap to create the array of panels. The curves themselves are plotted using a second data frame for which the conditional expectations were computed. $\endgroup$
    – whuber
    Commented Apr 3 at 20:30
  • $\begingroup$ the conditional expectation is the troubling part I cannot seem to replicate. thanks. $\endgroup$
    – Maximilian
    Commented Apr 4 at 15:35
  • $\begingroup$ @Maximilian I used the logistic link function explicitly, but one could always just resort to predict.glm. $\endgroup$
    – whuber
    Commented Apr 4 at 15:51
0
$\begingroup$

I think all you need to do is "score" (create a new column in your database that contains the predicted values for each record in your database) using the regression model coefficients and functional form of the model (for linear regression example, y = XB where y is the predicted value from the regression model, X is your database, and B is a vector with the model coefficients).

I'm not sure of the exact functional form of your regression model but in R you can write the equation from the command line:

y <- a + b*x

Hope this helps

$\endgroup$
4
  • $\begingroup$ statistics.ats.ucla.edu/stat/r/dae/logit.htm $\endgroup$ Commented Aug 17, 2016 at 17:09
  • $\begingroup$ Link to R code in previous comment $\endgroup$ Commented Aug 17, 2016 at 17:11
  • 1
    $\begingroup$ The question focuses on GLMs, for which this approach will not work. For linear regression your model is incomplete: you need to add random noise to y. Its variance is one more parameter of the model. $\endgroup$
    – whuber
    Commented Aug 17, 2016 at 17:30
  • 1
    $\begingroup$ Agree if goal is simulating draws from model distribution. Prior post predicts model mean for a given set of covariates. $\endgroup$ Commented Aug 18, 2016 at 19:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.