# Normalizing for Population in Regression

I am currently interested in looking at the relationship between the number of universities in a state and an economic index by state. It occurs to me though that the number of colleges in a state may be correlated with population. After investigating the correlation, it was found to be 0.82.

If I divide the number of colleges by population in that state, my concern is if it make sense for interpretation. Every value for this newly constructed variable is roughly 0.000001. Is there another method I should consider for normalizing for population?

• You could rescale to colleges per million population which might make the numbers about 1 Aug 14, 2022 at 0:08
• "Every value for this newly constructed variable is roughly 0.000001." Why do you perceive this as a problem? Is it a) because the values are very small, or b) because the values are approximately equal? Aug 15, 2022 at 20:36
• It is mainly for making scatterplots and eventually running regressions. That is the only reason I am wondering if it is a problem, the fact that the numbers are so small for interpretation. Aug 15, 2022 at 20:46

The fact that a variable - call it $$N_i$$ - defined as 'number of universities in a state divided by population of state' takes very small values should not be a problem. It will mean that the absolute differences in values between states will be very small, but that does not matter. What matters for a regression are the proportional differences in values. Suppose for example that for two states $$j$$ and $$k$$, $$N_j=0.0000010$$ and $$N_k=0.0000012$$. Then $$N_k$$ is 20% larger than $$N_j$$. If instead you adopted Henry's suggestion, so that the values became $$1.0$$ and $$1.2$$, the proportional difference would still be 20%.

Suppose now that you estimate a regression equation:

$$N_i = \beta_0 + \beta_1 E_i + \epsilon_i$$

where $$E_i$$ is an economic index, the $$\beta$$'s are parameters to be estimated, and $$\epsilon_i$$ is the error term. (I'm assuming here that $$N_i$$ is taken to be the dependent variable.) The small values of the $$N_i$$ won't prevent you from running such a regression. They will just mean that the $$\beta$$'s will also be small (though how small will also depend on the values of $$E_i$$). As ever with regressions, further questions will be whether the estimated value of $$\beta_1$$ is significantly different from zero, and if so whether the estimated relationship can properly be given a causal interpretation.

• While there is no issue theoretically, I think in practice rounding problems might occur with such small numbers - think about numerical matrix inversion.
– BrsG
Aug 16, 2022 at 11:46
• @BrsG Agreed, this could be a problem with particular regression software depending on thr point at which it rounded small values. Aug 16, 2022 at 12:13

Theoretically it does not matter as @Adam Bailey showed (as long as you ignore computer theory).

Realistically, it matters quite a bit for at least two reasons:

• how do you interpret your coefficient? As a human being, do you prefer to hear that there are 1.2 Universities per 1 Million inhabitants in a city, or do you find 0.000001 per person easier to graps? Personally, I find the the latter rather dismal. Most of the times (for example google hospital beds per capita), the reported values will be scaled to something that humans can grasp.

• As @BrsG noted, computers may be your enemy here. They predominantly use floating point math (and if you deal with large numbers, integer arithemtic will also bite you). This may be a bit esoteric without examples but can be quickly shown with code. I answered a question asking if the "future value function is broken" some time ago which is related to this. It is probably better to use an actual regression example though.

I will use the simplest possible and just copy paste from Wikipedia's OLS example. One may think OLS will be OLS, no matter what implementation or language you use but that is not generally true. Least squares estimates satisfy the condition that the residual is orthogonal to the columns of $$\mathbf{X}$$. $$\mathbf{X^\prime (y - X\widehat{\beta})} = \mathbf{0}$$ which can be re-written as $$\mathbf{X^\prime X}\widehat{\mathbf{\beta}}=\mathbf{X^\prime y}$$ Solving for $$\beta$$ provides the standard text book formulas for least squares estimate calculated as $$\widehat{\mathbf{\beta}}=\mathbf{X^\prime X}^{-1}\mathbf{X^\prime y}.$$

I decided to use Julia because its a very flexible language that is intuitive to read, even if you are not familiar with the language itself.

# manually create the data from the Wikipedia example
using DataFrames, GLM, LinearAlgebra
height = [1.47,1.50,1.52,1.55,1.57,1.60,1.63,1.65,1.68,1.70,1.73,1.75,1.78,1.80,1.83]
weight = [52.21,53.12,54.48,55.84,57.20,58.57,59.93,61.29,63.11,64.47,66.28,68.10,69.92,72.19,74.46]

data = DataFrame(weight=weight,height=height, height2 = height.^2)


The simple Wiki example is like this:

If you dont want a package but prefer the good old manual way, you could use several procedures. The usual way to define a "division" operation for matrices is to define $$A / B = 𝐴𝐵^{−1}$$, where $$𝐵^{−1}$$ denotes the inverse of B. Because matrix multiplication is not commutative, one can also define a left division or so-called backslash-division as $$A \backslash B = 𝐴^{−1}𝐵$$. Moore Penrose Inverse is just Julia's pinv(X). Basically, OLS is solving Ax=b where you add a column of 1s if you want an intercept, so matrix decompositions. (a.k.a. matrix factorizations) of A\b is OLS. That's why the examples below define X to contain a vector of ones for the intercept, and computes the regression result as X\weight.

