##################################################
## R code for section 4
## 2-19-2014
## Stephen Pettigrew
# Writing a Poisson log likelihood function into R:
loglike <- function(parameters, outcome, covariates){
# Add a column of 1's to your covariate matrix
# to account for the intercept
cov <- as.matrix(cbind(1, covariates))
# To make our code easier to read let's multiply
# the covariate matrix, X, by the parameters, beta
xb <- cov %*% parameters
sum(outcome * xb - exp(xb))
}
# Let's look test whether our function works with the
# curling data. The x covariate is whether the US
# had the hammer in that end
x <- c(0,1,0,1,0,1,0,1,0,1)
y <- c(1,0,1,2,0,1,0,2,0,1)
# Let's look at whether our function works by plugging
# the two parameters values into the function.
loglike(pars = c(1.2, .57), # these were arbitrarily picked
outcome = y,
covariates = x) # It doesn't give us an error so it works!
## Next week in section we'll talk about how to optimize functions,
## but just as a preview, let's find the MLE for the betas based on
## the data we generated
results <- optim(par = c(0,0),
outcome = y,
covariates = x,
fn = loglike,
method = "BFGS",
control = list(fnscale = -1, maxit=100000),
hessian = TRUE)
results$par ## There are our beta0 and beta1 estimates!