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

Question

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")

Answer 1

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)