Likelihood Function

Generating some data with a know probability of germination -> myP

myP<-0.67
myData<-rbinom(n=10,size=1,prob=myP)
myData

 [1] 1 1 1 1 0 1 1 0 1 0

mean(myData)

[1] 0.7

Binomial Likelihood Function

nllBinomial<-function(parmEstimate,seedlingDat){

  nllik<- -sum(dbinom(seedlingDat,size=1,parmEstimate,log=TRUE))

  cat("nllik= ",nllik,sep=" ",fill=T)
  return(nllik)
}

nllBinomial(0.6,myData)

nllik=  6.324652
[1] 6.324652

nllBinomial(0.1,myData)

nllik=  16.43418
[1] 16.43418

Minimizing the function using optim

parVec<-c(0.1) # Initial parameter values 
outLogReg<-optim(par=parVec,fn=nllBinomial,method="L-BFGS-B",lower=c(0.01),upper=c(0.99),seedlingDat=myData)

nllik=  16.43418
nllik=  16.36786
nllik=  16.5012
nllik=  13.88586
nllik=  13.88586
nllik=  13.60701
nllik=  10.07198
nllik=  10.05039
nllik=  10.09367
nllik=  8.576073
nllik=  8.560991
nllik=  8.591218
nllik=  6.609054
nllik=  6.602835
nllik=  6.61531
nllik=  6.11829
nllik=  6.117376
nllik=  6.119248
nllik=  6.108907
nllik=  6.109091
nllik=  6.108772
nllik=  6.108643
nllik=  6.108661
nllik=  6.108673
nllik=  6.108643
nllik=  6.108667
nllik=  6.108667
nllik=  6.108643
nllik=  6.108667
nllik=  6.108667

outLogReg$par # 

[1] 0.6999994

Plotting the likelihood function

library(purrr)
probabilityVector<-seq(from=0, to=1,length=100)
llikOut<-map_dbl(probabilityVector,nllLogBinomial,seedlingDat=myData)

nllik=  Inf
nllik=  27.61133
nllik=  23.49347
nllik=  21.10213
nllik=  19.41792
nllik=  18.12139
nllik=  17.07024
nllik=  16.18858
nllik=  15.43111
nllik=  14.76861
nllik=  14.18114
nllik=  13.65448
nllik=  13.17813
nllik=  12.74411
nllik=  12.34625
nllik=  11.97963
nllik=  11.6403
nllik=  11.32504
nllik=  11.03117
nllik=  10.75646
nllik=  10.49901
nllik=  10.25723
nllik=  10.02972
nllik=  9.8153
nllik=  9.612923
nllik=  9.421683
nllik=  9.240782
nllik=  9.069513
nllik=  8.907252
nllik=  8.753443
nllik=  8.607588
nllik=  8.469244
nllik=  8.338013
nllik=  8.213534
nllik=  8.095486
nllik=  7.983578
nllik=  7.877546
nllik=  7.777153
nllik=  7.682186
nllik=  7.59245
nllik=  7.507772
nllik=  7.427994
nllik=  7.352976
nllik=  7.282591
nllik=  7.216728
nllik=  7.155287
nllik=  7.098182
nllik=  7.045338
nllik=  6.99669
nllik=  6.952185
nllik=  6.911779
nllik=  6.875441
nllik=  6.843146
nllik=  6.814881
nllik=  6.790644
nllik=  6.770441
nllik=  6.754288
nllik=  6.742212
nllik=  6.734252
nllik=  6.730456
nllik=  6.730885
nllik=  6.735611
nllik=  6.744721
nllik=  6.758314
nllik=  6.776508
nllik=  6.799433
nllik=  6.82724
nllik=  6.860099
nllik=  6.898203
nllik=  6.94177
nllik=  6.991044
nllik=  7.046301
nllik=  7.107854
nllik=  7.176056
nllik=  7.251305
nllik=  7.334054
nllik=  7.424822
nllik=  7.524196
nllik=  7.632856
nllik=  7.751582
nllik=  7.881283
nllik=  8.023017
nllik=  8.17803
nllik=  8.3478
nllik=  8.534097
nllik=  8.739062
nllik=  8.965317
nllik=  9.216123
nllik=  9.495597
nllik=  9.80904
nllik=  10.16344
nllik=  10.56828
nllik=  11.03683
nllik=  11.58856
nllik=  12.25368
nllik=  13.08276
nllik=  14.17066
nllik=  15.73034
nllik=  18.44139
nllik=  Inf

