CODEWITHSUNDEEP
rvectorsimulation
Alison Meeth
Alison Meeth

Reputation: 95

Vectors for simulations

The code I have so far is written below. This is simulation so there is no actual data. I have two questions:

  1. I have two vectors (treat and cont) but I need to put them into one single vector which I did (vect), however, I need another vector that is coding for treatment vs. control. How do I do that?
  2. For my model (model) I need to fit a linear model testing for a treatment effect but I don't know how to add that effect into what I have or is that what it is testing in the code I have?
    library(car)
    
    treat=rnorm(3, mean = 460, sd = 110)
    treat
    cont=rnorm(3, mean = 415, sd = 110)
    cont
    vect=c(treat, cont)
    vect
    
    nsims = 1000
    
    p.value.saved = coeff.saved = vector()
    
    for (i in 1:nsims) {
      treat=rnorm(3, mean = 460, sd = 110)
      cont=rnorm(3, mean = 415, sd = 110)
      vect=c(treat, cont)
      
      model = glm(treat ~ cont, family = poisson)
      
      p.value.saved[i] = Anova(model)$P[1]
      coeff.saved[i] =  coef(model)
    }

Thank you!

Upvotes: 0

Views: 234

Answers (2)

Robin Gertenbach
Robin Gertenbach

Reputation: 10776

Something like this? (note that you'll get a bunch of warnings for running a poisson regression against continuous data.

n <- 3
nsims <- 10
do.call(
  rbind, 
  lapply(1:nsims, function(.) {
    treat <- rnorm(n, mean = 460, sd = 110)
    cont <- rnorm(n, mean = 415, sd = 110)

    # Instead of vect
    df <- data.frame(
      y = c(treat, cont), 
      x = rep(c("treat", "cont"), each = n)
    )

    # Model the values vs treatment indicator
    model <- glm(y ~ x, data = df, family = poisson)

    # Extract the model's p-value and coefficient of treatment.
    data.frame(p = car::Anova(model)$P, coef = coef(model)[2])
  })
)

Upvotes: 1

khoshnaw
khoshnaw

Reputation: 53

The first creates the string and the second bit will combine them. In your example they are both length 3, hence the 3 repetition in rep("trt",3)

treat_lab = c(rep("control", 3),rep("trt", 3))
treatment <- cbind(treat_lab,c(treat,cont))

Upvotes: 1

Related Questions