# R: How to create equation from logistic regression glm() model coefficient estimates? [closed]

I am trying to create a formula from the published glm() model coefficients estimates, i.e. I don't have raw data to reproduce the model. Further, I want to apply this formula for my dataframe of predictors.

Formula predicts forest damage based on tree age, species, height, etc. . I know that family is binomial with link-logit function. The model contains quantitative and categorical variables, and interaction between two specific variables, in total 17!! of them. This post shows how to convert estimates into formula, but contains just two variables. How to make it more reusable if I have more of them?

Let's say these are coefficients that I need to reproduce (dummy example, generated by hand):

                           Estimate Std. Error z value Pr(>|z|)
(Intercept)                 -28.747   20.31805  -1.224    0.071
speciesspruce                14.916   14.99081   0.878    0.320
speciesother                  5.642    8.43308   1.391    0.144
log(height)                   2.940    5.13024   1.259    0.208
age6-10                       8.639    4.24270   1.645    0.001
age>10                        0.084    1.24922   0.067    0.096
speciesspruce:log(height)    -5.074    3.57133  -0.860    0.090
speciesother:log(height)   -6.05681    3.00345  -1.514    0.030


How can I efficiently, and correctly get working formula out of estimates, and apply it to predict my data?

I am very beginner in gml() models, and the math behind seems pretty challenging. For example, simple linear regression formula seems pretty straightforward: y = alpha + beta*x.

But, in my glm(), it seems more complicated:

• y is log transformed (link=logit),
• reference class from categorical data is as pre-defined missing from output table (e.g. factor 'species' has three levels: 'pine, spruce, other' but in glm() estimates only 'speciesspruce', 'speciesother' are reported, but 'speciespine' is missing?)
• how to correctly specify interaction between two factors, species and height?

Here is my working example with same variables names and same interactions as desired output. Idea behind was to later replace manually my dummy coefficents by the real ones, get from desired one. I guess this approach is not the most straighforward one and maybe you have some more efficient and proper solution?

EDIT:

Can I simply replace my dummy estimate by the real ones, as suggested by @dave2e? how to get estimates for missing categorical variables?

mod$coefficients<-c(-28.747, 14.916, 5.642, 2.940, 8.639, 0.084, -5.074, -6.05681)  EDIT 2: I would like to use the REAL coefficient to predict new values from input data.frame table using predict.glm(updatedModel, df). How can I accomplish this? I have updated my code as well My dummy example: set.seed(42) row.num = 10 species <- factor(rep(c("pine", "spruce", "other"), each = row.num), levels = c("pine", "spruce", "other")) height <- c(runif(row.num, min = 10, max = 200), runif(row.num, min = 0, max = 100), runif(row.num, min = 30, max = 150)) age <- factor(sample(c("0-5", "6-10", ">10"), length(species), replace = TRUE), levels = c("0-5", "6-10", ">10")) damage <- rbinom(length(species), 1, 0.4) # put data together df<-data.frame(species, height, age, damage) # stringsAsFactors = FALSE # Create formula and specify interactions: mod<- glm(formula = damage ~ species + log(height) + age + log(height)*species, data = df, family = binomial(link = "logit")) summary(mod) Call: glm(formula = damage ~ species + log(height) + age + log(height) * species, family = binomial(link = "logit"), data = df) Deviance Residuals: Min 1Q Median 3Q Max -1.61759 -0.93763 -0.01525 0.83786 2.18666 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -18.74799 15.31805 -1.224 0.221 speciesspruce 14.91628 16.99081 0.878 0.380 speciesother 25.64298 18.43308 1.391 0.164 log(height) 3.94011 3.13024 1.259 0.208 age6-10 3.63905 2.21270 1.645 0.100 age>10 0.08413 1.24922 0.067 0.946 speciesspruce:log(height) -3.07414 3.57427 -0.860 0.390 speciesother:log(height) -6.05681 3.99977 -1.514 0.130 (Dispersion parameter for binomial family taken to be 1) Null deviance: 41.589 on 29 degrees of freedom Residual deviance: 29.140 on 22 degrees of freedom AIC: 45.14 Number of Fisher Scoring iterations: 5 # EDIT 2 # ---------------------------------------------------- # Get fake model coefficients coefficients(mod) # Create a new model to update the coefficients: mod.real <- mod # Manually create vector of REAL coefficients to update the parameters: mod.real$coefficients<-c(-28.747,
14.916,
5.642,
2.940,
8.639,
0.084,
-5.074,
-6.056)

# create new data frame to predict values
df.new<-df