plot(probabilityVector,-llikOut,typ='l')

Problem: Fit a linear regression to height vs weight data via maximum likelihood weight~Normal(mu,sigma) mu<-intercept + B*height

library(rethinking)
data("Howell1")
myData<-Howell1
myData<-myData[myData$age>18,]
head(myData)

plot(myData$height,myData$weight)

This likelihood function assumes a sd->1

nllRegression<-function(parmVector,weightDat,heightDat){

  intercept<-parmVector[1]
  B<-parmVector[2]
  
  w<-weightDat
  h<-heightDat
  mu<-intercept + B*h
  
  nllik<- -sum(dnorm(x=w,mean=mu,sd=1,log=TRUE))

  cat("nllik= ",nllik,sep=" ",fill=T)
  return(nllik)}

nllRegression(c(1,1),weightDat=myData$weight,heightDat=myData$height)

nllik=  2121050
[1] 2121050

Fitting the linear regression using maximum likelihood

parVec<-c(0.1,0.1) # Initial parameter values 
outReg<-optim(par=parVec,fn=nllRegression,method="L-BFGS-B",lower=c(-Inf,-Inf),upper=c(Inf,Inf), weightDat=myData$weight,heightDat=myData$height)

nllik=  156667.4
nllik=  156657.2
nllik=  156677.6
nllik=  155083.2
nllik=  158260
nllik=  2716113
nllik=  2716157
nllik=  2716070
nllik=  2722825
nllik=  2709410
nllik=  4595.789
nllik=  4595.834
nllik=  4595.744
nllik=  4599.937
nllik=  4599.937
nllik=  4595.789
nllik=  4595.834
nllik=  4595.744
nllik=  4599.936
nllik=  4599.937
nllik=  4595.788
nllik=  4595.833
nllik=  4595.743
nllik=  4599.935
nllik=  4599.936
nllik=  4595.784
nllik=  4595.829
nllik=  4595.739
nllik=  4599.931
nllik=  4599.932
nllik=  4595.768
nllik=  4595.814
nllik=  4595.724
nllik=  4599.916
nllik=  4599.916
nllik=  4595.706
nllik=  4595.751
nllik=  4595.661
nllik=  4599.853
nllik=  4599.854
nllik=  4595.456
nllik=  4595.501
nllik=  4595.411
nllik=  4599.604
nllik=  4599.604
nllik=  4594.458
nllik=  4594.503
nllik=  4594.413
nllik=  4598.605
nllik=  4598.606
nllik=  4590.467
nllik=  4590.512
nllik=  4590.422
nllik=  4594.615
nllik=  4594.615
nllik=  4574.574
nllik=  4574.619
nllik=  4574.53
nllik=  4578.721
nllik=  4578.722
nllik=  4512.1
nllik=  4512.143
nllik=  4512.056
nllik=  4516.247
nllik=  4516.247
nllik=  4279.77
nllik=  4279.809
nllik=  4279.732
nllik=  4283.918
nllik=  4283.918
nllik=  3432.287
nllik=  3432.287
nllik=  3432.287
nllik=  3436.435
nllik=  3436.435
nllik=  3432.287
nllik=  3432.287
nllik=  3432.287
nllik=  3436.435
nllik=  3436.435
nllik=  3432.287
nllik=  3432.287
nllik=  3432.287
nllik=  3436.435
nllik=  3436.435
nllik=  3432.287
nllik=  3432.287
nllik=  3432.287
nllik=  3436.435
nllik=  3436.435
nllik=  3432.287
nllik=  3432.287
nllik=  3432.287
nllik=  3436.435
nllik=  3436.435

