How to simulate Likert-scale data in R? I want to do a Monte-Carlo experiment to check my theoretical finding regarding Likert-scale data.
I want to have random test results with N participants, where every question has 5 answers: (0,1,2,3 and 4), given probabilities (0.1, 0.2, 0.4, 0.2, 0.1)
 A: A likert scale, as the term is typically used, is just an ordinal rating scale.  The phrase is often used for a single rating, which might have been called a likert item.  Traditionally, the idea was that you would have a set of likert items that all measure the same thing and have the same measurement properties.  The result is that you could sum (or average) the items and end up with a good measure of something that approximated a continuous, interval scale (see: levels of measurement).
On the other hand, to simulate data, you need to know the distribution that the data should have.  More generally, for simulation studies people generally want to have a data generating process for the resulting distribution.  A likert scale is a type of data gathering instrument, not a distribution and not a data generating process.
Thus, what you ultimately need is to specify a data generating process that you believe is appropriate for the eventual likert data that you want to simulate.  After that, there is just the trivial implementational details specific to the software you intend to use (in your case, R).  Because people conceptualize likert data as manifest data derived from a latent variable, the most common approach would be to simulate the latent variable according to the theorized distribution (perhaps a normal distribution), and then have a function that maps it to a small ordered set of numbers (e.g., $1, \ldots, 5$).  Note that moving from the latent to the manifest variable makes many of the parameters of the latent variable's distribution unidentifiable, so you often needn't bother worrying about them.  A simple approach would be to have just the two steps move directly to the final rating, but a more comprehensive approach could model each item with their own set of the two steps, and then have the likert scale combined from the items just as they would be in a real case.

Here is an example, coded in R.  I will imagine that there are 5 items that measure the same construct.  As such, they are moderately correlated.  Two items might be 'reverse scored', but I will assume that doesn't affect the result appreciably so I won't simulate that.  However, I will make some more strongly related to the underlying variable than others, and I will make some biased towards higher or lower ratings.
set.seed(8649)     # this makes the example exactly reproducible
N      = 10        # this is how much data I'll generate
latent = rnorm(N)  # this is the actual latent variable I want to be measureing

##### generate latent responses to items
item1 = latent + rnorm(N, mean=0, sd=0.2)  # the strongest correlate
item2 = latent + rnorm(N, mean=0, sd=0.3)
item3 = latent + rnorm(N, mean=0, sd=0.5)
item4 = latent + rnorm(N, mean=0, sd=1.0)
item5 = latent + rnorm(N, mean=0, sd=1.2)  # the weakest

##### convert latent responses to ordered categories
item1 = findInterval(item1, vec=c(-Inf,-2.5,-1, 1,2.5,Inf))  # fairly unbiased
item2 = findInterval(item2, vec=c(-Inf,-2.5,-1, 1,2.5,Inf))
item3 = findInterval(item3, vec=c(-Inf,-3,  -2, 2,3,  Inf))  # middle values typical
item4 = findInterval(item4, vec=c(-Inf,-3,  -2, 2,3,  Inf))
item5 = findInterval(item5, vec=c(-Inf,-3.5,-3,-1,0.5,Inf))  # high ratings typical

##### combined into final scale
manifest = round(rowMeans(cbind(item1, item2, item3, item4, item5)), 1)
manifest
# [1]  3.4  3.6  3.4  3.8  2.6  3.4  3.2  2.0  3.8  3.2
round(latent, 1)
# [1]  1.3  0.6  0.2  1.0 -1.5  0.1  0.4 -2.5  2.3 -0.3
cor(manifest, latent)
# [1] 0.9280074

A: To perform the simulation, here is a one line solution using the sample function:  
sample(0:4, N, replace = TRUE, prob = c(0.1, 0.2, 0.4, 0.2, 0.1))

#where:
# 0:4 is the sequence of values (0 to 4 in this case)
# N is the number of samples (participants)
# replace = TRUE for sampling with replacement 
# prob = c(0.1, 0.2, 0.4, 0.2, 0.1) is the probability of selection for each score.

A: One way to generate Likert data is according to a proportional odds model. Here, the underlying distribution of the (latent) response is a logistic random variable with a center $\mu$ that can vary as a function of one or more predictors. The latent variable is then thresholded into as many categories with an arbitrary number of cutpoints. Achieving a target number of response categories is quite difficult, either requiring advanced math or ( more likely) play it by ear.
set.seed(123)
n <- 1e6
beta <- 0.3
alpha <- sort(rnorm(5))
x <- seq(-3, 3, length.out = n)
z <- rlogis(n, beta*x)
y <- factor(findInterval(z, alpha))
library(MASS)
fit <- polr(formula = y ~ x)

Generates the association:
> coef(fit)
        x 
0.2982142 

a log-odds ratio for endorsing any higher response comparing two groups differing by 1 unit of $X$.
A: Original source is here: http://ravshansk.com/articles/likert.html
The following formula works not only for Likert-scale models, but for any categorically distributed variables.
Suppose, you want to generate a 5-category data (x1, x2, x3, x4, x5) for N participants with probabilities (1/10, 2/10, 4/10, 2/10, 1/10). The following formula will work:
distribution <- c(rep(x1,1),rep(x2,2),rep(x3,4),rep(x4,2),rep(x5,1))
potential_population <- rep(distribution, M) #M is any number greatN
likert_data <- sample(potential_population, N)

The main idea here is to write a list "distribution" with appropriate number of repetitions that would together satisfy the desired probabilities. Needless to say that you have to set common denominators and integer numerators for the probabilities.
