This is a variant of finding a total-degree based power series.

linearize[poly_] := Module[
  {vars = Variables[poly], t},
  Normal[Series[poly /. Thread[vars -> t*vars], {t, 0, 1}]] /. t -> 1
  ]

The example:

poly = 1 + x + y + z + x y + x y z;
linearize[poly]

(* Out[50]= 1 + x + y + z *)

I will give a suggestion for a procedure which (attempts) to form a Taylor series for a function of several variables in another thread.

Seems to be a bit tricky. Here a suggestion which , I think, has to be tested with other examples.

Say your function is

f = gamma + 3 x + a y + z + x y + x y z;

and you know the variables are

var = {x, y, z};

Get the constant term, if any, by setting all variables to 0

In[3]:= const = f /. (# -> 0 & /@ var)    
Out[3]= gamma

Make a new function f1 ( f without constant term )

In[4]:= f1 = f - const    
Out[4]= 3 x + a y + x y + z + x y z

Now for each variable set the other ones to zero and expand the (new) function in a series up to order 1

h[x_] := Normal[Series[f1 /. ( #1 -> 0 &) /@ Complement[var, {x}], {x, 0, 1}]]

Do this for every variable

In[8]:= h /@ var
Out[8]= {3 x, a y, z}

and finally get your linearized function f1L

In[9]:= f1L = const + Plus @@ (h /@ var)
Out[9]= gamma + 3 x + a y + z

Certainely you can combine all these steps in a single function

linF[f_, var_] := Module[{},
  con = f /. (# -> 0 & /@ var);
  f1 = f - con;
  con + Plus @@ (Function[x, Normal[Series[f1, {x, 0, 1}] /. (a_. #1 -> 0 &) /@ 
          Complement[var, {x}]]] /@ var)]

Then for example

In[17]:= linF[123 a + 4 x + beta y + 45 x z - 78 x y z^2, {z, y, x}]
Out[17]= 123 a + 4 x + beta y