Nonlinear regression and delta method example: --------------------------------------------- The file "child.dat" lists British census data on the percentage of women who are mothers (with at least one child), arranged by women's age and birthyear. The graph shows curves for each birthyear, that is, each curve represents one column of the table (Hamilton, 1990, pp. 168ff.). Women with the same birthyear are said to be in the same "birth cohort". > child _ read.table("child.dat",header=T) > child Age Y1920 Y1930 Y1940 Y1945 Y1950 Y1955 Y1960 1 15 0 0 0 0 0 0 0 2 20 7 9 13 17 19 18 13 3 25 39 48 59 60 53 45 39 4 30 67 75 82 82 75 68 NA 5 35 76 83 87 88 83 NA NA 6 40 78 86 89 90 NA NA NA 7 45 NA 86 89 NA NA NA NA > matplot(child$Age,as.matrix(child)[,-1],type="l") > legend(locator(),legend=names(child)[-1],lty=1:7) > title("Percent Mothers, by Women's Age and Birthyear") > title(xlab="Women's age",ylab="Percent") We will fit the Gompertz growth curve model to at least the 1940 data: y = f(x; alpha, beta, gamma) = alpha * exp ( -gamma * exp ( -beta * x )) where y = percent, x = age. Term of the log-likelihood: > nllterm _ formula ( ~ log(sig2)/2 + (Y1940 - alpha * exp( - gamma * + exp( - beta * Age)))^2/(2*sig2) ) Good starting values [after some trial and error]: alpha = 90 ( y ~ 90 as x -> infinity ) gamma = large, say 1000 ( y ~ 0 as x -> 0 ) beta = 0.3 ( e.g. 82 = f(30) ) sig2 = 0.1 ( trial and error with ms(), or see nls() solution below! ) ======================================================================== Finding MLE's with ms() ----------------------- > library(MASS) > nll3 _ deriv3(nllterm, c("alpha","beta","gamma","sig2"), + function(alpha, beta, gamma, sig2, Y1940, Age) NULL ) > ms1940 _ ms( ~ nll3(alpha,beta,gamma,sig2,Y1940,Age), data=child, + start=c(alpha=90,beta=0.3,gamma=1000,sig2=.1),trace=T) [ trace output omitted ] > summary(ms1940) Final value: -6.413099 Solution: Par. Grad. Hessian.alph Hessian.beta Hessian.gamm alpha 89.10386725 4.357045e-08 7.227430e+01 1.263514e+04 -5.097397e-01 beta 0.30956206 4.658676e-05 1.263514e+04 1.131245e+07 -5.189961e+02 gamma 941.97111362 -2.541274e-09 -5.097397e-01 -5.189961e+02 2.420364e-02 sig2 0.05887645 -2.127560e-06 -7.400318e-07 -7.912630e-04 4.316282e-08 Hessian.sig2 alpha -7.400318e-07 beta -7.912630e-04 gamma 4.316282e-08 sig2 1.009682e+03 Information: alph beta gamm sig2 alph 2.187500e-02 -2.030480e-04 -3.893241e+00 2.334103e-11 beta -2.030480e-04 7.329351e-06 1.528862e-01 -9.407097e-13 gamm -3.893241e+00 1.528862e-01 3.237646e+03 -2.144619e-08 sig2 2.334103e-11 -9.407097e-13 -2.144619e-08 9.904105e-04 Convergence: BOTH X- AND RELATIVE FUNCTION CONVERGENCE Computations done: Iterations Function Gradient 9 11 9 COMPARE WITH: > summary( ms(nllterm, data=child, + start=c(alpha=90,beta=0.3,gamma=1000,sig2=.1),trace=T) ) [output omitted] Standard errors of parameters: > sqrt(diag(summary(ms1940)$Information)) [1] 0.147901982 0.002707277 56.900320039 0.031470788 ======================================================================= Delta Method: Confidence Intervals, Etc. ---------------------------------------- Predicting a single new observation (at age 50): > fcn _ deriv3( ~ alpha * exp ( -gamma * exp ( -beta * x )), + c("alpha","beta","gamma"), function(alpha,beta,gamma,x) NULL ) > as.numeric(fcn(89.10387, 0.30956, 941.9711, 50)) [1] 89.08795 > print(grad _ attr(fcn(89.10387, 0.30956, 941.9711, 50),"gradient")) alpha beta gamma [1,] 0.9998214 0.7958231 -1.689698e-05 > summary(ms1940)$Information alph beta gamm sig2 alph 2.187500e-02 -2.030480e-04 -3.893241e+00 2.334103e-11 beta -2.030480e-04 7.329351e-06 1.528862e-01 -9.407097e-13 gamm -3.893241e+00 1.528862e-01 3.237646e+03 -2.144619e-08 sig2 2.334103e-11 -9.407097e-13 -2.144619e-08 9.904105e-04 > grad %*% summary(ms1940)$Information[1:3,1:3] %*% t(grad) [,1] [1,] 0.02167706 So an asymptotic 95% interval for the percent of women in the 1940 cohort who were mothers at age 50 would be > 89.08795 + c(-2,0,2)*sqrt(0.02167706) [1] 88.79349 89.08795 89.38241 Can do something similar to develop a "confidence envelope" for the curve: > x _ 15:45 > vals _ as.numeric(fcn(89.10387, 0.30956, 941.9711, x)) > grads _ attr(fcn(89.10387, 0.30956, 941.9711, x),"gradient") > varcov _ summary(ms1940)$Information[1:3,1:3] > vars _ grads %*% varcov %*% t(grads) > SEs _ sqrt(diag(vars)) > envelope _ vals + matrix(SEs,ncol=1)%*%matrix(c(-2,0,2),nrow=1) > matplot(x,envelope,type="l") Can't really see the envelope. Two explorations are worthwhile: 1. "Blow up" the picture by looking at just a piece of it. > matplot(x[15:30],envelope[15:30,],type="l") 2. Plot SE's along with data. > plot(x,SEs,type="l") > lines(x,as.numeric(fcn(89.10387, 0.30956, 941.9711, x)/360)) > points(child$Age,child$Y1940/360) [note: Ymax=360 and SEmax=.25, so divide y's by 4*90=360 to get on same scale...] =================================================================== NOTE: Can get to a similar place using the nls function, which directly solves nonlinear least squares problems using the Gauss-Newton method discussed in class [with a numerical approximation to the gradient]: > nls1940 _ nls(Y1940 ~ alpha * exp(-gamma*exp(-beta*Age)), + data=child,start=list(alpha=90,beta=0.30,gamma=1000),trace=T) > summary(nls1940) Formula: Y1940 ~ alpha * exp( - gamma * exp( - beta * Age)) Parameters: Value Std. Error t value alpha 89.103900 0.19578900 455.1010 beta 0.309562 0.00358636 86.3165 gamma 941.971000 75.35320000 12.5007 Residual standard error: 0.320989 on 4 degrees of freedom Correlation of Parameter Estimates: alpha beta beta -0.508 gamma -0.464 0.993 Now reconstruct the variance/covariance matrix from the correlation matrix and the SE's: > cormtx _ matrix(scan(),ncol=3) 1: 1 -0.508 -0.464 4: -0.508 1 0.993 7: -0.464 0.993 1 10: > SEnls _ c(0.19578900, 0.00358636, 75.35320000) > newvarcov _ diag(SEnls) %*% cormtx %*% diag(SEnls) > newvars _ grads %*% newvarcov %*% t(grads) > newSEs _ sqrt(diag(newvars)) > newenvelope _ vals + matrix(newSEs,ncol=1)%*%matrix(c(-2,0,2),nrow=1) > matplot(x,newenvelope,type="l") Compare the two sets of SE's: > plot(x,newSEs,type="l") > lines(x,SEs,type="l") =======================================================================