Upsampling with time offsets
Answer: You will see residual images of $X(f)$ at multiples $f_s$, $2f_s$ and $3f_s$, and distorted image of $X(f)$ at non-zero multiples of $4f_s$, when sampling in the manner you explained. Depending on value $e$, the size of residual will change. I have explained how in detail below. 

Ideally, sampling at $4f_s$ would have completely cancelled those images of $X(f)$ at multiples of $f_s$, $2f_s$ and $3f_s$ and you would've seen images of $X(f)$ only at multiples of $4f_s$ and the magnitude scaled by $4f_s$. 

(Explanation seems long only because I have included lots of pictures to show. Please follow through.)

I would like to give you an intuition on how to visualize sampling at any rate $f_s$. You probably have a pretty good idea about that. But then I would like to show pictorially what happens when you sample at $2f_s$ and then you can extend the idea to $4f_s$.




  
  Sampling a bandlimited signal $x(t)$ at sampling rate $f_s$ :
  




When you sample $x(t)$ at sampling rate $f_s$, you are basically, multiplying $x(t)$ with a periodic pulse train with period $T_s = \frac{1}{f_s}$ in time domain. Hence, in frequency domain you see a convolution of $X(f)$ with Fourier representation of that periodic pulse train.

Since the pulse train is periodic, it's Fourier representation will be obtained by computing Fourier Series. The pulse train can be represented by: $$\sum^{\infty}_{k=-\infty} \delta(t - kT_s)$$

 and it's Fourier transform as: $$\frac{1}{T_s}\sum^{\infty}_{k=-\infty} \delta(f - kf_s)$$

Notice that the magnitude of Fourier domain representation of the sampling pulse train has scaled up magnitude of $\frac{1}{T_s} = f_s$.


Assume, frequency representation of $x(t)$ as below :


Now, as described above, Sampling is nothing but multiplication of $x(t)$ with a $T_s$ periodic pulse train in time domain and hence a convolution of $X(f)$ with Fourier transform of the pulse train in frequency domain. Mathematically,
$$x(kT_s) = x(t).\sum^{\infty}_{k=-\infty}\delta(t-kT_s)$$
$$X_{sampled}(f) = X(f) * f_s. \sum^{\infty}_{k=-\infty}\delta(f-kf_s)$$





  
  Now, lets explore what happens when we sample at $2f_s$ or when the sampling pulse train becomes $\frac{T_s}{2}$ periodic :
  




The sampling pulse train becomes as shown below :


The sampling pulse train can be represented by a sum of two $T_s$ periodic pulse trains:


One pulse train as if $T_s$ periodic and centered at 0.
Second pulse train as if $T_s$ periodic but shifted by $\frac{Ts}{2}$.


Hence, mathematically it will be as follows:
$$\sum^{\infty}_{k=-\infty}\delta(t-k\frac{T_s}{2}) = \sum^{\infty}_{k=-\infty}\delta(t-kT_s) + \sum^{\infty}_{k=-\infty}\delta(t-kT_s - \frac{T_s}{2})$$


Hence, the Fourier Transform of the sampling pulse train will also sum of these 2 trains, as convolution is a linear operation. Also, use the time shift property of Fourier transform to get the result as follows:

$$\mathcal F \{ \sum^{\infty}_{k=-\infty}\delta(t-kT_s) + \sum^{\infty}_{k=-\infty}\delta(t-kT_s - \frac{T_s}{2}) \}$$ 
$$= f_s. \sum^{\infty}_{k=-\infty}\delta(f-kf_s) + f_s.e^{-j\pi \frac{f}{f_s}} \sum^{\infty}_{k=-\infty}\delta(f-kf_s)$$

$$= f_s. \sum^{\infty}_{k=-\infty}\delta(f-kf_s) + f_s.(cos(\pi\frac{f}{f_s})-\mathbb i.sin(\pi \frac{f}{f_s})). \sum^{\infty}_{k=-\infty}\delta(f-kf_s)$$

Notice that the sum is evaluated only at integral multiples of $f_s$, because of the $\delta(f - kf_s)$. What this means is that $sin(\pi \frac{f}{f_s})$ will always be $0$, so, no imaginary images of $X(f)$ will be seen, and $cos(\pi \frac{f}{f_s})$ will be $1$ at even multiples of $f_s$ and $-1$ at odd multiples of $f_s$. Pictorially, the fourier transform of pulse train which is used to sample at $2f_s$ will look like following :


So, the even multiples of $f_s$ will be doubled in magnitude to $2f_s$ and odd multiples of $f_s$ will cancel out each other to cancel the images of $X(f)$. 
This is the reason you see images of $X(f)$ only at multiples of $2f_s$ when sampling at double the rate, because images of $X(f)$ at odd multiples of $f_s$ cancel each other out.

Now, consider the case which you have explained in your question. When you break your pulse train into 4 pulse trains which are $T_s$ periodic individually, but shifted as below:


$\sum^{\infty}_{k=-\infty}\delta(t-kT_s)$
$\sum^{\infty}_{k=-\infty}\delta(t-kT_s -\frac{T_s}{4} - e)$
$\sum^{\infty}_{k=-\infty}\delta(t-kT_s -\frac{2T_s}{4} - 2e)$
$\sum^{\infty}_{k=-\infty}\delta(t-kT_s -\frac{3T_s}{4} - 3e)$


Conclusion:

When you check their Fourier transforms, you will find that images at multiples of $f_s$, $2f_s$ and $3f_s$ will not get cancelled completely because the negative impulses (both real and imaginary) are shifted by $e$, $2e$ and $3e$ respectively. And, image at multiples of $4f_s$ will also not be aligned exactly to give a scaling of $4f_s$ but they will be fudged around to give a distorted image of $X(f)$ except at $k=0$, that is the original image of $X(f)$ centered around DC.

Depending upon the value of e, the real and imaginary images of $X(f)$ will have residuals at $f_s$, $2f_s$ and $3f_s$, and images at non-zero multiples of $4f_s$ will be fudged around.