outReg$par # 

[1] -51.6295921   0.6251449

Plotting the fit

plot(myData$height,myData$weight)
abline(a=outReg$par[1],b=outReg$par[2])

PROBLEM 1: Modify the likelihood function to fit the sd to the data and replot the fit.

Solution:

library(rethinking)
data("Howell1")
myData<-Howell1
myData<-myData[myData$age>18,]
head(myData)

plot(myData$height,myData$weight)

Modify the likelihood to also fit sd

nllRegression<-function(parmVector,weightDat,heightDat,print=F){

  intercept<-parmVector[1]
  B<-parmVector[2]
  mySD<-parmVector[3]
  
  w<-weightDat
  h<-heightDat
  mu<-intercept + B*h
  
  nllik<- -sum(dnorm(x=w,mean=mu,sd=mySD,log=TRUE))

  if (print) cat("nllik= ",nllik,sep=" ",fill=T)
  return(nllik)}

Testing the likelihood function:

nllRegression(c(1,1,10),weightDat=myData$weight,heightDat=myData$height)

nllik=  22321.97
[1] 22321.97

Next, minimizing the negative log likelihood. By default optim performs minimization–that is why we use the NEGATIVE log likelihood

parVec<-c(0.1,0.1,10) # Initial parameter values 
outReg<-optim(par=parVec,fn=nllRegression,method="L-BFGS-B",lower=c(-Inf,-Inf, 0.01),upper=c(Inf,Inf, Inf), weightDat=myData$weight,heightDat=myData$height,print=F)
outReg$par # 

[1] -51.6295921   0.6251449   4.2428686

outReg

$par
[1] -51.6295921   0.6251449   4.2428686

$value
[1] 991.0056

$counts
function gradient 
      42       42 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

Plotting the fit

plot(myData$height,myData$weight)
abline(a=outReg$par[1],b=outReg$par[2])

PROBLEM 2: Fit the linear regression using the R lm function. Compare the parameter estimates.

Examining the data

head(myData)

lmFit<-lm(weight~height,data=myData)
summary(lmFit)

Call:
lm(formula = weight ~ height, data = myData)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.8344  -3.0747  -0.2775   2.8402  14.6456 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -51.62959    4.56320  -11.31   <2e-16 ***
height        0.62514    0.02947   21.21   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.255 on 344 degrees of freedom
Multiple R-squared:  0.5667,    Adjusted R-squared:  0.5655 
F-statistic:   450 on 1 and 344 DF,  p-value: < 2.2e-16

Question: How to the lm fits compare to the maximum likelhood fits? Note that lm uses ordinary least squares regression and not likelihood.

PROBLEM 3: Fit the linear regression using the Rethinking QUAP function. Compare the parameter estimates.

library(rethinking)
lm.qa <- quap(
    alist(
        weight ~ dnorm(mu,sd) ,  # normal likelihood
        sd ~ dexp(1),     # exponential prior
        mu<-intercept + beta*height,
        intercept~dnorm(0,10),
        beta~dnorm(0,10)
    ) ,
    data=list(weight=myData$weight,height=myData$height) )

precis(lm.qa)

Note the differences between the QUAP fit and the fits above. What do you think is causing this? Let’s look at a plot.

plot(myData$height,myData$weight)
abline(a=summary(lm.qa)$mean[2],b=summary(lm.qa)$mean[3])

Part of the problem is that the intercept and and beta are highly correlated so that different combinations of each fit the data well.

post <- extract.samples(lm.qa)
plot( intercept ~ beta , post , col=col.alpha(rangi2,0.1) , pch=16 )

Next steps. Could try to standardize data by subtracting mean of the DATA and dividing by the standard deviation of the DATA (see code below). These are Z scores. Compare fits of lm and quap using these new data. We will stop here and these would be exercises that you could pursue if interested.

weight<-myData$weight

Error: object 'myData' not found

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

plot(cars)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.