################################################## ## 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!