You wrote "I normally use the means of the bins as predictor (independent) variables, but I want to use the variances of the bins too". Maybe it would be interesting for you to use quantile regression to obtain such a predictor?

See :

https://community.wolfram.com/groups/-/m/t/1395719

Hi!

Notice that "TotalVariation" and "Tikhonv" are options for Fit. But I wonder if there is a way to determine the regularization parameter?

I couldn't find these options with an online search here for LinearModelFit

Hmm. Not sure how this is relevant to a linear regression?

Also, even though I could write code to find out the estimate on the regression coefficient from the formulas given, I also want all the other statistical "goodies" such as p values, confidence intervals, etc. Do I need to write code for all of these myself?

Actually, in this specific case the variances are not technically errors, and they are not necessarily small. I have binned data (e.g. per country), but in each bin there is a distribution of data. I normally use the means of the bins as predictor (independent) variables, but I want to use the variances of the bins too (higher moments of the distribution are not available or useful at this point) in the regression analysis. And of course, as is standard, I have the variances and means for the predicted (dependent) variables. So if the newest version of Mathematica could handle this (some hints in the other responses that this may be possible), I would be motivated to upgrade (I have not upgrades since version 8.2). Otherwise I will try to switch to R, which I really don't want to do. Any word from Mathematica developers?

No, the variance in each bin is not a new predictor, but would combine with the variance in the dependent variables to give a weight in the chisquare sum (in the way which I indicated above, maybe there should be a square root). This seems intuitively true, but not sure how rigorous it is as far as obtaining true parameters from the limit of infinite data sets. Is there a Mathematica function or option to LinearModelFit that can do this?

just some comment in defense of the mathematicas:

1. stepping stones in augmenting power of statistics functions were V8 and V10. statistics functionality was heavily revised there, so you're missing out on the augmentation of V10.
2. by now there are even more statistics functions, see the Repository launched in summer 2019.
3. it may or may not be true that $R$ or $Maple$ have even more of those specialized statistics functions but there is some chance that none has the exact special function that you're looking for. If you do find the exact function there which you are looking for, then please share a link and we will understand exactly what you're looking for. You're request is too vague, i doht understand it completely without an example, sorry.
4. $\text{Wolfram L}$ provides you all the functions and vocab items you need to create your own specialized statistics functions (utility functions), tailored to your needs. So if you need to compose your own function (because none of the competition has it), then you're in the right programming environment. And if statistics is all you do, day in day out, then a specialized comprehensive GUI-based statistics application would be the better choice.

Keep up posted how your quest is going!

I'm going to try to revive this post. If there are known variances in dependent variable, their inverse can be used as a weight in computing the chi square sum in the regression analysis. But what if we also have information about the independent variables' variances? Can we use that information somehow? Naively, I would say (for one independent variable xi and one dependent variable yi) to find the regression coefficient r that minimizes chisquare=Sumi[(r xi -yi)^2/(r^2 sigmaxi^2+sigmay_i^2)]. Is this correct? And if so, how can the Mathematica stat package implement it (I suppose I could write my own code to do it, if not available in stat package)