Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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]
- }
- > 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
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement