# Simulation with Dichotomous Random Variables in R

#### cw2453

##### New Member
Hi,

I am a bit new to R (I've used Stata before), but I am trying to figure out how to do a simulation in which I want dichotomous random variables whose expectation has a form like that of the distribution function of a logistic random variable.

So far, I've decided to use runif, and I want to use an indicator o show if a uniform is less than the expectation in order to simulate the dichotomous variables.

so far:

set.seed(12345)
runif (1000)

at this point, what should I do? I've generated 1000 "pseudo-random" numbers.. how would I go about with the second step?

#### Lazar

##### Phineas Packard
Could you not just simply use a single step like:
Code:
sample(0:1, 100, replace=TRUE, prob=c(.65,.35))

#### Dason

but I am trying to figure out how to do a simulation in which I want dichotomous random variables whose expectation has a form like that of the distribution function of a logistic random variable.
Can you explain what you mean by this? Do you mean you want to simulate data that is appropriate for a logistic regression? A logistic random variable is not the same thing. If you meant what I asked originally then do you have certain values for the covariate(s) that you're interested in?

#### Lazar

##### Phineas Packard
Is what you are suggest dason something like:

Code:
library(arm)
n <- 1000
a <- .25
b <- .56
x <- rnorm(n)
logprobs <- c(invlogit(a + b*x))
#So that the data frame
myData <- data.frame(x=x, y=sapply(1:n, function (x)sample(0:1, 1, replace=TRUE, prob=c(1-logprobs[x], logprobs[x])) ) )
#can be used in something like
summary(glm(y ~ x, family=binomial("logit"), data=myData))

#### Dason

Is what you are suggest dason something like:

Code:
library(arm)
n <- 100
a <- .25
b <- .56
x <- rnorm(n)
logprobs <- c(invlogit(a + b*x))
sapply(1:100, function (x)sample(0:1, 1, replace=TRUE, prob=c(1-logprobs[x], logprobs[x])) )
arm is quite a dependency just to use invlogit. Especially since
Code:
invlogit <- function (x){
1/(1 + exp(-x))
}
I guess I would write it like this:
Code:
n <- 100
x <- seq(-3, 3, length.out = n)
b0 <- .25
b1 <- .56
probs <- 1/(1 + exp(-(b0 + b1*x)))
out <- rbinom(n, 1, probs)
dat <- data.frame(x = x, outcome = out)