How to best analyze a dataset with only categorical variables?

I have a dataset with all multinomial (categorial) variables. The column names look like:

religionID, occupationID, sexID, ..., changedID

The independent variables are: "religionID" (a factor with 7 levels), "occupationID" (a factor with 3 levels), .... There are 15 total, all factors with multiple levels.

The dependent variable is: "changedID" (a factor with 4 levels).

In all variables, the ordering of the levels is arbitrary. (E.g., factor level 1 is not "better than" or "greater than" factor level 2)

I have about 500 rows of such data.

My goals:

1. Find out which independent variables effect the dependent variables (e.g., does religion matter at all to changedID?)

2. Find out which independent variables have more effect than other independent variables (e.g., is religion more important than sex?)

3. Find out which factor levels in each independent variable has the most effect (e.g., which religionID factor leads the most often to changedIDs of factor level "2"?)

Of course, I would also like to know if the results above are statistically significant.

I've thought of building contingency tables (using R's table() function), and then running chisq.test() on the result. I guess this could be used to satisfy goal #1, if I build one such table for each independent variable. But it doesn't help with goal #2 and goal #3.

Multinomial regression might be the answer, but I am confused as to whether it applies to my situation, and if it does, how to use it and interpret the results.

Bonus points if I can have an solution in R.

-