Edit: I added what I hope is a real answer (bad pun intended).
The linear model needs to be made explicit as there can be models that might look similar but are really parameterized very differently and would imply possibly very different analysis procedures and interpretation of results.
Using the notation from @whuber 's answer from CrossValidated, if we make the response variable ($z=x+i y$), predictor variable ($w=u+i v$), slope (b=$\gamma_1+i \delta_1$), intercept (a=$\gamma_0+i \delta_0$), and error ($\epsilon_r+i \epsilon_i$) all complex variables, then a simple linear model might look like the following:
$$z = a + b w + \epsilon$$
$$x + i y = \gamma_0 + i \delta_0 + (\gamma_1 + i \delta_1)*(y + i v) + \epsilon_r + i \epsilon_i$$
where $\epsilon_r$ and $\epsilon_i$ have a bivariate normal distribution.
Here is a Mathematica function using maximum likelihood to estimate the parameters:
(* X is the design matrix *)
(* z in the response variable *)
(* equalVariances is True if you want the real and imaginary components to have the same variance and False otherwise. *)
(* corr0 is True if you want to assume that the correlation between the real and imaginary error components to be zero.
Otherwise, the correlation will be estimated from the data. *)
complexLinearModelFit[{X_, z_}, equalVariances_, corr0_] := Module[
{model, modelR, modelI, logL, vR, vI, v, lmR, lmI, estimates,
initialValues, ρ0, σR0, σI0, mle, cov, aic, γ, δ, xR, xI, conditions},
(* Regression parameters *)
γ = Table[ToExpression["γ" <> ToString[i]], {i, Dimensions[X][[2]]}];
δ = Table[ToExpression["δ" <> ToString[i]], {i, Dimensions[X][[2]]}];
model = Sum[(γ[[j]] + I δ[[j]]) X[[All, j]], {j, Dimensions[X][[2]]}];
modelR = Re[ComplexExpand[model]] /. Im[h_] -> 0 /. Re[h_] -> h;
modelI = Im[ComplexExpand[model]] /. Im[h_] -> 0 /. Re[h_] -> h;
(* Get lists of variables *)
vR = Variables[modelR];
vI = Variables[modelI];
v = Variables[model];
(* Get regression coefficients associated with each variable *)
xR = Table[Coefficient[modelR, vR[[i]]], {i, Length[vR]}];
xI = Table[Coefficient[modelI, vI[[i]]], {i, Length[vI]}];
(* Perform linear models on real and imaginary components to get initial estimates of coefficients *)
lmR = LinearModelFit[{Transpose[xR], Re[z]}];
lmI = LinearModelFit[{Transpose[xI], Im[z]}];
(* Get the mean of parameter estimates estimate *)
estimates = Transpose[{Join[vR, vI], Join[lmR["BestFitParameters"], lmI["BestFitParameters"]]}];
initialValues = Table[Mean[Select[estimates, #[[1]] == v[[i]] &][[All, 2]]], {i, Length[v]}];
(* Now get initial values for the covariance matrix *)
ρ0 = Correlation[lmR["FitResiduals"], lmI["FitResiduals"]];
σR0 = lmR["EstimatedVariance"]^0.5;
σI0 = lmI["EstimatedVariance"]^0.5;
(* Put together all of the initial estimates,
log likelihood function, and determine conditions on parameters *)
logL = LogLikelihood[
BinormalDistribution[{0, 0}, {σR, σI}, ρ],
Transpose[{Re[z] - modelR, Im[z] - modelI}]];
If[equalVariances,
logL = logL /. σR -> σ /. σI -> σ;
conditions = {σ > 0};
v = Join[v, {σ}];
initialValues = Join[initialValues, {(σR0 + σI0)/2}],
conditions = {σR >= 0, σI >= 0};
v = Join[v, {σR, σI}];
initialValues = Join[initialValues, {σR0, σI0}]
];
If[corr0,
logL = logL /. ρ -> 0,
conditions = Join[conditions, {-1 <= ρ <= 1}];
v = Join[v, {ρ}];
initialValues = Join[initialValues, {ρ0}]
];
initialValues = Transpose[{v, initialValues}];
(* Find maximum likelihood estimates *)
mle = FindMaximum[{logL, conditions}, initialValues];
(* Estimates of standard errors *)
(cov = -Inverse[(D[logL, {v, 2}]) /. mle[[2]]]);
(* AIC *)
aic = -2 mle[[1]] + 2*Length[v];
(* Return results from FindMaximum, covariance matric estimate, and AIC *)
{mle, cov, aic}
]
Now for some examples. First consider the data posted:
data = {{1 + I, 0.985402 + 1.08528 I}, {1 + 2 I, 2.09444 + 1.00236 I},
{1 + 3 I, 3.29011 + 0.969815 I}, {1 + 4 I, 3.98937 + 0.969446 I},
{1 + 5 I, 4.71475 + 0.916196 I}, {1 + 6 I, 5.97777 + 0.994892 I},
{1 + 7 I, 6.35507 + 1.02171 I}, {1 + 8 I, 7.41285 + 0.9147 I},
{1 + 9 I, 8.73952 + 1.04088 I}, {1 + 10 I, 9.36015 + 0.947539 I}};
response = data[[All, 2]];
X = Transpose[{ConstantArray[1, Length[data]], data[[All, 1]]}];
lm = complexLinearModelFit[{X, response}, False, False];
(* Parameter estimates *)
clm[[1, 2]]
(* {γ1 -> 0.269328, γ2 -> -0.00612192, δ1 -> 1.93427, δ2 -> -0.914316,
σR -> 0.193555, σI -> 0.0483416, ρ -> 0.0186871} *)
(* Covariance matrix *)
clm[[2]] // TableForm

(* AIC *)
clm[[3]]
(* -22.6791 *)
Now for the example given on CrossValidated:
(* Predictor variable *)
w = {0 - 5 I, -3 - 4 I, -2 - 4 I, -1 - 4 I, 0 - 4 I, 1 - 4 I, 2 - 4 I, 3 - 4 I, -4 - 3 I, -3 - 3 I, -2 - 3 I, -1 - 3 I, 0 - 3 I, 1 - 3 I, 2 - 3 I, 3 - 3 I, 4 - 3 I, -4 - 2 I, -3 - 2 I, -2 - 2 I, -1 - 2 I, 0 - 2 I, 1 - 2 I, 2 - 2 I, 3 - 2 I, 4 - 2 I, -4 - 1 I, -3 - 1 I, -2 - 1 I, -1 - 1 I, 0 - 1 I, 1 - 1 I, 2 - 1 I, 3 - 1 I, 4 - 1 I, -5 + 0 I, -4 + 0 I, -3 + 0 I, -2 + 0 I, -1 + 0 I, 0 + 0 I, 1 + 0 I, 2 + 0 I, 3 + 0 I, 4 + 0 I, 5 + 0 I, -4 + 1 I, -3 + 1 I, -2 + 1 I, -1 + 1 I, 0 + 1 I, 1 + 1 I, 2 + 1 I, 3 + 1 I, 4 + 1 I, -4 + 2 I, -3 + 2 I, -2 + 2 I, -1 + 2 I, 0 + 2 I, 1 + 2 I, 2 + 2 I, 3 + 2 I, 4 + 2 I, -4 + 3 I, -3 + 3 I, -2 + 3 I, -1 + 3 I, 0 + 3 I, 1 + 3 I, 2 + 3 I, 3 + 3 I, 4 + 3 I, -3 + 4 I, -2 + 4 I, -1 + 4 I, 0 + 4 I, 1 + 4 I, 2 + 4 I, 3 + 4 I, 0 + 5 I};
(* Add in a "1" for the intercept *)
w1 = Transpose[{ConstantArray[1 + 0 I, Length[w]], w}];
z = {-15.83651 + 7.23001 I, -13.45474 + 4.70158 I, -13.63353 + 4.84748 I, -14.79109 + 4.33689 I, -13.63202 + 9.75805 I, -16.42506 + 9.54179 I, -14.54613 + 12.53215 I, -13.55975 + 14.91680 I, -12.64551 + 2.56503 I, -13.55825 + 4.44933 I, -11.28259 + 5.81240 I, -14.14497 + 7.18378 I, -13.45621 + 9.51873 I, -16.21694 + 8.62619 I, -14.95755 + 13.24094 I, -17.74017 + 10.32501 I, -17.23451 + 13.75955 I, -14.31768 + 1.82437 I, -13.68003 + 3.50632 I, -14.72750 + 5.13178 I, -15.00054 + 6.13389 I, -19.85013 + 6.36008 I, -19.79806 + 6.70061 I, -14.87031 + 11.41705 I, -21.51244 + 9.99690 I, -18.78360 + 14.47913 I, -15.19441 + 0.49289 I, -17.26867 + 3.65427 I, -16.34927 + 3.75119 I, -18.58678 + 2.38690 I, -20.11586 + 2.69634 I, -22.05726 + 6.01176 I, -22.94071 + 7.75243 I, -28.01594 + 3.21750 I, -24.60006 + 8.46907 I, -16.78006 - 2.66809 I, -18.23789 - 1.90286 I, -20.28243 + 0.47875 I, -18.37027 + 2.46888 I, -21.29372 + 3.40504 I, -19.80125 + 5.76661 I, -21.28269 + 5.57369 I, -22.05546 + 7.37060 I, -18.92492 + 10.18391 I, -18.13950 + 12.51550 I, -22.34471 + 10.37145 I, -15.05198 + 2.45401 I, -19.34279 - 0.23179 I, -17.37708 + 1.29222 I, -21.34378 - 0.00729 I, -20.84346 + 4.99178 I, -18.01642 + 10.78440 I, -23.08955 + 9.22452 I, -23.21163 + 7.69873 I, -26.54236 + 8.53687 I, -16.19653 - 0.36781 I, -23.49027 - 2.47554 I, -21.39397 - 0.05865 I, -20.02732 + 4.10250 I, -18.14814 + 7.36346 I, -23.70820 + 5.27508 I, -25.31022 + 4.32939 I, -24.04835 + 7.83235 I, -26.43708 + 6.19259 I, -21.58159 - 0.96734 I, -21.15339 - 1.06770 I, -21.88608 - 1.66252 I, -22.26280 + 4.00421 I, -22.37417 + 4.71425 I, -27.54631 + 4.83841 I, -24.39734 + 6.47424 I, -30.37850 + 4.07676 I, -30.30331 + 5.41201 I, -28.99194 - 8.45105 I, -24.05801 + 0.35091 I, -24.43580 - 0.69305 I, -29.71399 - 2.71735 I, -26.30489 + 4.93457 I, -27.16450 + 2.63608 I, -23.40265 + 8.76427 I, -29.56214 - 2.69087 I};
(* Estimation assuming ρ=0 (which is what the CrossValidated example assumes *)
complexLinearModelFit[{w1, z}, True, True][[1, 2]]
(* {γ1 -> -20.0172, γ2 -> -0.830797, δ1 -> 5.00968, δ2 -> 1.37827, σ -> 2.20038} *)
(* Now allow the estimation of ρ *)
complexLinearModelFit[{w1, z}, True, False][[1, 2]]
(* {γ1 -> -20.0172, γ2 -> -0.763237, δ1 -> 5.00968, δ2 -> 1.30859, σ -> 2.21424, ρ -> 0.810525} *)
The true values are γ1 -> -20, γ2 -> -0.75, δ1 -> 5, δ2 -> 1.299038, σ -> 2, and ρ -> 0.8 which are very close to the estimates when $\rho$ is allowed to be estimated.
The code should work for any number of complex predictors. But, again, that assumes that the models being fit with the code are the models that you need.