So sorry, I thought it was an identity but I typed lhs in place of rhs in one place. I probably wouldn't have answered otherwise. Now that I have, let me say that the linear fit does seem good. I can't say how unexpected that is, but perhaps not very. Here is a stab at an explanation along the lines of my previous silly answer. The short version is that the $y,z$ above, for the range given, are not *that* far off $\frac{4}{a+b+c}$ and $\frac{3}{a+b+c}$ so it is perhaps not amazing that $z \approx0.75y.$

You have $y=\frac{C(a)+C(b)+C(c)+C(a)C(b)C(c)}{S(a)+S(b)+S(c)+S(a)S(b)S(c)}$ and $z=\frac{C((a+b+c)/3)}{S((a+b+c)/3)}$ for $C(u)=\cos(u)$ and $S(u)=\sin(u)$ over the range $\frac{\pi}{18} \le u \le \frac{2\pi}{9} \lt 0.7.$

If we instead use the linear approximations $C(u)=1$ and $S(u)=u$ we are comparing $y=\frac{3+1}{a+b+c+abc}$ to $z=\frac{3}{a+b+c}.$ If you run your random experiments on those two I think that you will find the correlation coefficients similar and the residuals smaller. Those linear approximations are getting kind of off near the upper end of the range but at the middle $\theta=5\pi/36\approx 0.436,$ we have $\cos(\theta) \approx 0.906$ and $\sin(\theta) \approx 0.423.$

Keep that $S$ but use instead $C(u)=1-\frac{u^2}{2}$ to get $y=\frac{4-(a^2+b^2+c^2)+(a^2b^2+a^2c^2+b^2c^2)/4-a^2b^2c^2/8}{a+b+b+abc} \approx \frac{4-(a^2+b^2+c^2)}{a+b+b+abc}$ and $z=\frac{3-(a+b+c)^2/6}{a+b+c}.$ With those the correlation coefficients are again about the same and the residuals larger than in the linear case but smaller than for $\cos,\sin$

Since I worked it out: with $S(u)=u-u^3/6$, $C(u)=1-u^2/2$
$\frac{C(a)+C(b)+C(c)+C(a)C(b)C(c)}{S(a)+S(b)+S(c)+S(a)S(b)S(c)}=\frac{27(32-8(a^2+b^2+c^2)+2(a^2b^2+a^2c^2+b^2c^2)-a^2b^2c^2)}{216(a+b+c)+36(6abc-a^3-b^3-c^3)-36abc(a^2+b^2+c^2)+6(a^3+b^3+c^3)(abc)-a^3b^3c^3} $

and $\frac{C((a+b+c)/3)}{S((a+b+c)/3)}=\frac{9(18-(a+b+c)^2)}{(a+b+c)(54-(a+b+c)^2)}.$

I'll suggest that the higher order terms I dropped are not so significant both because the variables (with average value under $0.5$) occur to higher powers and because the coefficients are smaller. Higher order Taylor polynomials would have used $u^4/24$ and $u^5/120.$

I stuck to the form given although I might have preferred cross multiplying to have polynomials. It might be better to use polynomials centered at $u=\frac{5\pi}{36}$ which seem to look like $C(u) \approx 0.9986+0.0131u-0.5454u^2+0.0704u^4$ and $S(u)\approx -0.0005+1.0044u-0.0136u^2-0.1511u^3.$ That I leave to anyone who cares to check it out.

I'd also prefer to avoid the singularity at $0$ and work with the reciprocals. Then for $y=\frac{\sin a+\sin b +\sin c +\sin a\sin b \sin c}{\cos a+\cos b +\cos c +\cos a\cos b \cos c}$ and $z=\tan \frac{a+b+c}{3}$ the best linear fit for experiments like yours is $z \approx 1.021y+0.067.$

In that spirit I'll mentions that when $a=b=c=x$ we are comparing (the reciprocals) of $\tan(x)$ and $\frac{3\sin(x)+\sin^3(x)}{3\cos(x)+\cos^3(x)}=\tan(x)\frac{3+4\tan^2(x)}{4+3\tan^2(x)}.$ Those are equal at $x=0$ and $x=\frac{\pi}{4}.$ The derivatives differ more at those points than for the internal subinterval $(\frac{\pi}{18},\frac{2\pi}{9}).$ The exactly center interval of width $\frac{\pi}{6}$ is $(\frac{\pi}{24},\frac{5\pi}{24}).$ Perhaps the correlation would be even better, especially centered at $\frac{\pi}{8}$. This is only a downward shift of $\frac{\pi}{72}$ so it may not make much difference. Again, I don't plan to check.