1
$\begingroup$

I've created an example table (just in order to create a function) with:

ex<-data.frame(b=c(rep('A',50),rep('B',30), rep('C',20)), 
fl=round(runif(100,0,1),0),r=runif(100,0,0.5))
ex2<-cbind(ex,model.matrix(~b-1,ex))
lineal<-ex2$bB+ex2$bA*ex$fl+ex$fl
ex$clase<-round(1/(1+exp(-lineal)),0)

Then I run a logistic regression model (MASS library)

fm<-as.formula(clase~b+fl+r)
modT<-glm(clase~1, family=binomial, data = ex)
modT<-stepAIC(modT, scope = fm, family=binomial, data =ex, k = 4)
summary(modT)

As you can see coefficients are not significant, but I've created the class using them. So I don't understand why this is happening.

enter image description here

$\endgroup$
2
  • $\begingroup$ First, you should actually draw the outcome rather than round it up (not least to remove the numerical errors), e.g. as ex$clase <- rbinom(nrow(ex), prob=1/(1+exp(-lineal)), size=1). $\endgroup$ Apr 12, 2016 at 2:59
  • $\begingroup$ Second, you would not be surprised to fail to recover the original model and coefficients if you had defined your data as ex<-data.frame(b=c(rep('A',5),rep('B',3), rep('C',2)), fl=round(runif(10,0,1),0),r=runif(10,0,0.5)). Maybe ask yourself why that is. $\endgroup$ Apr 12, 2016 at 3:01

1 Answer 1

3
$\begingroup$

You have three problems, you are not randomizing when creating the dependent variable, and you are using the r variable in the GLM stepwise algorithm (and not when creating the dependent variable). In order to randomize, you need to create a variable that will be 1,0 with probability = ex$clase (according to your terminology). You are also running a model with intercept, and there is none defined in the simulation.

If you run this, you will see that the pvalues are almost 0. Some parameter combinations will yield a flatter likelihood, which will end up showing in greater pvalues ( with these ones you will get almost 0 for all of them) - since the stderrors will correspond to the inverse of the information matrix (the curvature of the likelihood)

ex <- data.frame(b=c(rep('A',50),rep('B',30), rep('C',20)), 
fl = round(runif(1000,0,1),0),r=runif(1000,0,1))
ex2 <- cbind(ex,model.matrix(~b-1,ex))
ex2$prob <- 1/(1+exp(-2*ex2$fl + ex2$bA + ex2$bB + ex2$bC))

ex2$clase <- ifelse((ex2$r>ex2$prob),1,0)

fm <- as.formula(clase~fl-1 + bA + bB + bC )
modT <- glm(clase~fl-1 + bA + bB + bC , family=binomial, data = ex2)
summary(modT)
modT <- stepAIC(modT, scope = fm , family=binomial, data =ex, k = 4)
summary(modT)
$\endgroup$
2
  • $\begingroup$ +1 for pointing out the absence of a random draw in the dgp. On the other hand the absence of an intercept in the dgp seems unproblematic; it should just be estimated to be zero... $\endgroup$ Apr 12, 2016 at 3:03
  • $\begingroup$ yep, I solved it with sapply(ex$prob, function (p){rbinom(1,1,p)}) $\endgroup$
    – GabyLP
    Apr 12, 2016 at 3:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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