Math 534 Assignment 4

(a) Use the R function nls() to apply the Gauss-Newton algorithm to obtain your solution. Make trace = TRUE to see iterative process.

data <- read.table("http://mathfaculty.fullerton.edu/mori/data/math534/serum_conc.txt", header = TRUE)
attach(data)


model <- function(alpha1, lam1, alpha2, lam2, t){
  E1 = exp(-lam1*t)
  E2 = exp(-lam2*t)
  
  f = alpha1*E1 + alpha2*E2 
  grad = cbind(E1, -alpha1*t*E1, E2, -alpha2*t*E2)
  attr(f,'gradient') <- grad
  f
}

Result = nls(Concentration ~ model(alpha1, lam1, alpha2, lam2, Time),
             start = list(alpha1 = 50, lam1 = .1, alpha2 = 100, lam2 = .05),
             trace = TRUE, nls.control(maxiter = 50, tol = 1e-5, minFactor = 1/1024,
                                       printEval=TRUE))

## 14435.8 :  5e+01 1e-01 1e+02 5e-02
## 14383.48 :   36.0690291   0.1143425 114.3875112   0.0523250
## 14319.21 :   21.53251262   0.14563304 129.85323976   0.05499193
## 13925.86 :   13.56656654   0.21934644 139.74011507   0.05741538
## 12501.35 :   16.22588400   0.37405375 141.14086727   0.05963966
## 9751.762 :   28.61423511   0.54118840 137.31989424   0.06317254
## 5761.923 :   50.59357224   0.66859816 131.77035350   0.07124101
## 1762.429 :   80.00157736   0.79046988 128.79114805   0.09037282
## 82.76046 :   95.6995999   1.0045958 143.1358196   0.1346073
## 44.19967 :   80.3406140   1.2315740 162.7383946   0.1617354
## 34.41079 :   80.9331983   1.3141072 163.0048707   0.1621803
## 34.38039 :   81.2533200   1.3050516 162.5684482   0.1617586
## 34.38026 :   81.2396522   1.3061827 162.6017218   0.1617947
## 34.38026 :   81.2414872   1.3060523 162.5975661   0.1617902
## 34.38026 :   81.2412775   1.3060674 162.5980454   0.1617907

(b) Plot the data and the fitted curve. Superimpose the graph by the fitted function that is obtained by a one-term clearance

# plot(Time,Concentration)
# t = seq(0,50,.1)
# lines(t,81.2412775*exp(-1.3060674*t) +162.5980454*exp(-0.1617907*t) ,col='green')
# lines(t,211.9203*exp(-0.2357*t),col='red', lty = 2)
# legend(40, 150, legend=c("One-Term", "Two-Term"),
#        col=c("red", "green"), lty=c(2,1), cex=0.6)]
time_df = data.frame(Time, Concentration)
t = seq(0,50,.1)
one_term = 211.9203*exp(-0.2357*t)
two_term = 81.2412775*exp(-1.3060674*t) +162.5980454*exp(-0.1617907*t)
lines_df = data.frame(t, one_term, two_term)

ggplot(time_df, aes(Time, Concentration)) + 
  geom_point() +
  geom_line(data = lines_df, aes(t,one_term, color = "One-Term")) + 
  geom_line(data = lines_df, aes(t, two_term, color = "Two-Term")) + 
  labs(main = "", y = "Concentration", x = "Time", color  = "")

It is clear that the two component model is significantly better at fitting the data, than the one component model, given that the fitted line for the two-component model(green) goes through every data point while the one component model does not.

(a) Write an R function to implement iteratively reweighted least squares. Use Exercise J-2.1 part (c) as a guide for your output and convergence criteria. Use Iterative reweighted least squares to find the maximum likelihood estimates for \(b_0,b_1,b_2\).

dose <- c(0,2,10,50,250) #  x
num.ani <- c(111,105,124,104,90) # n_i
num.tum <- c(4,4,11,13,60) # y

X <- cbind(matrix(1,5,1), dose, dose^2) # matrix of data including
beta0 <- c(.038,.01,.000002) # Starting Values 

JacfW_poiss = function(beta,X, num.ani){
  xb = X%*%beta # 5 by 1
  f = num.ani*(1 - exp(-xb)) # 5 by 1
  J = diag(as.vector(exp(-xb))) %*% X*num.ani 
  W = diag(1/as.vector(f*exp(-xb)))
  list(f=f, J=J, W=W)
}

