For an arbitrary $n\times n$ matrix $B$, define $$GL_B(n, \Bbb{R}) = \{A\in GL(n,\Bbb{R}) \,\vert\, A^TBA = B\}$$ Further, the Lie algebra of a Lie group, here, is defined as the tangent space at the identity.
The problem here is to show that $GL_B(n, \Bbb{R})$ is a Lie subgroup of $GL(n, \Bbb{R})$, and to find its Lie algebra.
For the first part, I tried to base it off a proof that $O_n$ is a Lie subgroup of the general linear group. It's obvious that $GL_B(n, \Bbb{R})$ is a subgroup of $GL(n, \Bbb{R})$. To show it's a submanifold, it suffices to find a smooth function $F$ so that $GL_B(n, \Bbb{R})$ is the inverse image of a regular value of $F$. Define $F: GL(n, \Bbb{R})\to M_B$, where $M_B$ is an appropriate vector subspace of $M^{n\times n}(\Bbb{R})$, such that $F(A) = A^TBA$. For example, if $B = I_n$, $M_n$ would be the set of symmetric matrices.
$GL(n, \Bbb{R})$ is identified with $\Bbb{R}^{n^2}$ in the natural manner, and so, given $X$, $Y$, we have
$$dF_X(Y) = \lim_{t\to 0} \frac{F(X+tY)-F(X)}{t} = X^TBY + Y^TBX$$ We wish to show that $B$ is a regular value of $F$; thus, assume $F(X) = B$, i.e. $X^TBX = B$. Now we need to show that the function $Y\mapsto X^TBY + Y^TBX$ is surjective. Thus, pick some $Q$ in $M_B$. At this point, I am unsure how to construct a $P$ so that $dF_X(P) = Q$. If $B=I_n$, I could choose $P = \frac{1}{2}XQ$. However, I am unsure how to construct such a $P$ (as well as what $M_B$ should be) for a general $B$.
As for the Lie algebra, I'm fairly confident it's the set of $n\times n$ matrices $X$ so that $BX = -X^TB$; however, I just got this from messing around, and I'm not clear as to how I can prove this is the Lie algebra.
To summarize, I have three questions:
What should my choice of $M_B$ be?
How can I show explicitly that $dF_X$ is surjective for $X\in F^{-1}(B)$?
How do you prove that the Lie algebra of $GL_B(n, \Bbb{R})$ is $\{X\in M^{n\times n}(\Bbb{R})\,\vert\,BX=-X^TB\}$?
Thank you for your time. If my approach is hopeless, then how would one go about this problem?
The efficient way to do this certainly is via the closed subgroup theorem. This says that a subgroup of $GL(n,\mathbb R)$ which at the same time is a closed subset for the topology induced from $\mathbb R^{n^2}$ automatically is a Lie subgroup. Now $\{X\in\mathbb R^{n^2}:X^TBX=B\}$ is a closed subset of $\mathbb R^{n^2}$ by continuity of matrix multiplication, so the theorem applies. This theorem also gives you two descriptions of the Lie algebra. The first is as the set of derivatives at $0$ of smooth curves from an open intervall containing $0$ to $\mathbb R^{n^2}$ which map $0$ to the identity matrix and have values in $GL_B(n,\mathbb R)$. This description readily shows that the Lie algebra of $GL_B(n,\mathbb R)$ has to be contained in $\{X:X^TB=-BX\}$. On the other hand, you can describe the Lie algebra as the set of those matrices $X$ such that $\exp(tX)\in GL_B(n,\mathbb R)$ for all $t\in\mathbb R$. Now using the good convergence property of exponential series and the binomial theorem, one easily verifies that $X^TB=-BX$ implies $\exp(tX)\in GL_B(n,\mathbb R)$ for all $t\in\mathbb R$, so you conclude that the Lie algebra of $GL_B(n,\mathbb R)$ actually equals $\{X:X^TB=-BX\}$.
If you are looking for an elementary proof along the lines you describe, things get a bit more complicated. One thing to notice is that the argument used in the case of $O(n)$ is a bit of an overkill for the problem. Having a regular value is actually more than you need in order that the pre-image is a smooth submanifold. By the constant rank theorem, it actually suffices to show (in your notation) that the rank of the linear map $Y\mapsto A^TBY+Y^TBA$ is the same for all $A\in GL_B(n,\mathbb R)$. If you call that rank $k$, then what you are trying to show is that there is a fixed $k$-dimensional subspace to which this map is a linear isomorphism for any $A\in GL_B(n,\mathbb R)$, which of course is much stronger (and need not be true in general).
A final possibility for an elementary approach would be to bring $B$ to some normal form. This also should be possible but slightly complicated, since you can decompose $B$ into a symmetric and a skew symmetric part, which both are preserved by each $A\in GL_B(n,\mathbb R)$, so there are quite a few possible normal forms that have to be taken into account. This also shows that the dimension of $GL_B(n,\mathbb R)$ hevaily depends on $B$, so the groups look much less uniform than the uniform simple description of their Lie albegras suggests.
Edits: For general $B$, I am pretty sure that your approach is hopeless (at least in the current form). As mentioned above, in order to get a submanifold, it suffices that, viewed as a map to the space of $n\times n$-matrices, the rank of the derivative $dF_X$ is the same for all $X\in GL_B(n,\mathbb R)$. I am pretty sure that it is not too difficult to verify this by elementary methods. But this of course is much waeker than the fact that $F$ itself has values in a linear subpsace and each $X\in GL_B(n,\mathbb R)$ is a regular point for $F$ as a function to that subspace. (Using the inverse function theorem, one proves that, given the condition on constant rank, one may locally change $F$ using a diffeomorphism in such a way that the stronger condition is satisfied, but there is no reason why it should work for the nice function $F$ you start with.)
The approach you use should work if $B$ is invertible and either symmetric or skew symmetric (which only works in even dimensions). In this case, you can take $M_B=\{X:(XB)^T=BX\}$, and the dimension of $G_B$ only depends on whether $B$ is symmetric or skew symmetric. (Up to isomorphism, you get the pseudo-orthogonal and symplectic groups, so indeed this are not too many different cases.)
In more general cases, I don't think that $F$ has values in a linear subspace of small enough dimension, since the maps $A^TBA$ have different kernels and images. Things remain reasonably easy if $B$ is either symmetric or skew symmetric, in which case the possible matrices are classified by rank and (only in the symmetric case) signature. In these cases you can determine the dimension in terms of the rank using linear algebra.
If $B$ is neither symmetric nor skew symmetric, things get complicated. Writing $B=B_s+B_a$ with $B_s=\frac12(B+B^T)$ and $B_a=\frac12(B-B^T)$ you see that any $A\in GL_B(n,\mathbb R)$ preserves both $B_s$ and $B_a$. In addition to the properties of these two forms, you now have an issue of "relative position". For example, if $B_s$ is positive definite it defines an inner product and one can "diagonalize" $B_a$ with respect to this inner product, leading to eigenvalues of the form $\pm s_i$ with $s_i\in\mathbb R$. Now the size of the stabilizer depends on the multiplicities among these eigenvalues. The case of general $B_s$ is even more complicated ...