3
$\begingroup$

I am trying to generate a data frame of fake data for exploratory purposes. Specifically, I am trying to produce data with a binary dependent variable (say, failure/success), and a categorical independent variable called 'picture' with 5 levels (pict1, pict2, etc.). I am following the answer provided here, which allows me to successfully generate the data. However, I need each level of 'picture' to occur the same number of times (i.e. 11 repetitions of each level = 55 total observations per subject).

Here is a reproducible example of what has worked so far (code from user: ocram):

library(dummies)

#------ parameters ------
n <- 1000 
beta0 <- 0.07
betaB <- 0.1
betaC <- -0.15
betaD <- -0.03
betaE <- 0.9
#------------------------

#------ initialisation ------
beta0Hat <- rep(NA, 1000)
betaBHat <- rep(NA, 1000)
betaCHat <- rep(NA, 1000)
betaDHat <- rep(NA, 1000)
betaEHat <- rep(NA, 1000)
#----------------------------

#------ simulations ------
for(i in 1:1000)
{
  #data generation
  x <- sample(x=c("pict1","pict2", "pict3", "pict4", "pict5"), 
              size=n, replace=TRUE, prob=rep(1/5, 5))  #(a)
  linpred <- cbind(1, dummy(x)[, -1]) %*% c(beta0, betaB, betaC, betaD, betaE)  #(b)
  pi <- exp(linpred) / (1 + exp(linpred))  #(c)
  y <- rbinom(n=n, size=1, prob=pi)  #(d)
  data <- data.frame(picture=x, choice=y)

  #fit the logistic model
  mod <- glm(choice ~ picture, family="binomial", data=data)

  #save the estimates
  beta0Hat[i] <- mod$coef[1]
  betaBHat[i] <- mod$coef[2]
  betaCHat[i] <- mod$coef[3]
  betaDHat[i] <- mod$coef[4]
  betaEHat[i] <- mod$coef[5]
}

However, as you can see from the output, each level of the factor 'picture' does not occur the same number of times (i.e. 200 times each).

> summary(data)
picture     choice     
pict1:200   Min.   :0.000  
pict2:207   1st Qu.:0.000  
pict3:217   Median :1.000  
pict4:163   Mean   :0.559  
pict5:213   3rd Qu.:1.000  
            Max.   :1.000 

Moreover, it is not entirely clear to me how to manipulate the initial beta values as to determine the probability of success/failure for each level of 'picture'. I cannot comment the original question because I do not yet have the necessary reputation points.

$\endgroup$
1
  • $\begingroup$ I'm still working on this, just in case anybody has any more suggestions. $\endgroup$ – babylinguist Sep 25 '14 at 15:28
3
$\begingroup$
  1. If you want 200 copies of each of 5 levels in a polytomous variable in random order, then do this instead:

    x <- sample(rep(paste0('pict', 1:5), 200))
    
  2. If you want to control for overall prevalence of a specific outcome, then you must choose which beta you will fudge. I usually do beta0.

    MM         <- model.matrix(~x)
    betas      <- rnorm(4)
    prevTarget <- 0.3
    prevDiff   <- function(beta0)  prevTarget - 
                                   mean(binomial()$linkinv(MM%*%c(beta0, betas)))
    beta0      <- uniroot(prevDiff, c(-100, 100))$root
    mean(binomial()$linkinv(MM%*%c(beta0, betas)))
    
$\endgroup$
8
  • $\begingroup$ Could you give a little more detail as to how I can use your answer in conjunction with the code I provided above? If it is not to be used with the code above, how exactly can I generate a data frame from your code? $\endgroup$ – babylinguist Sep 22 '14 at 18:38
  • $\begingroup$ @jvcasill not following your question. There's a lot of replication, e.g. using model.matrix instead of dummy and binomial()$linkinv instead of a by-hand expit. The argument to mean in the last line is the fitted values, whence you can call rbinom. May I suggest you run the two coding examples line-by-line and figure out what you're trying to do? $\endgroup$ – AdamO Sep 22 '14 at 18:48
  • $\begingroup$ I am trying to generate a data frame similar to the one the original code generates. However I want an equal amount of observations for each level of the factor 'picture'. I also want to be able to manipulate the probability of success/failure at each level. It is clear to me that your code does something, but I don't understand how to use it for what I am asking. Does that make sense? $\endgroup$ – babylinguist Sep 22 '14 at 18:49
  • $\begingroup$ Append y = rbinom(1000,1,binomial()$linkinv(MM%*%c(beta0, betas))), and data <- data.frame(picture=x, choice=y). $\endgroup$ – AdamO Sep 22 '14 at 19:22
  • 1
    $\begingroup$ @jvcasill Please look here. github.com/aomidpanah/simulations/blob/master/logregsim.r Given the work you are trying to do, I think it would be good for you to explore this in some depth and better understand how these functions are working. As I suggested, line-by-line analysis and scansion will be a very good way to do this. Also examine the output from proc.time. $\endgroup$ – AdamO Sep 22 '14 at 19:58
0
$\begingroup$

I eventually found an acceptable answer to my question. This may or may not be the best/most sophisticated way to handle this, but it enabled me to get what I was after.

set.seed(1)

intercept1 = -6.0
beta1      = 2.75
xtest1     = rnorm(1300, 1, 3)
linpred1   = intercept1 + xtest1 * beta1
prob1      = exp(linpred1)/(1 + exp(linpred1))
runis1     = runif(1300, 0, 1)
ytest1     = ifelse(runis1 < prob1, 1, 0)

intercept2 = -7.0
beta2      = 2.75
xtest2     = rnorm(1300, 1, 3)
linpred2   = intercept2 + xtest2 * beta2
prob2      = exp(linpred2)/(1 + exp(linpred2))
runis2     = runif(1300, 0, 1)
ytest2     = ifelse(runis2 < prob2, 1, 0)

intercept3 = -7.0
beta3      = 2.75
xtest3     = rnorm(1300, 1, 3)
linpred3   = intercept3 + xtest3 * beta3
prob3      = exp(linpred3)/(1 + exp(linpred3))
runis3     = runif(1300, 0, 1)
ytest3     = ifelse(runis3 < prob3, 1, 0)


newdf1 <- data.frame(prop = ytest1, x = xtest1, week = "w1")
newdf1 <- newdf1[with(newdf1, order(week, x)), ]
newdf1$stim <- rep(seq(-60, 60, by = 10), each = 100)

newdf2 <- data.frame(prop = ytest2, x = xtest2, week = "w3")
newdf2 <- newdf2[with(newdf2, order(week, x)), ]
newdf2$stim <- rep(seq(-60, 60, by = 10), each = 100)

newdf3 <- data.frame(prop = ytest3, x = xtest3, week = "w6")
newdf3 <- newdf3[with(newdf3, order(week, x)), ]
newdf3$stim <- rep(seq(-60, 60, by = 10), each = 100)

temp <- rbind(newdf1, newdf2)
df <- rbind(temp, newdf3)

I am now able to manipulate the data frame, perform logistic regression, and produce plots. Thank you to everybody who provided input.

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ You don't need the runif() & ifelse() calls, you can just use rbinom(). You may also want to define a new function lo2p = function(lo){ exp(lo)/(1+exp(lo)) } for clearer code & to save you some typing. $\endgroup$ – gung - Reinstate Monica Nov 30 '14 at 3:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.