How to perform chi-square test for more than 2 variables? I want to test with the chi-squared test if there is a difference between the number of stars of a hotel for each country and the gender of a booker.
An example of my dataset, bookings, is:
country     gender  numbofStars
Germany       male            0
Netherlands female            0
Belgium       male            5
Netherlands   male            3
Germany       male            3
Belgium     female            2
Germany     female            1
Netherlands   male            4

The number of stars can be variable between 0 and 5 stars. Do you have any suggestions to use the chi-squared test on the different groups?
 A: First, your issue is not about multiple levels but more than 2 variables.  It is possible to do chi-square on a 3 dimensional table, but I don't know how to do it in R. However:
Second: Why are you trying to use chi-square here? This seems like a regression problem. Number of stars is the dependent variable, gender and country are independent variables (and you may want their interaction, too).  Since stars can only be a non-negative integer, you may consider a count regression model (e.g. Poisson or negative binomial). 
A: You could compute manually the expectated numbers for country, gender and score, $E_{ijk}$, based on the group totals (with independence as hypothesis): 
$$E_{ijk} = T \frac{T_{i}}{T}*\frac{T_{j}}{T}*\frac{T_{k}}{T}$$
with T the total of individuals, $T_{i}$ the total number of individuals of country $i$, $T_{j}$ the total individuals of gender $j$, $T_{k}$ the total number of individuals of score $k$.
Then together with the observed numbers $O_{ijk}$ and the typical formula $$\chi^2 = \sum (O_{ijk}-E_{ijk})^2/E_{ijk}$$ you get a $\chi^2$ statistic for which you can calculate the cumulative probability pchisq by using $i j k -(i-1)-(j-1)-(k-1)-1 = ijk -i-j-k+2$ degrees of freedom.

Note that the R function chisq.test uses the check for small expectation values (when the normal approximation to the binomial distribution starts to deviate too much)
if (any(E < 5) && is.finite(PARAMETER)) 
    warning("Chi-squared approximation may be incorrect")

and this is not performed when you do these calculations manually.
A: I think I would treat numbofStars as an ordinal variable.  It seems to me to be more akin to a Likert rating than an actual count.
Fortunately, R has the ordinal package, which allows for ordinal regression as flexible as standard OLS regression.
In addition to the analysis of deviance below, Objects from this package are supported by the emmeans package, which you can use for post-hoc testing.
The webpage listed in the example below has a few other tricks, like getting a pseudo r-square for the model.
A couple of notes on the following example:
1) The dependent variable must be an ordered factor in R.
2) The RVAideMemoire package modifies car::Anova to handle the clm object correctly.  Be careful you don't use car::Anova without RVAideMemoire or you may get wildly incorrect results.  If you don't like this, you can use anova to compare models, or emmeans to set up contrasts.
### This example adopted from
### http://rcompanion.org/handbook/G_11.html

if(!require(ordinal)){install.packages("ordinal")}
if(!require(car)){install.packages("car")}
if(!require(RVAideMemoire)){install.packages("RVAideMemoire")}

Data = read.table(header=TRUE, text=
"country     gender  numbofStars
Germany       male            0
Netherlands female            0
Belgium       male            5
Netherlands   male            3
Germany       male            3
Belgium     female            2
Germany     female            1
Netherlands   male            4
")

Data$numbofStars = factor(Data$numbofStars, ordered=TRUE)

levels(Data$numbofStars)

   ### [1] "0" "1" "2" "3" "4" "5"

library(ordinal)

Model = clm(numbofStars ~ gender + country, data=Data)

library(car)

library(RVAideMemoire)

Anova(Model, type = "II")

   ### Analysis of Deviance Table (Type II tests)
   ###
   ### Response: numbofStars
   ###        LR Chisq Df Pr(>Chisq)  
   ### gender    5.5818  1    0.01815 *
   ### country   4.3381  2    0.11429 

A: For this specific case, I like the suggestion in my other answer better.  But this one addresses the question asked:  How to do a test of association on a table of nominal counts with more than two dimensions.
~
There are a couple of relatively easy approaches to determine if there's an association on a contingency table with three dimensions.
The Cochran–Mantel–Haenszel test can be used for 3-dimensions only.  One dimension needs to be identified as the stratifying dimension, so it tests the independence of the variables of the m x n tables within k strata.  Some versions of the test allow for only 2 x 2 tables within the strata.  In your case it makes sense to test for the association of gender and numbofStars, and treat country as the stratifying variable.
Another approach is to use log-linear models, which test for association among nominal variables.  They can test for mutual independence,conditional independence, or partial independence. The references listed here may be helpful in understanding these options.
See the links listed in the code for ideas for post-hoc testing and testing assumptions.
### These examples adopted from
### http://rcompanion.org/handbook/H_06.html
### http://rcompanion.org/handbook/H_08.html

if(!require(MASS)){install.packages("MASS")}

I'll need to simulate some data for these examples.
Note that in constructing the table, country is the stratifying variable, and is listed third in the xtabs function.
set.seed(12345)

gender = c(rep("Male", 36), rep("Female", 36))
country  = rep(c("Belgium", "Germany", "Netherlands"), 24)
numbofStars = c(sample(c(0:4), 36, replace=TRUE), sample(c(1:5), 36, replace=T))

Data = data.frame(country, gender, numbofStars)

Data$numbofStars = factor(Data$numbofStars)

levels(Data$numbofStars)

Table = xtabs( ~ gender + numbofStars + country,
              data=Data)

ftable(Table)

Cochran–Mantel–Haenszel test:
mantelhaen.test(Table)

### Cochran-Mantel-Haenszel test
###
### data:  Table
### Cochran-Mantel-Haenszel M^2 = 19.765, df = 5, p-value = 0.001383

Log-linear model:
library(MASS)

loglm( ~ gender + numbofStars + country,
       Table)

### Call:
### loglm(formula = ~gender + numbofStars + country, data = Table)
###
### Statistics:
###                       X^2 df    P(> X^2)
### Likelihood Ratio 48.69355 27 0.006417156
### Pearson          40.27253 27 0.048305047

