######################################################################################################################################### ##Question 1 - Who gets preventative care for colon cancer? ######################################################################################################################################### #################################Import the data######################################################################################### bfrss2010 = read.table( "http://www.montana.edu/ebelasco/ecns562/homework/bfrss2010.dat", header = T) names(bfrss2010) # [1] "COLON2" "HLTHPLAN_BIN" "white" "black" "AI" # [6] "Asian" "BMI4" "CollegeGrad" "HSGrad" "LessThanHS" #[11] "LowIncome" "MidIncome" "HighIncome" "AGE" "Rural" #[16] "rand" #################################Import two functions for OLS and Logit################################################################## source(file="L:/ECNS562/R/ols.R") source(file="L:/ECNS562/R/logit.R") #################################Set X matrix for use#################################################################################### X = as.matrix(cbind(1, bfrss2010[, c('AGE', 'HLTHPLAN_BIN', 'black', 'AI', 'Asian', 'BMI4', 'CollegeGrad', 'HSGrad', 'LowIncome', 'HighIncome', 'Rural')])) ##################################OLS Estimation######################################################################################### my_ols = ols(X, bfrss2010$COLON2) cbind(my_ols$beta, my_ols$se, my_ols$t) my_ols$R_sq #################################Logit Estimation using ################################################################################# b0 <- my_ols$beta; my_logit <- optim(b0, logit, x=X, y=bfrss2010$COLON2, hessian=TRUE, control=list(maxit=10000)) ################################Compute standard errors, t-values, and LL################################################################ my_logit_se <- sqrt(diag(solve(my_logit$hessian))) my_logit_t <- my_logit$par/my_logit_se LL <- -my_logit$value ################################Restricted model and LRI################################################################################# my_logit0 <- optim(as.matrix(.01), logit, x=as.matrix(X[,1]), y=bfrss2010$COLON2, hessian=TRUE, method="BFGS", control=list(maxit=10000)) LL0 = -my_logit0$value LRI = 1 - (LL/LL0) ################################Compute Marginal effects################################################################################# lambda = exp(X%*%my_logit$par)/(1+exp(X%*%my_logit$par)) #lambda = plogis(X%*%my_logit$par) mean.lambda = mean(lambda) marg = mean.lambda*(1-mean.lambda)*my_logit$par marg ######################################################################################################################################### ##Question 2 - Discrete choice model regarding travel alternatives ######################################################################################################################################### ###############################Import the data########################################################################################### data("TravelMode",package="AER") choice1 <- as.numeric(TravelMode$choice)-1 choice2 <- t(matrix(choice1,4,210)) Ttme <- t(matrix(TravelMode$wait,4,210)) Invc <- t(matrix(TravelMode$vcost,4,210)) Invt <- t(matrix(TravelMode$travel,4,210)) GC <- t(matrix(TravelMode$gcost,4,210)) Hinc <- t(matrix(TravelMode$income,4,210)) PSize <- t(matrix(TravelMode$size,4,210)) ones <- matrix(1,210,1) const.air <- cbind(ones,0,0,0); const.train <- cbind(0,ones,0,0); const.bus <- cbind(0,0,ones,0); ###############################construct X matrix and inputs################################################################################# x <- cbind(GC,Ttme,const.air,const.train,const.bus) n <- nrow(x); j <- 4; k <- ncol(x)/j; ###############################Conditional Logit function##################################################################################### source(file="L:/ECNS562/R/mnlogit.R") ###############################Optimize CL function and obtain statistics#################################################################### my_mnlogit <- optim(c(.1,.1,.1,.1,.1), mnlogit, x=x, y=choice2, hessian=TRUE, control=c(maxit=10000)) my_mnlogit_se <- sqrt(diag(solve(my_mnlogit$hessian))) my_mnlogit_t <- my_mnlogit$par/my_mnlogit_se cbind(my_mnlogit$par, my_mnlogit_se, my_mnlogit_t) #################################Version 2 of CL function to obtain Expected Utility########################################################## mnlogita <- function(b,x,y) { V <- matrix(0,n,j) for (h in 1:j) { V[,h] <- b[1]*x[,h] + b[2]*x[,(j+h)] + b[3]*x[,(2*j+h)] + b[4]*x[,(3*j+h)] + b[5]*x[,(4*j+h)] } P <- exp(V)/apply(exp(V),1,sum) LL <- -sum(y*log(P)) return(V) } ##################################Compute Baseline EU########################################################################################## V0 <- mnlogita(my_mnlogit$par,x=x,y=choice2) P0 <- exp(V0)/apply(exp(V0),1,sum) expV0 <- exp(V0) EV0 <- log(apply(expV0,1,sum)) alpha <- -my_mnlogit$par[1] ##################################Marginal impacts############################################################################################## marg_at_air <- mean(P0[,1]*(1-P0[,1])*my_mnlogit$par[2])*100 marg_at_train <- mean(P0[,2]*(0-P0[,1])*my_mnlogit$par[2])*100 marg_at_bus <- mean(P0[,3]*(0-P0[,1])*my_mnlogit$par[2])*100 marg_at_car <- mean(P0[,4]*(0-P0[,1])*my_mnlogit$par[2])*100 cbind(marg_at_air, marg_at_train, marg_at_bus, marg_at_car) #################################Compute new EU with additional driving costs##################################################################### new.x <- x; new.x[,4] <- x[,4]*1.25; V1 <- mnlogita(my_mnlogit$par,x=new.x,y=choice2) P1 <- exp(V1)/apply(exp(V1),1,sum) expV1 <- exp(V1) EV1 <- log(apply(expV1,1,sum)) #################################How many are predicted to choose each############################################################################## car.choice <- matrix(0,n,2) max0 <- apply(P0,1,max); max1 <- apply(P1,1,max); for(i in 1:n){ if(P0[i,4]