1
$\begingroup$

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.

enter image description here

$\endgroup$
  • 2
    $\begingroup$ I tried to replicate in R, but my actual and estimated coefficients lay along the $y=x$ line, as expected. $\endgroup$ – Matthew Drury Oct 1 '15 at 1:32
  • $\begingroup$ @MatthewDrury can you please post your R code? $\endgroup$ – DankMasterDan Oct 1 '15 at 5:35
  • $\begingroup$ Could you please explain your code? - not everyone reads Octave/Matlab. $\endgroup$ – Scortchi - Reinstate Monica Oct 1 '15 at 21:11
  • $\begingroup$ @Scortchi, I have added comments to the code. Let me know if you have any questions? $\endgroup$ – DankMasterDan Oct 1 '15 at 22:50
  • 2
    $\begingroup$ If your C = randn(D,1) is generating 10 data from a standard normal, why do you have a value >12? $\endgroup$ – gung - Reinstate Monica Oct 1 '15 at 23:38
1
$\begingroup$

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)

enter image description here

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

enter image description here

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.

$\endgroup$
  • $\begingroup$ The only difference in what I did what to use glm.fit on the matrix directly. I had essentially exactly the same results. Sorry I didn't get to posting my code before going to work this morning. $\endgroup$ – Matthew Drury Oct 1 '15 at 23:20
  • $\begingroup$ Hmmm, interesting. I've never really used glm.fit, @MatthewDrury. $\endgroup$ – gung - Reinstate Monica Oct 1 '15 at 23:22
  • $\begingroup$ @MatthewDrury, no problem. I can delete this, if you want to post your own code. $\endgroup$ – gung - Reinstate Monica Oct 1 '15 at 23:35
  • $\begingroup$ No, it's not worth the trouble. glm.fit is just called on the constructed model matrix from glm after the formula is applied. $\endgroup$ – Matthew Drury Oct 2 '15 at 1:04
  • $\begingroup$ Ok, so could this be an issue with Matlab? I just tried the same code on an older version of Matlab and got the same thing... If anyone has Matlab could they try it to see if they get the same thing? $\endgroup$ – DankMasterDan Oct 2 '15 at 18:51
1
$\begingroup$

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...

$\endgroup$
  • 1
    $\begingroup$ I almost commented on that, but figured it was a matlab quirk. I'm sorry, you could have had it solved earlier if i did. Lesson learned for me. $\endgroup$ – Matthew Drury Oct 2 '15 at 20:56
  • $\begingroup$ I was actually looking for answers to the same question posted, which I think is still valid independet of the code and solution to its bug you posted. The standardize flag in glmnet has the following text (in the package description): "The coefficients are always returned on the original scale." $\endgroup$ – hirschme Jul 22 '16 at 10:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.