df.new$predicted<-predict.glm(mod.new,type="response")  • If I understand your question. Since it looks like your equation is correct, you can substitute in the coefficients in from the original model: mod$coefficients<-c(-28.747, 14.916, 5.642, 2.940, 8.639, 0.084, -5.074, -6.05681) and then use the predict function to calculate the odds: predict(mod, df) Commented Dec 16, 2019 at 21:14
• @dave2e in this case, should I just be carefull about the order of coefficients to don't mess up the specific estimates*parameters in a formula? Should I stransform the y from logaritms, or does the predict(mod, df) make it for me? Commented Dec 17, 2019 at 9:43
• Yes, the order of the coefficients do matter. And based on your equation the model will take the log of the height, so the just the raw linear value is the input. Commented Dec 17, 2019 at 13:30

Full model (with significant and non-significant parameters) is given in Suvanto et al. 2019. In detail it consists from 17 estimates (last 2 of them are second-order interactions). Although one solution to your question would be to use some internal R function, here I show you some internal details of how it really works (i.e. GLM with logit and interactions).

Suvanto, S., Peltoniemi, M., Tuominen, S., Strandström, M., & Lehtonen, A. (2019). High-resolution mapping of forest vulnerability to wind for disturbance-aware forestry and climate change adaptation. bioRxiv, 666305.

Note: you need this mini function logit2prob to convert value given by model equation (logit) to odds and finaly into probability (p).

Creating function according to model output:

library(car)
logit2prob <- function(logit){
odds <- exp(logit)
prob <- odds / (1 + odds)
return(prob)
}

# Function starts HERE...
forest.damage.function <- function(spec.spruce.BIN,
spec.other.BIN,
height.CONT,
last_thinning.6.10.BIN,
last_thinning.over.10.BIN,
wind.CONT,
open_stand_border.BIN,
soil_min.fine.BIN,
soil_organic.BIN,
soil_depth.more.30cm.BIN,
site_fertility.BIN,
temperature_sum.CONT,
damage_density.2to3.BIN,
damage_density.morethan3.BIN) {
temp.value = (-14.690) +
(- 8.494 * spec.spruce.BIN) +
(- 9.314 * spec.other.BIN) +
log(1.661 * height.CONT) +
(- 0.298 * last_thinning.6.10.BIN) +
(- 0.844 * last_thinning.over.10.BIN) +
log(0.749 * wind.CONT) +
(+ 0.310 * open_stand_border.BIN) +
(- 0.356 * soil_min.fine.BIN) +
(- 0.216 * soil_organic.BIN) +
(+ 0.214 * soil_depth.more.30cm.BIN) +
(- 0.425 * site_fertility.BIN) +
(+ 0.096 * temperature_sum.CONT) +
(+ 1.104 * damage_density.2to3.BIN) +
(+ 1.898 * damage_density.morethan3.BIN) +
(+ 1.634 * spec.spruce.BIN * log(height.CONT)) +
(+ 1.625 * spec.other.BIN * log(height.CONT))
final.value <- logit2prob(car::logit(temp.value))
return(final.value)
}
# ...and ends HERE.


Notice, that variables which enter the function ends either with .BIN or .CONT. This is for you to know what to insert. In case of .BIN use either 0 or 1. In cace of .CONT, use exact value of continual variable (f.e. temperature).

Another important thing. In model details you can find, that tree species could be given in 3 categories (pine, spruce, other), but as you can see only estimates for spruce and other care listed. This is not TRUE. Estimate for pine is in fact the value of Intercept. So in case that you are estimating for pine, the input looks like this:

forest.damage.function(spec.spruce.BIN = **0**,
spec.other.BIN = **0**,...etc.)


If you are estimating for spruce, the input looks like this:

forest.damage.function(spec.spruce.BIN = **1**,
spec.other.BIN = **0**,...etc.)


The same ("the Intercept thing" listed above) is true for all other categorical variables.

# example uf use...
forest.damage.function(spec.spruce.BIN = 1,
spec.other.BIN = 0,
height.CONT = 5,
last_thinning.6.10.BIN = 0,
last_thinning.over.10.BIN = 1,
wind.CONT = 35,
open_stand_border.BIN = 1,
soil_min.fine.BIN = 1,
soil_organic.BIN = 0,
soil_depth.more.30cm.BIN = 1,
site_fertility.BIN = 1,
temperature_sum.CONT = 350,
damage_density.2to3.BIN = 0,
damage_density.morethan3.BIN = 1)

• wow, @ladislavNado!! thank you for all the work! But, please, my new predictors will be in a dataframe in a long format, i.e. I not have categorical variables coded as 1 (present) or 0 (no present); ratehr they will be as df\$spruce <- c("spruce", "spruce", "pine", "other", "spruce",. etc..) How can I account for this in your approach? Commented Dec 19, 2019 at 11:27
• also, I guess we are applying log value of the variables, not of the coefficient? ie. not log(1.661 * height.CONT) but 1.661 * log(height.CONT)`? Commented Jan 7, 2020 at 15:10