Why does GLM/logistic regression rescale estimated parameters?

Question

I am trying to interpret my estimated paramter values with real data. In the process I made a very simple logistic regression example (GLM w/ logit link) where ground truth is available. I found that the estimated parameters were scaled versions of the ground truth parameters - furthermore the outlying probabilities (very close to 0/1) where poorly calibrated. Why is this & is there anyway to correct it?

Here is my Matlab code & example:

N = 100000; D = 10;   % N=# of observations. D=# of feautures
logit = @(x) 1./(1+exp(-x));    % defining logit/sigmoid link function
V = randn(N,D);   % randomly generated observation/design matrix from std. normal distr.
C = randn(D,1);   % randomly generated ground truth coefficients
S = logit(V*C);   % ground truth probability of firing
y = S>rand(length(N),1);   % observed binary sequence done by 'biased coin flipping'

Cp = glmfit(V,y,'binomial','link','logit','constant','off'); % Cp is est coef vector
yP = logit(V*Cp);     % esimted prob. of firing
clf; 
subplot(321); plot(C,Cp); 
subplot(312); I = [1:300]+1000; plot([S(I),yP(I)]);

And the produced figures: Notice that in the top plot of the true vs estimated paramters, the line passes through the origin (no shift) but has a slope different than 1, as indicated by the dashed red y=x line (ie there is a scale factor). Notice in the bottom plot the outlying probabilites are slightly off.

Answer 1

Like @MatthewDrury, I don't notice much amiss here when I run this in R:

set.seed(2568)  # this makes the example exactly reproducible

N = 100000; D = 10
logit = function(x){ 1/(1+exp(-x)) }
V = matrix(rnorm(N*D), nrow=N, ncol=D)
C = rnorm(D);  C  # these are the true parameter values:
# [1] -0.8242989  0.1945421 -0.4708643 -0.4535035 -0.4737520 -0.7723543
# [7] -0.6969630 -0.2676406 -0.7744012 -0.4767246
S = logit(V%*%C)  # NB, R uses "%*%" for matrix multiplication
y = as.numeric(S>runif(N))   # I convert the T & F into 1 & 0
d = data.frame(cbind(y, V))  # I put y & V into a data frame for convenience

m = glm(y~., data=d, family=binomial)
Cp = coef(m);  Cp  # these are the estimated parameter values
# (Intercept)           V2           V3           V4           V5           V6 
# 0.004374693 -0.828212861  0.194661608 -0.472036943 -0.450215381 -0.487840814 
#           V7           V8           V9          V10          V11 
# -0.768954798 -0.694885635 -0.286258326 -0.766237973 -0.485269940
yP = logit(cbind(1,V)%*%Cp)

windows()
  plot(c(0,C), Cp)  # I add a 0 at the front for the absent intercept
  abline(h=0, col="blue", lty=3)
  abline(v=0, col="blue", lty=3)
  abline(0, 1, col="red", lty=3)

p = cbind(S, yP)[order(S),]  # I sort the probability values for a cleaner figure
windows()
  plot(p[,1], p[,2])

Edit: Perusing glmfit, I see that the last two arguments in your code were to suppress the intercept. I reran the above with m = glm(y~.-1, data=d, family=binomial), but get the same good results.

Answer 2

Slightly embarrasing, but I just realized my bug. I had:

y = S>rand(length(N),1);

But here N=10000 is a scalar not a vector, so length(N)=1. Essentially I was giving the signal a hard threshold. The correct code is:

y = S>rand(N,1);

Thank you to everyone for their help & sorry to have wasted your time. Not sure what the stachexchange policy is, but I will be glad to delete the question if that is appropriate...