When rendering a cube, everything looks proportional straight on, however, when I start to move my camera up so that it's looking down at my cube, everything starts to squish. I've attached how the cube looks straight on (correct) and how it looks when the camera's above it looking down (squished).

Also, here's my projection matrix code:

Matrix4 Matrix4::Perspective(float fov, float aspectRatio, float zNear, float zFar)
{
	Matrix4 resultMat;
	float zDelta = zFar - zNear;
	float halfRadians = Math::ToRadians(fov / 2.0f);
	float halfTanRad = tanf(halfRadians);
	
	// setup the projection matrix
	resultMat.m[0] = 1.0f / (halfTanRad * aspectRatio);
	resultMat.m[5] = 1.0f / halfTanRad;
	resultMat.m[10] = -(zFar + zNear) / zDelta;
	resultMat.m[11] = -1.0f;
	resultMat.m[14] = -2.0f * zNear * zFar / zDelta;
	resultMat.m[15] = 0.0f;
	return resultMat;
}

Also, my aspect ratio is always set correctly whenever the window resizes (window.x / window.y), so there shouldn't be any issues there.

Your projection matrix doesn't change depending on the view-angle. It stands to reason that the projection matrix is an unlikely culprit.

Check your camera view matrix or cube transform. Looks like the pre-projection coordinate system axis aren't perpendicular, or have different lengths.

I was also leaning towards the view matrix myself after checking the projection. Here's how I set my view matrix for this example:

viewMat = Matrix4::LookAt(Vector3(cosf(panOffset.x) * 3.0f, panOffset.y, -sinf(panOffset.x) * 3.0f), Vector3::Zero());
viewMat.Inverse();

Matrix4::LookAt() is your typical LookAt method which takes the camera's position, target position and up vector (set to <0.0, 1.0, 0.0>, if non specified which is why it's omitted). Here's the source for that:

Matrix4 Matrix4::LookAt(Vector3 position, Vector3 target, Vector3 up)
{
	Matrix4 resultMat;
	Vector3 forward = target - position;
	forward.Normalize();
	Vector3 right = Vector3::Cross(forward, up);
	right.Normalize();
	
	resultMat.Set(    right.x,	up.x,	-forward.x,	position.x, 
			  right.y,	up.y,	-forward.y,	position.y,
			  right.z,	up.z,	-forward.z,	position.z,
			  0.0f,		0.0f,	0.0f,		1.0f);
	return resultMat;
}

So, if this looks correct, then it's most-likely the Inverse() method for Matrix4:

