At the distance $r = R$ the electric field is $E = Q / 8 \pi \epsilon_0 R^2$, which is half way between the field values immediately to either side. As you say this is a 'mathematical shell', and this surface field value is calculated within the mathematical model that comes out of the physical model of electrostatics, but this physical model reaches its limits on considering a surface of zero thickness - it would be unlikely to be possible to measure such a surface field value in an experiment or to detect it in any way.
One way to calculate this value of $E$ is given towards the end of the answer [PSE1]. The answer [PSE2] describes an issue related to your question of how Gauss's Law no longer applies when a Gaussian surface contains point charges or a surface charge distribution - in this case the Generalized Gauss Theorem applies instead. In the latter theorem a factor of $1/2$ appears for charges on the surface. Volume distributions of charge don't affect the Gaussian surface in the same way since on a surface of zero thickness they produce no charge.
A second method of calculating $E$ at $r = R$ is to carry out the integration manually using the spherical formula for surface integrals described in [MSE]. This formula applies to a surface $S$ which does not intersect the $z$-axis.
Consider a point $P$ on the spherical shell of charge $S$. Arrange axes so that $P$ is at $(0, 0, R)$. Let $\sigma = Q / 4 \pi R^2$ be the surface charge density. The contribution to electric field at $P$ by the charge element at $P$ will be infinite. However we can arrive at a meaningful mathematical value for $E$ at $P$ by using a limiting process. Cut out a small circle of $S$ centered at $P$, of polar angular radius $\delta$, to leave an approximation $S_{\delta}$ to $S$. Then mathematically the electric field $E$ at $P$ due to $S$ can be considered to be the limit of the field $E_{\delta}$ at $P$ due to $S_{\delta}$ as $\delta \rightarrow 0^+$. By symmetry, $E_{\delta}$ is in the radial direction as the tangential components of the field cancel. We will also cut out a corresponding circle at the diametrically opposite side of the sphere from $P$ so $S_{\delta}$ does not intersect the $z$-axis and the spherical formula [MSE] applies.
If $\underline{\mathbf{\hat{r}}} = \sin \theta \cos \varphi \: \underline{\mathbf{i}} + \sin \theta \sin \varphi \: \underline{\mathbf{j}} + \cos \theta \: \underline{\mathbf{k}}$ is the radial unit vector at general position $(R, \theta, \varphi)$ on $S$ (using Physics convention [SPH] for spherical coordinates $(r, \theta, \varphi)$), then since the radial direction at $P$ is in the $\underline{\mathbf{k}}$ direction we have:
$$
E_{\delta} = \int_{S_{\delta}} \frac{1}{4 \pi \epsilon_0} \cdot \frac{-R \; \underline{\mathbf{\hat{r}}} + R \; \underline{\mathbf{k}}}{|-R \; \underline{\mathbf{\hat{r}}} + R \; \underline{\mathbf{k}}|^3} \cdot \underline{\mathbf{k}} \; \sigma \; dS = \int_{S_{\delta}} \frac{1}{4 \pi \epsilon_0 R^2} \cdot \frac{-\underline{\mathbf{\hat{r}}} + \underline{\mathbf{k}}}{|-\underline{\mathbf{\hat{r}}} + \underline{\mathbf{k}}|^3} \cdot \underline{\mathbf{k}} \; \sigma \; dS
$$
noting the denominator does not go to zero as the point $P$ is avoided on the surface $S_{\delta}$.
But $|-\underline{\mathbf{\hat{r}}} + \underline{\mathbf{k}}|^2 = \sin^2 \theta + (1 - \cos \theta)^2 = 2 - 2 \cos \theta$ and $(-\underline{\mathbf{\hat{r}}} + \underline{\mathbf{k}}) \cdot \underline{\mathbf{k}} = 1 - \cos \theta$, so:
$$
E_{\delta} = \int_{S_{\delta}} \frac{\sigma}{8 \sqrt{2} \pi \epsilon_0 R^2} \cdot \frac{1}{(1 - \cos \theta)^{1/2}} \; dS \hspace{2em} \mbox{(surface integral over $S_{\delta}$)}
$$
Using the formula for surface integration in spherical coordinates described in [MSE], with $r(\theta, \varphi) = R, \; \forall \; (\theta, \varphi) \in D_{\delta} = [\delta, \pi - \delta] \times [0, 2 \pi]$ and $\displaystyle \frac{\partial r}{\partial \theta} = \frac{\partial r}{\partial \varphi} = 0$ throughout $S_{\delta}$, we then have :
\begin{eqnarray*}
E_{\delta} & = & \int_{D_{\delta}} \frac{\sigma}{8 \sqrt{2} \pi \epsilon_0 R^2} \cdot \frac{R^2 \sin \theta}{(1 - \cos \theta)^{1/2}} \; dA \hspace{2em} \mbox{(area integral over region $D_{\delta} \subseteq \mathbb{R}^2$)} \\
& = & \int_{\delta}^{\pi - \delta} \int_0^{2\pi} \frac{\sigma}{8 \sqrt{2} \pi \epsilon_0} \cdot \frac{\sin \theta}{(1 - \cos \theta)^{1/2}} \; d\varphi \; d\theta \hspace{2em} \mbox{(converting to iterated integral)} \\
& = & \int_{\delta}^{\pi - \delta} \frac{\sigma}{4 \sqrt{2} \epsilon_0} \cdot \frac{\sin \theta}{(1 - \cos \theta)^{1/2}} \; d\theta = \frac{\sigma}{4 \sqrt{2} \epsilon_0} \left[ 2(1 - \cos \theta)^{1/2} \right]_{\delta}^{\pi - \delta} \\
\mbox{thus} \hspace{2em} E & = & \lim_{\delta \rightarrow 0^+} E_{\delta} = \frac{\sigma}{2 \epsilon_0} = \frac{Q}{8 \pi \epsilon_0 R^2}, \hspace{1em} \mbox{as required}.
\end{eqnarray*}
References
[PSE1] Physics StackExchange Answer : Using Gauss's law when point charges lie exactly on the Gaussian surface, https://physics.stackexchange.com/a/805337/111652
[PSE2] Physics StackExchange Answer : Using Gauss's law when point charges lie exactly on the Gaussian surface, https://physics.stackexchange.com/a/544481/111652
[MSE] Maths StackExchange Answer : Surface integrals in spherical coordinates, https://math.stackexchange.com/a/4867164/417024
[SPH] Spherical coordinate system, https://en.wikipedia.org/wiki/Spherical_coordinate_system
