The standard errors of your estimated coefficients
The standard errors of your estimated coefficients, $\hat\beta$, are estimated according to:
$$\widehat{s.e.}(\hat\beta_j
) = \sqrt{s^2(X^TX)_{jj}^{-1}}$$
where $s^2$ is an estimate of the deviation/variance of the noise and $X$ is the matrix of your regressors.
In other words the standard error of the coefficient $\beta_j$ scales with the inverse of the root mean square of the column/regressor $X_j$. (you should be able to verify this easily which helps you with debugging)
Note that these are marginal errors and the formula can be generalized to provide a covariance matrix of the variance of the estimated coefficients. If there is correlation then the marginal errors might be very large but will not reflect the distribution of error correctly
Example of direct (manual) computation of standard error of coefficients
Computational example (in R, but I hope that the annotations make it clear for you to turn it into python code):
###
### creating data
###
set.seed(1)
### matrix with three columns
x <- cbind(c(1,1,1,1,1,1),
c(1,2,3,4,5,6),
c(0,0,0,1,1,1))
### some output based on the second column
### with added normal distributed noise with sigma=0.1
y <- x[,2] + rnorm(6,0,0.1)
### computing the model with the standard lm function
mod <- lm(y~0+x)
summary(mod)
### computing the std.err manually based on
### sd of the residuals
### the inverse (using the function solve) of the matrix product t(x)*x
### a correction sqrt(5/3) because the residuals should be based on n-p degrees of freedom and not n-1
### (and we have p=3, three parameters)
sd(residuals(mod))*sqrt(solve(t(x) %*% x)[c(1,5,9)])*sqrt(5/3)
output lm function:
> summary(mod)
Call:
lm(formula = y ~ 0 + x)
Residuals:
1 2 3 4 5 6
-0.08565 0.06098 0.02467 0.05709 -0.00386 -0.05323
Coefficients:
Estimate Std. Error t value Pr(>|t|)
x1 0.08863 0.08888 0.997 0.392139
x2 0.93438 0.03848 24.279 0.000153 ***
x3 0.27629 0.13145 2.102 0.126334
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.07697 on 3 degrees of freedom
Multiple R-squared: 0.9998, Adjusted R-squared: 0.9996
F-statistic: 5126 on 3 and 3 DF, p-value: 4.625e-06
output manual:
> sd(residuals(mod))*sqrt(solve(t(x) %*% x)[c(1,5,9)])*sqrt(5/3)
[1] 0.08887551 0.03848422 0.13144865
This above code shows where those standard errors 0.089, 0.038, 0.131 come from.
Check the root mean square size of your regressors
So in your situation you probably have some very large (root mean squared) size for several of the columns/regressors $\vert X_j \vert$. And that is why the parameters (and the estimated error) have such low values.
(I would also suspect that there might be some computational errors, however this is difficult to inspect with the given information. Anyway, at least you can now verify with this answer whether it is the standard case, ie. the relation with $X^TX$)
In addition: Possibly you might want to do a mixed effects model? Those fixed effects per user do they make sense to be considered fixed effects, or can we more specifically regard them as following a normal distribution?