partial fractions expansion inverse Z-transform, help
First convert your input-output relation to Z-domain to see:

$$ Y(z) = 0.25 \left( 1 + z^{-1} + z^{-2} + z^{-3} \right) X(z)$$

From which you find the transfer function of the LTI system as:

$$ H(z) = \frac{Y(z)}{X(z)} = 0.25 \left( 1 + z^{-1} + z^{-2} + z^{-3} \right) $$

Now this is an FIR system with 4 taps. As MattL. tried (without luck) to convince you that for such a system a partial fraction expansion is not the proper way to compute the output. Indeed for any FIR system, the transfer function is of the form of a $H(z) = \sum h[k]z^{-k}$ and for such a system the output is defined through a convolution as the following: 

$$ Y(z) = H(z) X(z) =  \left( \sum h[k]z^{-k} \right) X(z) \implies y[n] = \sum h[k]x[n-k] $$

so for your example the proper way to describe the output for the input $x[n]= u[n]$ is

$$ \boxed{ y[n] = 0.25 ~( u[n] + u[n-1] + u[n-2] + u[n-3]) } $$

Nevertheless, assuming you want to insist on a partial fraction expansion to compute the output $y[n]$ then you would do the following:

Given $x[n] = u[n]$ and its Z-transform $X(z) = \frac{1}{1-z^{-1}} $, then the output is:

$$Y(z) = \frac{1}{4} ~\frac{ 1 + z^{-1} + z^{-2} + z^{-3} } {1 -z^{-1} } $$

Now, ignoring the linear scale factor $1/4$ to the end, since the degree of numerator polynomial is larger (in negative) than that of the denominator polynomial, you should first apply a long division of numerator polynomial into denominator polynomial to get :

$$ \frac{ 1 + z^{-1} + z^{-2} + z^{-3} } {1 -z^{-1} }  = -(3 + 2 z^{-1} + z^{-2}) + \frac{ 4} {1 -z^{-1} } $$

hence you have this (including the scale factor):
$$ Y(z) =  -(\frac{3}{4} + \frac{2}{4} z^{-1} +  \frac{1}{4} z^{-2}) + \frac{1} {1 -z^{-1} } $$

And finally you get the output from the PFE method via inverse Z-transformas:
$$ \boxed{ y[n] =  u[n] - \frac{3}{4} \delta[n] - \frac{2}{4} \delta[n-1] -  \frac{1}{4} \delta[n-2] } $$