rm(list=ls())
n = 20; beta0 = 0; beta1 = 1
x = rnorm(n) # Standard normal
x = round(x,2)
xb = beta0 + beta1*x
p = exp(xb)/(1+exp(xb))
cbind(x,xb,p)
y = rbinom(n,1,p) # Bernoulli
cbind(x,xb,p,y)


n = 200; beta0 = 0; beta1 = 1
x = rnorm(n) # Standard normal
x = round(x,2)
xb = beta0 + beta1*x
p = exp(xb)/(1+exp(xb))
y = rbinom(n,1,p) # Bernoulli

summary(glm(y ~ x, family=binomial))
summary(lm(y ~ x))



nsim = 2000
pvals = matrix(NA,nsim,2)
colnames(pvals) = c("Logistic","Linear")

for(j in 1:nsim)
    {
    x = rnorm(n) # Standard normal
    x = round(x,2)
    xb = beta0 + beta1*x
    p = exp(xb)/(1+exp(xb))
    y = rbinom(n,1,p) # Bernoulli
    sumlogr = summary(glm(y ~ x, family=binomial))
    pvals[j,1] = sumlogr$coefficients[2,4]
    sumlinr = summary(lm(y ~ x))
     pvals[j,2] = sumlinr$coefficients[2,4]
    }