HW8comp

gradient and hessian matrix are used in calculating high-dimentional Newton-Raphson method, these are all route for us to find the best fit for the functions. Gradient is a column vector with first derivatives of f(x). Hessian matrix is square matrix filled with second-order partial derivatives of f(x). we want our model stable, and the uncertainty of the model and the correlation between the variables can be represented through variance-covariance matrix.

Ex 20-1:

1. Please use the Newton-Raphson method to find the maximum likelihood estimate (MLE) of the regression coefficients of logistic regression.

resp <- read.csv("resp.csv")

# A->1, P->0
resp$treatment2 <- ifelse(resp$treatment == "A", 1,
                          ifelse(resp$treatment == "P", 0, NA))


X <- cbind(rep(1,length(resp$outcome)),resp$treatment2, resp$age, resp$baseline)
Y <- resp$outcome

ftn <- function(betacoef){
        pi1 <- exp(X%*%betacoef)/(1+exp(X%*%betacoef))
        gradient <- t(X)%*%(Y-pi1)
        hessian <- -t(X)%*%diag(c(pi1*(1-pi1)),length(resp$outcome))%*%X
        Loglike <- sum(Y*log(pi1/(1-pi1))+log(1-pi1))
        return(list(gradient, hessian, Loglike))
}

(beta <- newtonraphson(ftn, c(0,0,0,0)))

## Algorithm converged

##             [,1]
## [1,] -0.79805332
## [2,]  1.23475884
## [3,] -0.01140389
## [4,]  1.98241179

hh <- glm(outcome ~ treatment2+age+baseline, data= resp, family = binomial)
hh$coefficients

## (Intercept)  treatment2         age    baseline 
## -0.79805332  1.23475884 -0.01140389  1.98241179

2.Please find the variance-covariance matrix for the betas .

solve(-ftn(beta)[[2]])

##              [,1]          [,2]          [,3]          [,4]
## [1,]  0.109437745 -3.062892e-02 -2.201856e-03 -2.640573e-02
## [2,] -0.030628915  5.133126e-02  3.292548e-05  1.351023e-02
## [3,] -0.002201856  3.292548e-05  6.630252e-05 -2.084783e-05
## [4,] -0.026405731  1.351023e-02 -2.084783e-05  5.426992e-02

vcov(hh)

##              (Intercept)    treatment2           age      baseline
## (Intercept)  0.109437743 -3.062891e-02 -2.201856e-03 -2.640573e-02
## treatment2  -0.030628914  5.133126e-02  3.292548e-05  1.351023e-02
## age         -0.002201856  3.292548e-05  6.630252e-05 -2.084783e-05
## baseline    -0.026405731  1.351023e-02 -2.084783e-05  5.426992e-02

3. Please find the log likelihood at the betas .

ftn(beta)[[3]]

## [1] -247.9434

logLik(hh)

## 'log Lik.' -247.9434 (df=4)