bool Matrix4::Inverse()
{
	float invDet = Determinant();
	float values[16] = {
		0.0f, 0.0f, 0.0f, 0.0f,
		0.0f, 0.0f, 0.0f, 0.0f,
		0.0f, 0.0f, 0.0f, 0.0f,
		0.0f, 0.0f, 0.0f, 0.0f
	};
	
	// get the determinant
	if(invDet == 0.0f)
	{
		printf("det = 0.0f\n");
		return false;
	}
	invDet = 1.0f / invDet;
	
	// process the first column
	values[ 0] = m[ 5] * m[10] * m[15] - m[ 5] * m[11] * m[14] - m[ 9] * m[ 6] * m[15] +
				 m[ 9] * m[ 7] * m[14] + m[13] * m[ 6] * m[11] - m[13] * m[ 7] * m[10];
	values[ 1] =-m[ 1] * m[10] * m[15] + m[ 1] * m[11] * m[14] + m[ 9] * m[ 2] * m[15]
				-m[ 9] * m[ 3] * m[14] - m[13] * m[ 2] * m[11] + m[13] * m[ 3] * m[10];
	values[ 2] = m[ 1] * m[ 6] * m[15] - m[ 1] * m[ 7] * m[14] - m[ 5] * m[ 2] * m[15] +
				 m[ 5] * m[ 3] * m[14] + m[13] * m[ 2] * m[ 7] - m[13] * m[ 3] * m[ 6];
	values[ 3] =-m[ 1] * m[ 6] * m[11] + m[ 1] * m[ 7] * m[10] + m[ 5] * m[ 2] * m[11] -
				 m[ 5] * m[ 3] * m[10] - m[ 9] * m[ 2] * m[ 7] + m[ 9] * m[ 3] * m[ 6];
	
	// process the second column
	values[ 4] =-m[ 4] * m[10] * m[15] + m[ 4] * m[11] * m[14] + m[ 8] * m[ 6] * m[15] -
				 m[ 8] * m[ 7] * m[14] - m[12] * m[ 6] * m[11] + m[12] * m[ 7] * m[10];
	values[ 5] = m[ 0] * m[10] * m[15] - m[ 0] * m[11] * m[14] - m[ 8] * m[ 2] * m[15] +
				 m[ 8] * m[ 3] * m[14] + m[12] * m[ 2] * m[11] - m[12] * m[ 3] * m[10];
	values[ 6] =-m[ 0] * m[ 6] * m[15] + m[ 0] * m[ 7] * m[14] + m[ 4] * m[ 2] * m[15] -
				 m[ 4] * m[ 3] * m[14] - m[12] * m[ 2] * m[ 7] + m[12] * m[ 3] * m[ 6];
	values[ 7] = m[ 0] * m[ 6] * m[11] - m[ 0] * m[ 7] * m[10] - m[ 4] * m[ 2] * m[11] +
				 m[ 4] * m[ 3] * m[10] + m[ 8] * m[ 2] * m[ 7] - m[ 8] * m[ 3] * m[ 6];
	
	// process the third column
	values[ 8] = m[ 4] * m[ 9] * m[15] - m[ 4] * m[11] * m[13] - m[ 8] * m[ 5] * m[15] +
				 m[ 8] * m[ 7] * m[13] + m[12] * m[ 5] * m[11] - m[12] * m[ 7] * m[ 9];
	values[ 9] =-m[ 0] * m[ 9] * m[15] + m[ 0] * m[11] * m[13] + m[ 8] * m[ 1] * m[15] -
				 m[ 8] * m[ 3] * m[13] - m[12] * m[ 1] * m[11] + m[12] * m[ 3] * m[ 9];
	values[10] = m[ 0] * m[ 5] * m[15] - m[ 0] * m[ 7] * m[13] - m[ 4] * m[ 1] * m[15] +
				 m[ 4] * m[ 3] * m[13] + m[12] * m[ 1] * m[ 7] - m[12] * m[ 3] * m[ 5];
	values[11] =-m[ 0] * m[ 5] * m[11] + m[ 0] * m[ 7] * m[ 9] + m[ 4] * m[ 1] * m[11] -
				 m[ 4] * m[ 3] * m[ 9] - m[ 8] * m[ 1] * m[ 7] + m[ 8] * m[ 3] * m[ 5];
	
	// process the fourth column
	values[12] =-m[ 4] * m[ 9] * m[14] + m[ 4] * m[10] * m[13] + m[ 8] * m[ 5] * m[14] -
				 m[ 8] * m[ 6] * m[13] - m[12] * m[ 5] * m[10] + m[12] * m[ 6] * m[ 9];
	values[13] = m[ 0] * m[ 9] * m[14] - m[ 0] * m[10] * m[13] - m[ 8] * m[ 1] * m[14] +
				 m[ 8] * m[ 2] * m[13] + m[12] * m[ 1] * m[10] - m[12] * m[ 2] * m[ 9];
	values[14] =-m[ 0] * m[ 5] * m[14] + m[ 0] * m[ 6] * m[13] + m[ 4] * m[ 1] * m[14] -
				 m[ 4] * m[ 2] * m[13] - m[12] * m[ 1] * m[ 6] + m[12] * m[ 2] * m[ 5];
	values[15] = m[ 0] * m[ 5] * m[10] - m[ 0] * m[ 6] * m[ 9] - m[ 4] * m[ 1] * m[10] +
				 m[ 4] * m[ 2] * m[ 9] + m[ 8] * m[ 1] * m[ 6] - m[ 8] * m[ 2] * m[ 5];
	
	// apply computed values to the matrix
	for(int i=0;i<16;++i)
		m[i] = values[i] * invDet;
	return true;
}

My matrix is laid out in memory column-major, as you'd see in most tutorials (I calculate my matrices as P*V*M rather than M*V*P due to this):

// col1  col2  col3  col4
   m[ 0] m[ 4] m[ 8] m[12] // row 1
   m[ 1] m[ 5] m[ 9] m[13] // row 2
   m[ 2] m[ 6] m[10] m[14] // row 3
   m[ 3] m[ 7] m[11] m[15] // row 4

Now, a lot of this math is over my head. I think my inverse method is based off of proofs I found online on a math website.

LookAt looks incorrect as forward and up will not be perpendicular, so after you have calculated right you must re-calculate up from cross(right, forward). (Or cross(forward, right) perhaps..)

You were completely right! I took your advise, and also checked the source to gluLookAt(). Looks like I forgot a cross product after I got the right and forward vectors to recompute the up vector. The code's exactly the same, except I added that cross product after right.normalize() like so:

Matrix4 Matrix4::LookAt(Vector3 position, Vector3 target, Vector3 up)
{
	Matrix4 resultMat;
	Vector3 forward = target - position;
	forward.Normalize();
	Vector3 right = Vector3::Cross(forward, up);
	right.Normalize();
	up = Vector3::Cross(right, forward);
	
	resultMat.Set(right.x,	up.x,	-forward.x,	position.x,
				  right.y,	up.y,	-forward.y,	position.y,
				  right.z,	up.z,	-forward.z,	position.z,
				  0.0f,		0.0f,	0.0f,		1.0f);
	return resultMat;
} 

It appears to work now.