GN = function(y,X,num.ani, beta0, Jac, Wt = 1, maxit,IRLS = TRUE){
values <- matrix(0, ncol = 6, nrow = maxit, 
                   dimnames = list(rep("", maxit), c("IT", "b0", "b1", "b2", "MRE", "||Gradient||")))
it = 0
condition = FALSE
while (it < maxit & condition == FALSE){
  it = it +1
  a <- do.call(Jac, list(beta0, X, num.ani))
  J = a$J
  f = a$f
  if (IRLS==TRUE){
    Wt = a$W
  }
  JW = t(J) %*% Wt
  dir = solve(JW%*%J)%*%JW %*% (y-f)
  beta1 = beta0 + dir
  mre = norm((beta1-beta0)/pmax(1,abs(beta1)), type = "I")
  gradnorm = norm(JW %*% (y-f), type = "2")
  
  values[it,] <-c(format(it, digits = 0), formatC(beta0[1],digits = 12,format = "f"),
                     formatC(beta0[2],digits = 12,format = "f"), formatC(beta0[3],digits = 12,format = "f"),
                     format(x = mre, scientific = TRUE, digits = 2),
                     format(x = gradnorm, scientific = TRUE, digits = 2) ) 
  # if IRLS = TRUE, no step-having
  beta0 = beta1
  if(gradnorm < 10^(-9) && mre < 10^(-6)){
    condition = TRUE # If both conditions are met, tol <- FALSE breaks the while loop 
  }
  
}
#print(values[1:it,], quote = FALSE)


rownames(beta0) = c("b0", "b1", "b2")
colnames(beta0) = c("MLE")
return(list("values" = values[1:it,], "beta0" = beta0)) 
}

a <- GN(num.tum,X,num.ani,beta0,'JacfW_poiss', Wt=1, maxit = 30, IRLS = TRUE)

(b) Plot the probability of positive response for doses ranging from 0 to 250. Superimpose this plot by the observed proportions # of incidents/number tested.

prop <- num.tum/num.ani # proportion #incidents/#tested
a <- matrix(a$beta0,3,1)
x <- seq(0,250)
s <- cbind(rep(1,251), x, x^2)
p <- 1 - exp(-s%*%a)

(a) Use iteratively reweighted least squares to obtain the maximum likelihood estimates for the parameters.

data <- read.table("http://mathfaculty.fullerton.edu/mori/math534/examdata/blowBF.txt", header = TRUE)
attach(data)
y = data$y

incpt = matrix(1,length(y), 1)
X = cbind(incpt,log(data$D), data$S)

JacfW_poiss = function(beta,X){
  xb = X%*%beta # 659 by 1
  f = exp(xb)/(1+exp(xb)) # 659 by 1
  J = diag(as.vector(exp(xb)/(1+exp(xb))^2)) %*% X 
  W = diag(1/as.vector(exp(xb)/((1+exp(xb))^2)))
  list(f=f, J=J, W=W)
}

GN = function(y,X, beta0, Jac, Wt = 1, maxit,IRLS = TRUE){
  values <- matrix(0, ncol = 6, nrow = maxit, 
                   dimnames = list(rep("", maxit), c("IT", "b0", "b1", "b2", "MRE", "||Gradient||")))
  it = 0
  condition = FALSE
  while (it < maxit & condition == FALSE){
    it = it +1
    a <- do.call(Jac, list(beta0, X))
    J = a$J
    f = a$f
    if (IRLS==TRUE){
      Wt = a$W
    }
    JW = t(J) %*% Wt
    dir = solve(JW%*%J)%*%JW %*% (y-f)
    beta1 = beta0 + dir
    mre = norm((beta1-beta0)/pmax(1,abs(beta1)), type = "I")
    gradnorm = norm(JW %*% (y-f), type = "2")
    
    values[it,] <-c(format(it, digits = 0), formatC(beta0[1],digits = 12,format = "f"),
                    formatC(beta0[2],digits = 12,format = "f"), formatC(beta0[3],digits = 12,format = "f"),
                    format(x = mre, scientific = TRUE, digits = 2),
                    format(x = gradnorm, scientific = TRUE, digits = 2) ) 
    # if IRLS = TRUE, no step-having
    beta0 = beta1
    if(gradnorm < 10^(-9) && mre < 10^(-6)){
      condition = TRUE # If both conditions are met, tol <- FALSE breaks the while loop 
    }
    
  }
  #print(values[1:it,], quote = FALSE)
  rownames(beta0) = c("b0", "b1", "b2")
  colnames(beta0) = c("MLE")
  return(list('values' = values[1:it,], 'beta0' = beta0))
}

a <- GN(y,X,beta0,'JacfW_poiss', Wt=1, maxit = 30, IRLS = TRUE)

(b) Plot the three-dimensional graph of \(\pi(\beta)\), using your maximum likelihood estimates and as a function of \(X_1 = S\) and \(X_2 = \log(D)\), For your plot, \(X_1\) should range in the interval (0,1) and \(X_2\) should range in the range of observed value for \(\log(D)\).

# 3-D plot
d = seq(5,30) # range of diameter
s = seq(0,1,length.out = length(d)) # range of severity 
prob <- function(d,s){
  len = length(d)
  r = numeric(length(len))
  mleB <- c(-9.562085, 3.197564, 4.508593) # mleB obtained from GN
  for(i in 1:len){
    x = matrix(c(1,log(d[i]),s[i]),1,3)
    ex <- exp(x%*%mleB)
    r[i] <- ex/(1+ex)
  }
  r
}

threed <- outer(d,s,prob)


#install.packages("plot3D")
library(plot3D)

persp3D(x = s, y = log(d) , z = threed, col = "red", alpha = .5, xlab = "Severity", ylab = "Diameter", zlab = "Probability Y = 1", main = "")

(c) What does your model predict for the probability that a tree with a diameter of 10 located in an area with wind severity of 0.3 is blown down?

# Using the prob() function defined above, let d = 10, s = .3
p = prob(10,.3)