## @knitr lecture14prep library(emdbook) library(mvtnorm) library(ggplot2) library(bbmle) #load the wolf data wolves <- read.csv("./data/16e2InbreedingWolves.csv") ## @knitr likelihoodPrep set.seed(697) counts <- rpois(10, 15) lambdaVals <- 0:50 lik <- sapply(lambdaVals, function(x) prod( dpois(counts, lambda=x) ) ) # ## @knitr likelihoodMLE colvec<-rep("white", 50) colvec[counts]<-"red" mle<-lambdaVals[which(lik==max(lik))] poisCurvemle <- dpois(lambdaVals, mle) barplot(poisCurvemle, ylab="Probability Density", xlab="lambda", col=colvec, width=1, xlim=c(0,50), main=paste("MLE of Lambda = ", mle, sep="")) ## @knitr LoglikelihoodPlot lambdaVals<-lambdaVals[-1] lik <- sapply(lambdaVals, function(x) prod( dpois(counts, lambda=x) ) ) # ll <- sapply(lambdaVals, function(x) sum( dpois(counts, lambda=x, log=TRUE) ) ) plot(lambdaVals, ll, ylab="Log-Likelihood", xlab="lambda", pch=19, ylim=c(-60,-25)) ## @knitr LoglikelihoodConfPlot plot(lambdaVals, ll, ylab="Log-Likelihood", xlab="lambda", pch=19, ylim=c(-60,-25)) abline(h=max(ll)-1.92, col="red", lwd=2, lty=2) ## @knitr likelihoodCI lConf <- max(ll) - 1.92 lConf #Get which values of lambda that have LL >= 95% cutoff confVals <- lambdaVals[which(ll >= lConf)] confVals[c(1, length(confVals))] ## @knitr LRTest #compare our estimated fit to the hypothesis that lambda = 12 G12 <- 2*(ll[17] - ll[12]) pchisq(G12, 1, lower.tail=F) #compare our estimated fit to the hypothesis that lambda = 15 G15 <- 2*(ll[17] - ll[15]) pchisq(G15, 1, lower.tail=F) ## @knitr wolfMeanLL #load the wolf data wolves <- read.csv("./data/16e2InbreedingWolves.csv") meanVals <- seq(0, 0.5, .001) wolfLL <- sapply(meanVals, function(x) sum(dnorm(wolves$inbreeding.coefficient, mean=x, sd=0.1, log=TRUE))) ## @knitr wolfMeanLLMaxEst maxIDX <- which(wolfLL== max(wolfLL)) #the LL estimate of the mean meanVals[maxIDX] #comparison to mean estimate mean(wolves$inbreeding.coefficient) ## @knitr wolfCILLImage plot(wolfLL ~ meanVals, type="l", lwd=1.2, xlim=c(0.15, 0.3), ylim=c(0,22)) ciIDX <- which(abs(wolfLL-wolfLL[maxIDX]) <= 1.92) points(meanVals[maxIDX], wolfLL[maxIDX], col="blue", pch=19, cex=2) abline(h=wolfLL[maxIDX]-1.92, v=0, col="red", lty=2, lwd=1.4) ## @knitr wolfCILL ciIDX <- which(abs(wolfLL-wolfLL[maxIDX]) <= 1.92) #lowerCI meanVals[ciIDX[1]] #upperCI meanVals[ciIDX[length(ciIDX)]] ## @knitr wolfCompare #Compare to Mean = 0.3 idx3 <- which(meanVals==0.3) #calculate the Chisq LR G <- 2*(wolfLL[maxIDX] - wolfLL[idx3]) #put it to the test pchisq(G, 1, lower.tail=F) ## @knitr wolfTwoBrute sdVals <- seq(0.01, 0.4, 0.001) meanVals <- seq(0.01, 0.4, .01) vals <- as.data.frame(expand.grid(meanVals=meanVals, sdVals = sdVals)) getWolfL <- function(msVec) sum(dnorm(wolves$inbreeding.coefficient, mean=msVec[1], sd=msVec[2], log=TRUE)) vals$LL <- apply(vals, 1, getWolfL) library(ggplot2) ggplot(data=vals, mapping=aes(x=meanVals, y=sdVals, fill=LL)) + geom_raster(interpolate=T) + scale_fill_gradient(low="blue", high="red", limits=c(-10,22)) ## @knitr mleWolfMeanBruteForce sdVals <- seq(0.01, 0.4, 0.001) meanVals <- seq(-0.01, 0.5, .01) llmat <- matrix(rep(NA, length(sdVals) * length(meanVals)), nrow=length(meanVals)) for (i in 1:length(meanVals)) { for(j in 1:length(sdVals)) { llmat[i,j] <- sum(dnorm(wolves$inbreeding.coefficient, mean=meanVals[i], sd=sdVals[j], log=TRUE)) } } contour(meanVals, sdVals, llmat, levels=seq(-3,21,3), xlab="Mean", ylab="SD") ## @knitr mleWolfMeanBruteForceSurf persp(meanVals, sdVals, llmat, xlab="Mean", ylab="SD\n", zlim=c(0,22), zlab="Log-Likelihood") ## @knitr profileBrute profileFromMat <- function(rowCol) data.frame(t(apply(llmat, rowCol, function(x) c(ll = max(x), idx = which(x == max(x)))))) meanProfile <- profileFromMat(1) meanProfile$sdVal <- sdVals[meanProfile$idx] contour(meanVals, sdVals, llmat, levels=seq(-3,21,3), xlab="Mean", ylab="SD") lines(meanVals, meanProfile$sdVal, lwd=2, col="red") ## @knitr profileBrute2 sdProfile <- profileFromMat(2) sdProfile$meanVal <- meanVals[sdProfile$idx] contour(meanVals, sdVals, llmat, levels=seq(-3,21,3), xlab="Mean", ylab="SD") lines(meanVals, meanProfile$sdVal, lwd=2, col="red") lines(sdProfile$meanVal, sdVals, lwd=2, col="blue") ## @knitr mleWolfFunction #first write a function that you want to minimize wolfLL <- function(m, s) -sum(dnorm(inbreeding.coefficient, mean=m, sd=s, log=TRUE)) #now feed the function to an optimizer wolfLL_Fit <- mle2(wolfLL, data=wolves, start=list(m=0.2, s=0.2)) ## @knitr mleWolfFunctionWarn wolfLL_Fit <- mle2(wolfLL, data=wolves, start=list(m=0.2, s=0.2)) warnings() ## @knitr mleWolfFunctionFix #first write a function that you want to minimize wolfLL <- function(m, s){ #if an impossible value for a param #return NaN if(s <= 0) return(NaN) -sum(dnorm(inbreeding.coefficient, mean=m, sd=s, log=TRUE)) } #now feed the function to an optimizer wolfLL_Fit <- mle2(wolfLL, data=wolves, start=list(m=0.2, s=0.2)) ## @knitr mleWolvesOutput summary(wolfLL_Fit) ## @knitr mleWolvesMeanProfile mleWolvesNames <- mle2(inbreeding.coefficient ~ dnorm(mean=mean, sd=sd), data=wolves, start=list(mean=0.1,sd=0.1)) plot(profile(mleWolvesNames)) ## @knitr mleWolvesMeanCI confint(wolfLL_Fit) # Wald CIs assume quadratic relationship at peak confint(wolfLL_Fit, method="quad") ## @knitr mleWolvesNorm mleWolves <- mle2(inbreeding.coefficient ~ dnorm(mean=m, sd=s), data=wolves, start=list(m=0.1,s=0.1)) ## @knitr mleWolves2 mleWolvesPups <- mle2(pups ~ dpois(lambda=lambda), data=wolves, start=list(lambda=3)) ## @knitr mleWolvesRegression wolfRegLL <- function(slope, int, sdResid){ #in case of non-possible SD value, NaN if(sdResid <= 0) return(NaN) #fitted values as means fittedPups <- slope * inbreeding.coefficient + int -sum(dnorm(pups, mean = fittedPups, sd = sdResid, log=T)) } mleWolfLM <- mle2(wolfRegLL, data=wolves, start=list(slope=-10, int=6, sdResid =1.2)) ## @knitr mleWolvesRegressionCompare pup_lm <- lm(pups ~ inbreeding.coefficient, data=wolves) coef(mleWolfLM) coef(pup_lm) ## @knitr mleWolvesNull mleWolfNull <- mle2(wolfRegLL, data=wolves, start=list(int=6, sdResid =1.2), fixed=list(slope=0)) anova(mleWolfLM, mleWolfNull) ## @knitr PufferExercise #load the pufferfish data puffer <- read.csv("./data/16q11PufferfishMimicry Caley & Schluter 2003.csv") pufferFun <- function(slope, intercept, residSD){ #in case of non-possible SD value, NaN if(residSD <= 0) return(NaN) predatorsFit <- intercept + slope*resemblance -sum(dnorm(predators, predatorsFit, sd=residSD, log=T)) } puffer_mle <- mle2(pufferFun, start=list(slope=1, intercept=1, residSD=2), data=puffer) ## @knitr PufferExercise2 puffer_mle2 <- mle2(predators ~ dnorm(slope*resemblance + intercept, sd=residSD), start=list(slope=1, intercept=1, residSD=2), data=puffer) ## @knitr PufferExerciseSummary coef(summary(puffer_mle)) coef(summary(lm(predators ~ resemblance, data=puffer))) ## @knitr PufferExerciseCI confint(puffer_mle) ## @knitr PufferExerciseNHST puffer_mle_null <- mle2(predators ~ dnorm(intercept, sd=residSD), start=list(intercept=9, residSD=4), data=puffer) anova(puffer_mle, puffer_mle_null)