len = length("Normal Equation result 1")
println(len)
X= [ones(15) height height2] # adding the intercept

#b = (x'x)\x'y
println(lpad("Normal Equation result 1",len) *"$((X'X)^-1*X'*weight)") # using “normal equations” approach # which is the following println(lpad("Normal Equation result 2",len) *"$(inv(transpose(X)*X)*transpose(X)*weight)" )

β = X \ weight # Uses QR decomposition
println(lpad("QR result",len) *"$β") β = pinv(X) * weight # uses SVG println(lpad("SVD result",len) *"$β")

β = cholesky(Hermitian(X'X))\(X'weight) # Cholesky decomposition
println(lpad("Cholesky result",len) *"$$β") #println("$$(lpad(T,7)):")


So far so good and just what Wikipedia computes as well. Now suppose we measure height in cm (*100).

Still working well and as expected just $$/100$$ and $$/100^2$$ respectively.

Here is the problem

• If we get a bit fancier, and use nanometers, we get

Now, this does not look very good anymore. The official package shows everything as zero and struggles with SE and confidence intervals. QR and SVD clearly start to deteriorate and essentially you may be scratching your head if the results make sense.

• What definitely does not make sense is to have perfect multicollinearity. This is easy to show if we use height and height*2 (similar to the example in cross validated that I linked).

The official package is not exactly great in telling you that you cannot invert this matrix and simply excludes one of the two in computing by default. To be fair, the documentation is pretty clear about this though and if you include dropcollinear=false you get the warning that the matrix is not positive definite and that the Cholesky factorization failed. On the other hand, the manual formulas are giving you very explicit warnings. However, this only works if there is perfect multicollinearity (at machine precision).

Using the same logic, but making height2 very small (I chose ./100000000000000), you will see the following!

BANG! Completely different values depending on the implementation, but all provide a value despite there being perfect multicollinearity.

Explanation:

Floating-Point representation is by far the most common way of representing real numbers in computing. However, squeezing infinitely many real numbers into a finite number of bits requires an approximate representation because most calculations with real numbers will produce quantities that cannot be exactly represented.

The distance between two adjacent representable floating-point numbers is called machine epsilon.This distance is not constant, and the representable floating-point numbers are densest in the real number line near zero, and grow sparser exponentially as one moves farther away from zero.

theme(:juno)
range = -10.0^4:1:10.0^4
plot(range, [eps(i) for i in range],  legend = false, title =  "Machine epsilon for various values")
yaxis!("Machine Epsilon")
xaxis!("Values")


With large values (especially integers) you have a problem with overflow, which occurs when an arithmetic operation attempts to create a numeric value that is outside of the range that can be represented with a given number of digits – either higher than the maximum or lower than the minimum representable value. The most common result of an overflow is to wrap around the maximum. Computer arithmetic of integers is implemented that way. You can read some more here. This may sound crazy, but is exactly how car odometers worked, except that signed integers display negative values as well. It is also how the 12-hour clock works which wraps around every 12 hours (if it is 6:00 right now, 9 hours later it will be 3:00 and not 6+9 = 15).

Also, floating-point operations are not associative.

This post by Stefan Karpinski about array ordering and floating point math not being associative illustrates this very well. He wrote the below function to check how many different sums you can get by reordering a set of floating-point numbers. He defines a function called sumsto which is always returning the exact same 2046 floating-point numbers but returns them in a different order based on x. My version is slightly modified from the original code to suit the newer Julia version I am using.

realmax = reinterpret(Float64, reinterpret(UInt64, Inf)-1)
function sumsto(x::Float64)
0 <= x < exp2(970) || throw(ArgumentError("sum must be in [0,2^970)"))
n, p₀ = Base.decompose(x) # integers such that n*exp2(p₀) == x
[realmax; [exp2(p) for p in -1074:969 if iseven(n >> (p-p₀))]
-realmax; [exp2(p) for p in -1074:969 if isodd(n >> (p-p₀))]]
end


If you check out the same function with less values, say 13, you can see the pattern.

x = 0.0
println("For $$x the result is$$(foldl(+, sumsto(x)))")
println("For $$x the result is$$(foldl(+, sumsto(eps(x))))")
x = 1.23
println("For $$x the result is$$(foldl(+, sumsto(x)))")
x = ℯ
println("For $$x the result is$$(foldl(+, sumsto(x+0)))")
x=π
println("For $$x the result is$$(foldl(+, sumsto(x+0)))")
x = 6.0221409e23
println("For $$x the result is$$(foldl(+, sumsto(x)))")
x = 9.979201547673598e291
println("For $$x the result is$$(foldl(+, sumsto(x)))")


For any 64-bit float value of x from 0 up to (but not including) 2^970, the naive left-to-right sum of the vector returned by sumsto(x) is precisely x. To show it is not just some vodoo, I use the left-to-right sum with "normal" values. It clearly works here, but even with the 13 element vector from above, it produces the value 0.0 (although it is not always exact to the decimal in this case).

Summary

Long story short, I do not recommend to use very off the scale variables, especially if mixed together. The way computers operate is not very favourable in these cases.

Moreover, intuitively it makes more sense to me to read there are ~10 ICU beds per 100,000 as opposed to comparing 0.000083 beds per capita in Russia vs 0.0003870 in Germany.