The standard statistical tools are the correlation coefficient (see Michael Chernick's answer) which is a value between [-1,1] and is unit independent. Related to the correlation coefficient is the covariance. The covariance is affected by units but may be easier to interpret. However, I don't like either of those options in the general case. I don't like them because they are not conformal transformation independent. Consider that a straight horizontal or vertical line is deemed to be non-linear by both of those measures.
A better unitless option is to use a singular value decomposition (SVD). The SVD breaks data up into component pieces ranked by magnitude of their contribution to the whole. The ratio of the largest singular number to the second largest singular number is hence a metric of linearity. Note, that to use this method you must first centralize the data (make average X, Y, Z, etc coordinates equal zero).
SVD, D matrix:
Ratio of singular Values
The ratio above can be roughly interpreted as the data being three times as long in the direction of the best fit line as it is cross-linear.
For a solution with units that has units and doesn't require an SVD.
Do some line fit that has the center of the line as one of the parameters. Using the centralized data above this is simple:
line pt = 0 0 (always the case for centralized data)
line direction = -0.999956849 -0.009289783
Vectors from the center of the line to each point are the centralized coordinates of the points. Find the lengths of the of the projection of these vectors onto the line (absolute value of the vector dot the line direction), and the length of the perpendicular vector component (length of vector cross line direction).
Length Parallel, Length Perpendicular
The maximum of parallel projections is the stretch of data along the line. The maximum length of the perpendicular projection is a measure of the non-linearity. The ratio of the two is an approximation of the singular values ratio above.
1. Affine invariance in linearity is not possible. Consider, in an affine transformation we could scale all but one of the coordinate axes to near zero (making any set of points linear). So conformal invariance is the best we can do.
2. These methods are NOT ROBUST to outlier data.
3. Examples are 2D but generalized to N-dimensional.