# SCS 2012: Mixed Models with R/Lab 3

### From Wiki1

### ### ### Lab Session Section 3: Generalized linear mixed models and related topics ### ## ## Related Topics: ## #* Accelerated Longitudinal Designs and age-period-cohort linear confounding #* Modeling seasonal and periodic effects with Fourier Analysis #* Using general splines to model effect of age or time #* Linear, quadratic, cubic and natural cubic splines #* General spline generator: splines with arbitrary degrees and smoothness #* Defining hypothesis matrices and using Wald tests to explore splines #* Plotting splines and spline features with confidence bounds #* Plotting log-odds or probabilities #* Interpreting hypothesis tests using confidence bounds #* Bonferroni and Scheffe confidence bound adjustment factors #* Testing non-linear cohort effects #* Alternatives to glmmPQL: lmer, glmmML,GLMMGibbs, library( spidadev ) library( p3ddev ) # Also need later: # # library( GLMMGibbs ) # library( lme4 ) # library( glmmML ) # # # Data: Respiratory infections among pre-school children in Indonesia # - 274 children were oberved at community clinics up to 6 times 3 months apart # - Baseline age ranges from 4 to 80 months # Question: Experience suggests that boys get more infections than # girls especially around 3 to 4 years of age. # # This data set is particularly interesting because it has a structure very # similar to Statistics Canada longitudinal surveys like the NPHS and the # NLSCY. There are 6 waves of data with subjects entering at different ages. # # This kind of design has been called an "Accelerated Longitudinal Design" # -- also a cross-sequential design or a cohort design (Miyazaki and # Raudenbush, 2000). Its competitors are the cross-sectional design in which # individuals of different ages are observed at one point in time, and a # single cohort design in which a single cohort is followed over a very # long time. # # In the cross-sectional design, age effects and cohort effects are perfectly # confounded. # # In the single cohort design, age effects and period effects are perfectly # confounded. # # In the ALD (accelerated longitudinal design) no pair of age, cohort.birth or period # are perfectly confounded but the LINEAR effect of any one is confounded with # a combination of the linear effects of the other two, in that the predictors # obey: # # period = cohort.birth + age # # so that there are only two degrees of freedom to estimate three parameters. # Note that we let period = time of measurement, cohort = time of birth. # # If we use all three variables, the model will be singular. Thus we drop one # of the three. If we drop 'cohort' on the basis that the effect of cohort # is more 'distal' than that of 'period' or 'age', we get the following # where we use 'beta[x]' for the 'true' but unknowable effect of 'x' and # 'beta.estimated[x]' for its estimated value in a particular model. # # combined effect of period, cohort and age: # # = beta[period] x period + beta[cohort.birth] x cohort.birth + beta[age] x age # # = beta[period] x period + beta[cohort.birth] x (period - age) + beta[age] x age # # = (beta[period] + beta[cohort.birth]) x period + (beta[age] - beta[cohort.birth]) x age # # So # beta.estimated[period] = beta[period] + beta[cohort.birth] # and # beta.estimated[age] = beta[age] - beta[cohort.birth] # # So depending on our assumption about beta[cohort.birth] we get: # # beta[period] = beta.estimated[period] - beta[cohort.birth] # and # beta[age] = beta.estimated[age] + beta[cohort.birth] # # In other words, the same data can be explained by a 0 effect of cohort # or by a non-zero effect of cohort with corresponding adjustments to the # estimated effect of age and period. # # Note that if cohort is identified by 'age at baseline' then the # sign of the cohort effect is the reverse of the sign if it is measured # by data of birth. So: # # beta[period] = beta.estimated[period] + beta[cohort.baseline.age] # and # beta[age] = beta.estimated[age] - beta[cohort.baseline.age] # # # It may seem paradoxical that the separate effects of age, period and cohort # cannot be estimated and the reason is that they, in fact, can be. # 'Non-linear' effects can be estimated separately under the assumption # that there are no interactions, i.e. that the separate effects are # additive. # It is only the linear effects that are fully confounded. # # Miyazaki and Raudenbush (2000) propose an approach using mixed models to # test wether cohort effects can be neglected. Essentially their tests ask # whether each cohorts trajectory 'splices' into those of adjoing cohorts. # If so the entire age trajectory can be explained by age effects only. # Their approach, however, assumes that there are no period effect, an effect # that can be irregular and potentially strong in my limited experience. # # We propose later to test a model that includes period effects and test # for possible effect of cohort beyond linear. If higher order cohort # can be ignored, then the estimation of an age trajectory is simplified. # # # # A first look at the data: # data(Indonesia) ?Indonesia ind <- Indonesia xqplot( ind ) scatterplot.matrix( ind , ell = TRUE ) indu <- up( ind, ~ id ) xqplot( indu ) # number of observations per child: ind$n <- capply( ind, ~ id, nrow ) indu <- up( ind, ~ id ) xqplot( indu ) scatterplot.matrix( indu, ell = TRUE ) inda <- up( ind, ~ id , all = T ) xqplot( inda ) scatterplot.matrix( inda, ell = T ) # # QUESTIONS: # How would you interpret # 1) the correlation between n and time # 2) the tapering of the season vs n plot # # Notes: # Although there are 1154 observations on 274 subjects, there are only tab( ind, ~ resp) # 107 instances of respiratory infections in tab( up( ind, ~ id, all=T)$resp > 0 ) # 81 subjects. For a dichotomous response, the 'effective sample size' # (Harrell, 2001) is the number of cases in the smaller group, here 107. # # Despite the seemingly large number of observations, we need to # interpret results cautiously. # # # Looking at marginal relationships between resp and age # # # Note: I do not carry out the following sequence of plots when # I'm analyzing data unless I'm preparing a graphic for publication # or I'm trying to impress someone. For exploratory analysis # most plots are created with 1 or 2 lines of code. # # The following sequence shows how you might go # from a simple plot for quick visualization to a more complex plot. # xyplot( resp ~ age, ind ) xyplot( resp ~ age, ind , panel = function(x,y,...) { panel.xyplot( x, y, ...) panel.loess( x, y, ...) } ) # fails for a subtle reason: the 'default' in panel.loess # family = "symmetric" which tries to fit a ROBUST estimate # of location. Here, however, a robust estimator is precisely # what we don't want. We want to fit the curve with an # ordinary mean to estimate the probability of 'resp'. # We get that by specifying family = "gaussian" # ?panel.loess xyplot( resp ~ age, ind , panel = function(x,y,...) { panel.xyplot( x, y, ...) panel.loess( x, y, ..., family = 'gaussian') } ) # focus scale of y to see predicted proportion of resp xyplot( resp ~ age, ind , ylim = c(0, .2), panel = function(x,y,...) { panel.xyplot( x, y, ...) panel.loess( x, y, ..., family = 'gaussian') } ) # create a factor for sex ind$Sex <- factor( ifelse( ind$sex == 0, "Male", "Female")) xyplot( resp ~ age | Sex, ind , ylim = c(0, .2), panel = function(x,y,...) { panel.xyplot( x, y, ...) panel.loess( x, y, ..., family = 'gaussian') } ) td( lwd = 2 ) # by specifying 'lwd' here instead of as an argument to # to xyplot, the change in lwd will be applied to # the key as well as to the main graph. # But it will persist for subsequent graphs. # xyplot( resp ~ age , ind , ylim = c(0, .2), groups = Sex, panel = panel.superpose, panel.groups = function(x,y,...) { panel.loess( x, y, ..., family = 'gaussian') }, auto.key = list( columns = 2, lines = T, points = F) ) # # Analysis # # Issues: # # -- data are longitudinal clustered in children so obs. are not independent # # -- the trajectories are not linear or quadratic, each might have a # a different degree as a polynomial # # -- age is confounded with time, hence possible period effects # # -- age and period are confounded with cohort, i.e. 'baseline_age' # # # We used a generalized linear mixed model fitted by penalized quasi-likelihood # # Essentially, we fit a sequence of 'lme' models with changing weights # much as a 'glm' involves fitting a sequence of 'lm' models with changing # weights. # # With true 'glm's the result is equivalent to maximizing # the likelihood of the the glm model. # # With 'glmm' models it is not, # hence PQL = penalized QUASI likelihood. # We cannot do LRTs to compare models but can use Wald tests. # # For binary outcomes with few observations per cluster, this approach # can yield biased estimates, particularly of variance parameters. # So we interpret results cautiously. # # The advantage of glmmPQL is that it is relatively fast and convergent # and its approach is not too difficult to understand since it parallels # lm and glm. Also glmmPQL allows correlation structures and realively complex # structures for the G matrix, allowing, for example, the fitting of # smoothing splines. # # One alternative is the 'lmer' function in 'lme4' that has a superior # fitting algorithm but lacks some of the other features of glmmPQL. # Currently it also lacks many methods. # library(MASS) # for glmmPQL fit <- glmmPQL ( resp ~ time + season + age + sex + height_for_age , data = ind, family = binomial, random = ~ 1 | id) summary( fit ) # # QUESTION: anything wrong? # fit <- glmmPQL ( resp ~ time + factor(season) + age * sex + height_for_age , data = ind, family = binomial, random = ~ 1 | id) summary(fit) wald( fit, "season") # straddling the Equator doesn't eliminate seasonal effects!! # QUESTIONS: # 1) What happened to season? it went from not significant to highly significant. # 2) Note that sex and its interaction with age are not significant. Do we quit # looking for sex effects? Or why not? # ind$Season <- factor( c("Spring","Summer","Fall","Winter") [ ind$season ] ) ind$Season <- reorder( ind$Season, ind$season) # QUESTION: Why are we doing this? tab( ~ Season, ind) # QUESTION: Why so imbalanced? tab( ~ Season + time, ind ) # We have 5 degrees of freedom for time (WHY?) and we've used 4: # 1 for linear time and # 3 for season. # We could try one more before moving on: fit2 <- glmmPQL ( resp ~ time + I(time^2) + Season + age * Sex + height_for_age , data = ind, family = binomial, random = ~ 1 | id) summary( fit2 ) wald( fit2, 'time' ) # # QUESTION: Should we drop quadratic time? # I would lean to drop because I'm concerned about the size of # the eventual model. # # So our provisional model is: fit3 <- glmmPQL ( resp ~ time + Season + age * Sex + height_for_age , data = ind, family = binomial, random = ~ 1 | id) summary( fit3 ) wald( fit3, 'Season' ) wald( fit3, 'Sex' ) wald( fit3, list( Interaction = ":")) # # Modeling Seasonal effects # # # In this example time is discrete as if all subjects were observed # nearly at the same time on each occasion. Often time is more # continuous, measured to the day or to the month. # # When the time scale of a study is a substantial part of a year # or stretches over many years, and the phenomemon might include # a periodic (seasonal) pattern that is repeated # each year, you can consider modeling periodic effects. # # Another common context for periodic effects is modeling # diurnal periodic patterns. Circadian cycles that are not # locked into a fixed period can sometimes be approximated # by diurnal models or modeled through autoregressive models # on the R side of a mixed model. # # The basic tool for periodic patterns is Fourier Analysis. # Fourier Analysis is to periodic effects what # exponential models are to asymptotic effects and what # lines and polynomials are to general smooth effects. # # Fourier analysis sounds extremely advanced but it's only # slightly more complicated than fitting polynomials. # # Fourier analysis involves fitting a sine 'wave' to the data, a wave # whose period is equal to the basic period (e.g. year, day). # This fundamental wave corresponds to a line in general analysis. # If a single sine wave doesn't do the job we consider adding # higher harmonics: sine waves whose periods are increasingly # smaller whole fractions of the basic period, i.e. 1/2 year, # 1/3 year, 1/4 year, 1/5 year and so on. As in fitting polynomials # we can keep going until the is no significant improvement. # # Just like it takes at least 2 points to fit a line, 3 points to fit # to fit a quadratic,4 points to fit a cubic and so on, # it takes takes at least # # 3 points within each period to fit a fundamental wave # 5 points ... to fit a fundamental + a first harmonic ( 1/2 year) # 7 points ... to fit up to the second harmonic (1/3 year) # 9 points ... to fit up to the third harmonic (1/4 year) # ... etc # # In our data we have 4 points so we can't go higher than the # fundamental. The fundamental is not a 'saturated' model for the # effect and we can test whether time used as a factor (the # saturated model) gives a better fit. # # The equivalent of terms in a polynomial (i.e. x, x^2, x^3, etc. # are sine,cosine pairs. So additional terms come in twos # # Suppose we have data measured in days stretching from -100 to 600, # almost two years. days <- -100:600 # Since we have 365 points # per year we could go very high with our harmonics but we # just look at the kinds of patterns generated by the first # few harmonics. # # For a first look, we won't fit a model but we will look at # what models look like. A bit like looking at a straight line # then a quadratic, then a cubic, etc. # # Note that our time points DO NOT have to be evenly spaced NOR # at the same times in each period. # # The period of the sin or cos function is 2*pi 2*pi # The period for our data is 365.25 days. (The .25 is only # necessary if our data stretches over such a long time that # ignoring leap years would cause the model to get out of phase, # but there's no harm in illustrating how easy it is to # incorporate a fractional period. # # To produce a sine curve with a period of 365.25 days, we need # to change our time scale from that of the data to that of the # the sin function. We do this by dividing by the period of # the data and multiplying by the period of the sin function. plot( days, sin( 2*pi* days / 365.25) , type = 'l') abline( h = 0, col = 'gray' ) abline( v = 365.25*(-1:2), col = 'gray') # When we fit a model, sin( 2*pi*days/365.25) will be # a regressor. Fitting the model will fit a coeffiecient for # this regressor, for example, 0.345 lines( days, 0.345 * sin( 2*pi* days / 365.25), col = 'green3') # A single sin curve will always have value 0 at 0, i.e. # its 'phase' must be zero. How could we fit other possibilities # for the phase (the first time when the curve is equal to 0). # Obviously we're talking about two parameters here: # the amplitude and the phase # We could use a formally non-linear model but the magic # of Fourier analysis is that we can accomplish this by # combining our sin curve with a cos curve of the same period # but possibly different amplitudes; lines ( days, 0.789 * cos( 2*pi* days / 365.25), col = 'blue') # Note that the cos curve is a sin curve with phase # equal to 1/4 year. # # When we add them we get another 'sine' curve with some other # phase: lines( days, 0.345 * sin( 2*pi* days / 365.25) + 0.789 * cos( 2*pi* days / 365.25), col = 'red') # The amplitude of the resulting curve is sqrt( 0.345^2 + 0.789^2) # and its (a) phase is 365.25 * atan2( -0.345, 0.789) # Fitting sine curves ('sine' generic for 'sin' or 'cos' also # called trigonometric curves) fit2.sin <- glmmPQL ( resp ~ time + sin( 2*3.1416 * (time + 6) / 4) + cos( 2*3.1416 * (time + 6) / 4) + age * Sex + height_for_age , data = ind, family = binomial, random = ~ 1 | id) summary( fit2.sin ) # QUESTIONS: # 1 ) Why did we add 6 to time? What happens if we don't. # 2 ) The 'sin' term is highly significant but the 'cos' term is # not. Should we drop the 'cos' term. # # The answer to the second question is too important to leave up # in the air: # If we drop the 'cos' term, then the 'sine' curve is force to # have a phase of 0, i.e. the curve = 0 at time = 0. # Thus, the model would not be invariant if we add a constant # to the time variable and the model would violate a principle of # invariance. # # Another way of thinking about it is that each term is # 'marginal' to the other and dropping one without the other # violates the principle of marginality. Except in very rare # circumstances, should be added or dropped as a pair. This # means that we need to carry out a wald test with 2 numDFs: wald ( fit2.sin, "sin|cos") # Of course we knew this would be significant. # # If we had more time points, we could try adding higher harmonics, # the next being: lines( days, sin(2 * 2 * 3.1416 * days / 365.25), lty = 2, col = 'green3') lines( days, cos(2 * 2 * 3.1416 * days / 365.25), lty = 2, col = 'blue') # Note how easy it is to generate each harmonic. You just multiply # the whole argument of the 'sin' or 'cos' by the appropriate integer: # 2 for the first harmonic # 3 for the second, etc. # # Alas we have too few points per year to illustrate this here. # If we did, we could keep adding other harmonics until we decide # to stop. Note that lower harmonics are treated as 'marginal' to # higher harmonics, just like polynomials. # # To make the process much cleaner and easier, we can define a # 'Sine' function that generates the matrix with our cos,sin pair # incorporating phase shifts and period scaling: Sine <- function( x ) cbind( sin = sin( 2 * 3.1416 * (x + 6) / 4), cos = cos( 2 * 3.1416 * (x + 6) / 4)) # to see what this does: time.ex <- seq( -1, 6, .1) head( Sine( time.ex ) ) matplot( time.ex, Sine( time.ex ), type = 'l') # generating higher harmonics is a cinch: matlines( time.ex , Sine( 2 * time.ex ), col = 'blue') matlines( time.ex , Sine( 3 * time.ex ), col = 'red') # combining all this with random coeffients produces white noise: try it a few times: plot( time.ex, cbind( Sine(time.ex) , Sine(2*time.ex), Sine(3 * time.ex)) %*% rnorm( 6 ), type = 'l') # Back to Indonesia: (we can fit the same model with:) Sine <- function( x ) cbind( sin = sin( 2 * 3.1416 * (x + 6) / 4), cos = cos( 2 * 3.1416 * (x + 6) / 4)) fit2.sin <- glmmPQL ( resp ~ time + Sine ( time ) + age * Sex + height_for_age , data = ind, family = binomial, random = ~ 1 | id) summary( fit2.sin ) wald( fit2.sin, 'Sine') # QUESTIONS: # 1) Sine( time ) is a function of 'time' NOT 'Season'. # Is is appropriate to think of it as modeling a # a Season effect. # 2) With this model, what is the interpretation of the # coefficient for the intercept term. # # There are 3 degrees of freedom for Season (4 levels - 1). # Our Sine function uses 2. # We'd like to test the next harmonic by adding the term: # # ... + Sine(2*time) + ..... # # but this would require 4 degrees of freedom. # How can we test whether we need to include the omitted # degree of freedom in this case? Just include the 'cos' # term of the next highest harmonic. Note: very rarely # the cos term won't work and you should try again with # the sin term. # fit2.sinS <- update( fit2.sin , . ~ . + Sine(2*time)[,2]) summary(fit2.sinS) # Turns out to be highly significant so we need all available # degrees of freedom to model Season. # # So we might as well just use the factor Season. # # This suggests that if we had more time points we would have # needed at least the first harmonic in our model. # # With the same number of subjects and the same number of # occasions for each subject BUT subjects coming in at different # dates, we could have used the variability in dates to # fit and test higher harmonics. # # So we are back where we started BUT much wiser! # EXERCISE: # # Try fitting a periodic curve to the sunspot data # in the 20th century # Note: we will use part of the same data later for splines. # data(sun) # monthly sunspot data in 20th century plot( spots ~ year , sun) # This is not strictly periodic data because the phase is # not fixed but it wanders randomly. # Trying a model with a period of 9.5 years plus a few harmonics Sine <- function( year ) cbind( sin = sin( 2*3.1416*year), cos = cos( 2.*3.1416*year)) # NOTE: By using a period of 1 in the Sine function we set things up so # that when we use the Sine function # we DIVIDE BY THE PERIOD and MULTIPLY FOR HARMONICS: fitsun <- lm( spots ~ Sine( year / 9.5) + Sine( 2 * year / 9.5) + Sine( 3 * year / 9.5), sun) summary( fitsun ) lines( predict( fitsun) ~ year, sun) # EXERCISE: # # Look at the graph, particularly at the peaks of the fitted curve # Is the fit too slow or too fast? # Adjust the period to see whether you can find one that fits better> # # # Modeling age with splines # # Splines: a very simple idea # # Problem with fitting a single polynomial: # # - Might need a very high degree for a good fit # # - But high degree polynomials go wild away from data # (Runge's phenomenon) # # - Points have influence at a distance -- generally not reasonable # # Spline: # # - split the range into disjoint intervals (e.g. 3 intervals) # # - the points at the end points are called 'knots' # ( 3 intervals => 2 knots ; p intervala => p-1 knots) # # - fit a separate polynomial in each interval # # - BUT do it in such a way that the polynomials join up # at the knots. # # - How? different 'degrees' of joining up: # - 0: just continuous but perhaps a kink # - change in slope = 1st derivative # - 1: smooth: no kink but perhaps a sudden change in curvature # - change in 2nd derivative = curvature = acceleration # - 2: smooth of degree 2: no sudden change in curvature but # possible sudden change 3rd derivative = jerk # - 3, 4, 5: snap, crackle and pop # # - Difficult: generating regressors that do the right job # BUT if the regressors are generated for you, nothing's # easier than using splines. # # gsp is a simple function to generate the columns of the X matrix # for arbitrary splines # # # Primer for splines # data(sun) # monthly sunspot data in 20th century spex <- sun[ sun$year > 1980,] xqplot( spex ) plot( spots ~ year , spex) # # Linear spline: lines continuous at knots: choose knots where slope changes: # # General Spline Generator: # Preliminary step: # - transform time so it has a resonable range # - max( abs(time) )^(highest degree) should not be much greater than 10000 # - you can keep the original time scale for plotting spex$ys <- spex$year - 1980 ?gsp # define a function to generate the spline: # linear spline splin <- function( x ) { gsp( x, knots = c( 1986, 1990, 1996 )-1980, # 3 knots => 4 intervals degree = c( 1, 1, 1, 1 ) , # linear in each interval smooth = c( 0, 0, 0 ) # continuous at each knot ) } fitlin <- lm( spots ~ splin(ys), spex ) summary( fitlin ) lines( predict(fitlin) ~ year, spex) abline( v = 1980 + splin(NULL)$knots , col = 'gray') # quadratic spline: quadratic in each interval and smooth at knots spquad <- function( x ) { gsp( x, knots = c( 1989, 1993 )-1980, # 2 knots => 3 intervals degree = c( 2, 2, 2 ) , # quadratic in each interval smooth = c( 1, 1 ) # smooth at each knot ) } fitquad <- lm( spots ~ spquad(ys), spex ) summary( fitquad ) lines( predict(fitquad) ~ year, spex, col = 'red') abline( v = 1980 + spquad(NULL)$knots , col = 'pink') # note the appearance of the spline where it goes over the knot # cubic spline spcubic <- function( x ) gsp( x, 1987 - 1980, c(3, 3), 2) fitcubic <- lm( spots ~ spcubic(ys), spex ) summary( fitcubic ) lines( predict(fitcubic) ~ year, spex, col = 'blue') abline( v = 1980 + spcubic(NULL)$knots , col = 'blue') AIC ( fitlin, fitquad, fitcubic ) # fit linear was best for 6 dfs # but we might have done better # by finding better places for knots # natural cubic spline # add 'boundary knots' at the extreme values of the data and # extrapolate with linear term # spnat <- function( x ) gsp( x, c(1980, 1987, 1998) - 1980, c(1, 3, 3, 1), c(1,2,1)) fitnat <- lm( spots ~ spnat(ys), spex ) summary( fitnat) lines( predict(fitnat) ~ year, spex, col = 'green') abline( v = 1980 + spnat(NULL)$knots , col = 'green') # # Summary: # The right spline depends on the phenomenon being modeled # # Cubic splines are very popular but for many situations they # don't bend fast enough. # # AIC ( fitlin, fitquad, fitcubic, fitnat ) # # Using splines for age # quantile(ind$age, seq(0,100,5)/100) ind$age.y <- ind$age / 12 # transform to age in years -- good to limit max(abs(x)) -- # when using raw polynomial and splines to control condition # number of X matrix # try a cubic fit.3 <- glmmPQL ( resp ~ time + Season + ( age.y + I(age.y^2) + I(age.y^3)) * Sex + height_for_age , data = ind, family = binomial, random = ~ 1 | id) summary( fit.3 ) wald( fit.3, 'Sex' ) # does not reveal much # Try a natural cubic spline qs <- quantile( ind$age.y, c(0,1,2,3,4,5,6)/6) qs # natural cubic spline puts knots at both ends with linear terms to extrapolate if necessary sp <- function( x ) gsp( x, qs, c(1,3,3,3,3,3,3,1), c(2,2,2,2,2,2,2)) # natural cubic spline using 5 interior # knots separating intervals # witn equal proportions of the data # and two boundary knots fitsp <- glmmPQL( resp ~ sp(age.y) * Sex + time + Season + height_for_age, ind, random = ~ 1 | id, family = binomial ) summary( fitsp ) pred <- expand.grid( age = 4:84, Sex = levels( ind$Sex), height_for_age = 0, Season = levels(ind$Season), time = 1:6) pred$age.y <- pred$age / 12 pred$resp.lo <- predict( fitsp, pred, level = 0) wald( fitsp, 'Sex') # No strong evidence for an overall effect of Sex # But what if we target the question we wanted to ask # Effect of Season: xyplot( resp.lo ~ Season, pred, type = 'l', subset = Sex == "Male" & time == 1 & age == 50) wald( fitsp, Ldiff( fitsp, "Season", reflevel = "Spring")) # Overall period effect: wald( fitsp, list( "period effects" = "time|Season") ) # QUESTION: What does this say? Recall that we are # controlling for age at the time of observation xyplot( resp.lo ~ age, pred, groups = Sex, type = 'l', subset = Season == "Summer" & time == 1, auto.key = list( columns = 2, lines = T, points =F)) # # Using Wald tests to explore splines # # We need an estimate of the difference between the two curves wald( fitsp ) # # We see that the columns of the model matrix consist of: # 1: a column of 1s for the intercept, # 2-7: a block of 6 columns for the spline, # 8: male indicator # 9: time # 10-12: a block of 3 columns for Season # 13: height_for_age # 14-19: another block of 6 columns for the interaction between # age and Sex # # To estimate the difference between the two sexes at each age # we need to generate the correct L matrix. This would seem formidable # without a detailed understanding of the parametrization of the # spline. # # However, there is a function to help with this. It generates the # portions of the L matrix for the blocks that model the spline. # # # To generate the portion of the matrix to evaluate the spline at # at age.y = 2, 4, and 6, use sc( sp , c(2,4,6) ) # # This is the portion of the L matrix that would estimate the # value of the response offset from the intercept. # # To evaluate the derivative (slope) of the spline at the same # points: sc( sp, c(2,4,6), D = 1) # first derivative (slope) of spline # the second, third and fourth derivatives: sc( sp, c(2,4,6), D = 2) sc( sp, c(2,4,6), D = 3) sc( sp, c(2,4,6), D = 4) # # To handle possible discontinuites at knots, a fourth argument # 'type' that can be set to 0, 1 or 2 signifies that the # derivative should be evaluated from the left, from the right # or that the saltus should be evaluated. # To create an L matrix to estimate Female response at # each month, given values of other predictors. The values # are eventually arbitrary and we choose: # # Season == "Spring" & time == 2 # # We need: # # Lfemale = the L matrix is Season == "Summer" & time == 1, # Lfemale <- cbind( 1, sc( sp, seq(4,84,1)/12), 0, 1, 1, 0, 0, 0, 0*sc(sp,seq(4,84,1)/12) ) some( Lfemale ) zapsmall( some( Lfemale )) Lmale <- cbind( 1, sc( sp, seq(4,84,1)/12), 1, 1, 1, 0, 0, 0, 1*sc(sp,seq(4,84,1)/12) ) zapsmall( some( Lmale )) Ldiff <- Lmale - Lfemale zapsmall( some( Ldiff) ) (ww <- wald( fitsp, Ldiff )) # estimates gap at each age wald( fitsp, "Sex") # QUESTION: How is this output similar to the previous wald output? # Turn the output of the 'wald' function into a data frame for plotting # include eta.hat +/- 2 SE df.diff <- as.data.frame( wald( fitsp, Ldiff), se = 2) df.diff$age <- seq(4,84,1) some( df.diff ) # estimated coef with Upper and Lower confidence Bound xyplot( coef + U2 + L2 ~ age, df.diff, type = 'l', panel = function(x,y,...) { panel.superpose(x,y,...) panel.abline( h = 0) }) # Better labels td( lty = c(1,2,2), col = c('blue','blue','blue'), lwd = 2) xyplot( coef + U2 + L2 ~ age, df.diff, type = 'l', sub = "Boy-Girl gap in log odds of respiratory infection", ylab = "Log odds difference +/- 2 SEs (single prior hypothesis)", xlab = "age (in months)", panel = function(x,y,...) { panel.superpose(x,y,...) panel.abline( h = 0) }) # To make this more interpretable, we would like to transform # graphs to a probability scale where possible. # Fitted responses are easy: probability = 1/(1+exp(-log.odds)) # Let's write it as a function Invlogit <- function ( lo ) 1/(1+exp(-lo)) Logit <- function( p ) log( p / (1-p) ) # i.e. log(odds(p)) # So the graph showing each sex is, instead of: xyplot(resp.lo ~ age, pred, groups = Sex, type = 'l', subset = Season == "Summer" & time == 1, auto.key = list( columns = 2, lines = T, points =F)) # we get: xyplot( Invlogit(resp.lo) ~ age, pred, groups = Sex, type = 'l', ylab = "Probability of respiratory infection", subset = Season == "Summer" & time == 1, auto.key = list( columns = 2, lines = T, points =F)) # QUESTION: Are the bumps and valleys artifacts of our knots? # # EXERCISE LATER: redo this playing with different knots, # possibly knots chosen on the basis of the initial plot # What price should you pay? The knot adjustments are # non-linear parameters and you could pay 1 DF per knot # to be safe. # # # Plotting the difference: # # There isn't a direct way to turn the difference in log odds into # a difference of probabilities. Probabilities are 1-1 functions of log odds # but probability differences are not 1-1 functions of differences in log odds # # However: # # the ratio of odds = exp( diff of log odds ) # and if the probabilities are relatively small, then the ratio of odds # is approximately the proportion of probabilities which we can # wrap our minds around. # # td( lty = c(1,2,2), col = c('blue','blue','blue'), lwd = 2) xyplot( exp(coef) + exp(U2) + exp(L2) ~ age, df.diff, type = 'l', sub = "Boy/Girl odds ratio for respiratory infection", ylab = "Odds ratio +/- 2 SEs (single prior hypothesis)", xlab = "age (in months)", panel = function(x,y,...) { panel.superpose(x,y,...) panel.abline( h = 0) }) # oops xyplot( exp(coef) + exp(U2) + exp(L2) ~ age, df.diff, type = 'l', ylim = c(.1,10), sub = "Boy/Girl odds ratio for respiratory infection", ylab = "Odds ratio +/- 2 SEs (single prior hypothesis)", xlab = "age (in months)", panel = function(x,y,...) { panel.superpose(x,y,...) panel.abline( h = 1) }) # not pretty. If the changing the plotted values does not work, then # just change the labels on the axes! xyplot( coef + U2 + L2 ~ age, df.diff, type = 'l', ylim = log( 2^c(-8,8)), scale = list( y = list( at = log( 2^seq(-7,7,1)), labels = fractions( 2^seq(-7,7,1))), x = list( at = seq( 0,90,10))), sub = "Boy/Girl odds ratio for respiratory infection (on log odds scale)", ylab = "Odds ratio +/- 2 SEs (single prior hypothesis)", xlab = "age (in months)", panel = function(x,y,...) { panel.superpose(x,y,...) panel.abline( h = log( 2^seq(-7,7,1)), col = 'gray') panel.abline( h = 0) panel.abline( v = seq(0,90,5), col = 'gray') }) # # So ... at 35 months it seems that boys are estimated to be over 4 times # as likely as girls to get respiratory infections and # they are at least twice as likely as girls using a 95% CI. # # BUT this strictly works only if we have one age as a prior hypothesis. # What if we have selected 35 months because it is significant. # # We can adjust the CIs to take multiple prior hypotheses into # account using Bonferroni adjustment or to protect against any # number by using a Scheffe adjustment. # # # Buying a fishing license # # # Rules for fishing licenses: # # The analogy doesn't quite work because # you are not paying for the number of fish you catch, # you are paying for the particular fish you plan to # look at whether you decide to keep them or not. # If you've paid to look at the fish, # there's no extra cost for keeping it! # # This subtlety is critical: it's your intentions that # count, not your actions!! You pay for the privilege # to look at fish, not for keeping them. # # 1) No window shopping. You have to pay for every fish you PLAN to look at. # # 2) The first fish is free. # # 3) Bonferroni's franchise: if you plan to look at a fixed number of fish # you get your licence from Mr. Bonferroni who will charge you per # fish you plan to look at. # # 4) Scheffe's franchise: You can get an unlimited fishing licence for a # SUBSPACE of hypotheses from Mr. Scheffe. He charges according to # the dimension of the subspace. # # 5) If you find that a Scheffe license is cheaper for your purposes # than a Bonferroni license, you can get a full refund on your # Bonferroni license and use a Scheffe license instead. # # Note: there are many other purveyors of fishing licences but # they are special purpose licences. They tend to be cheaper when # they are just right for you application. But they aren't as general # as Bonferroni and Scheffe. # # How do we buy fishing licences? # # It's extremely easy: # When we use # # estimated.value +/- 2 x SE # # the '2' comes from the fact that the interval -2 to +2 contains # approximately 95% of the probability for a standard normal random # variable. # # All we need to do is to use something other than a '2' # # For Bonferroni to test H hypotheses, we want to choose K (instead of 2) # so that the probability that H normals (independent or not) will all # fall in the interval (-K , K) is at least 95%. # # For Scheffe to fish in a SUBSPACE of dimension D, we want to choose K # so that any linear combination chosen from a particular subspace of # dimension D, will fall in the interval (-K , K) # # The following are two functions that will produce the right factor # to use instead of '2': Bonferroni.factor <- function( nhypotheses, alpha = .95, denDF = Inf) qt( 1 - (1-alpha)/(2*nhypotheses), denDF) Scheffe.factor <- function( dimension , alpha = .95, denDF = Inf) sqrt( dimension * qf( alpha, dimension, denDF)) # Note: qnorm( .975 ) Bonferroni.factor( 1 ) Scheffe.factor( 1 ) cat( "\n\nCost of protection:\n") mat <- cbind( Bonf = Bonferroni.factor( 1:20 )/qnorm(.975) , Sch = Scheffe.factor( 1:20)/qnorm(.975)) rownames( mat ) <- 1:20 round( mat, 2) # # Note that costs are not really comparable because you aren't buying the # same thing. # # # Some scenarios: # # Scenario 1: # # Suppose we had preselected (credibly on some basis independent of the data) # three ages we wanted to test, say 36,42 and 48 months. This could use # a Bonferroni adjustment with nhypotheses = 3. # # Scenario 2: # # Suppose we want to test 36 months but we played with two of the knots to # make the gap as big as possible. We could use a Scheffe adjustment # with dimension 3 (2 more than 1) # # Scenario 3: # # Suppose we're just looking for any Sex differences at all, we have # no preconceived idea where to look. We could use s Scheffe adjustment # with dimension equal to the dimension of the Sex comparisons. This is # numDF in the wald output for Ldiff wald(fitsp, Ldiff) # i.e. 7. # # QUESTION: How many prior hypotheses do you think we would need to test # before a Scheffe adjustment is cheaper than a Bonferroni adjustment. # # Last question first: uniroot( function(K) Bonferroni.factor(K) - Scheffe.factor(7), c(1,1000)) # Scheffe is, surprisingly, quite expensive. # # Now the first 3 questions: ses = c(std = qnorm(.975), Bonf.3 = Bonferroni.factor( 3 ), Sch.3 = Scheffe.factor( 3 ), Sch.7 = Scheffe.factor( 7 )) ses dmult <- as.data.frame( wald( fitsp, Ldiff), se = ses) some(dmult) dmult$age <- seq(4,84,1) td( col = c('blue','blue','green3','red','black','blue','green3','red','black'), lty = c( 1,2,3,4,5,2,3,4,5), lwd = 2) xyplot( coef + Ustd + UBonf.3 + USch.3 + USch.7 + Lstd + LBonf.3 + LSch.3 + LSch.7 ~ age, dmult, type = 'l', ylim = log( 2^c(-8,8)), scale = list( y = list( at = log( 2^seq(-7,7,1)), labels = fractions( 2^seq(-7,7,1))), x = list( at = seq( 0,90,10))), sub = "Boy/Girl odds ratio for respiratory infection (on log odds scale)", ylab = "Odds ratio with confidence intervals", xlab = "age (in months)", panel = function(x,y,...) { panel.superpose(x,y,...) panel.abline( h = log( 2^seq(-7,7,1)), col = 'gray') panel.abline( h = 0) panel.abline( v = seq(0,90,5), col = 'gray') }, auto.key = list( text = c('Estimated','Ordinary 95% CI', "Bonferroni k=3", "Scheffe d=3", "Scheffe d=7"), lines = T, points = F, columns = 2) ) # # EXERCISE: # What happens if, instead of using Scheffe(7) we resort to # Bonferroni on a fine grid of ages, e.g. every whole month. # Do we get significant results? Is this a reasonable approach? # # QUESTION: # Is there a connection between the p-value of wald( fitsp, Ldiff)[[1]]$anova$`p-value` # of 0.06 with 7 dfs in the overall F-test for Sex differences and the fact that # the specific Sex differences at each age are never significant when # we use the Scheffe( 7 ) adjustment? # # EXERCISES: # 1) Experimenting with knots: # See what happens if you choose different knots. # Coordinate with participants near you to try # different experiments and then compare results. # Fit a model with different knots, then compare # results for testing sex differences: # Things to try: # 1) Try 7 interior knots instead of 5 # 2) Try 3 interior knots and move them around # to increase the significance (decrease the # p-value) for sex comparisons # 3) Try a natural quadratic spline with 5 # interior knots # # e.g. Try using 3 interior # possibly knots chosen on the basis of the initial plot # What price should you pay? The knot adjustments are # non-linear parameters and you could pay 1 DF per knot # to be safe. # # # Age - period = cohort # # We performed the analysis above ignoring cohort (baseline_age), i.e. # assuming there is no effect. # # Consider: # # If 20-year =-olds in 2000 are much healthier than 20-year-olds in 1990 # this could have two explanations: # - General health improved a lot from 1990 to 2000 (period effect), or # - People born in 1980 are much healthier than people born in # 1970 (cohort effect) # What if people of all ages saw steadily improving health from 1990 to # to 2000. Again this could be the period or it could be that # people born later are healthier than people born earlier. # Of course, if this is the explanation, we would expect the # the phenomenon to be persistent. If some cohorts improve steadily # from 1990 to 2000 and others do not, then the explanation that # it's only a period effect is harder to support. # # The problem is that steadily (linearly) improving health over a # long period of time looks exactly the same as lineary improving # cohorts over a long period of time. Age, period and cohort are # collinear. The only way to break the collinearity would be to find # a 20-year-old in 2000 who was born in 1970. # # So disentangling linear effects can't be done. However, non-linear # effects are a different thing. For example if everyone suddenly # becomes less stressed between 1994 and 1995, it could be a period # effect (the referendum?) or it could be that a characteristic of # people born in 1970 was that they would become less stressed at the age of 25, # people born in 1971 would become less stressed at the age of 24, # people born in 1972 would become less stressed at the age of 23, etc. # The only way to blame a sudden change in one calendar year without # accepting that it's a period effect is to postulate a very complex # and implausible interaction between cohort and age. # # In our analysis, we did not include cohort effects. If we had tried # by including baseline_age, we would have had a collinearity in the model. # Is there a way of verifying whether it's reasonable to ignore cohort # effect and to attribute changes to age and period? One way is to test # for cohort effects that are beyond linear. We can do this by # generating a spline and then removing its linear component. If # such cohort effects are negligible, we gain some confidence about # ignoring them. If they are significant, we will wish to explore # them. # # # Since our model is already burdened with parameter, we will simplify # it by removing all but time and age effects before adding higher-order # cohort effects. # # In addition, the fact that we can never disentangle linear age # effects from possibly compensating linear period and # linear cohort effects implies that the non-linear effects # of age should be particularly interesting since they # are the only age effects that are identifiable. # # summary( ind$baseline_age) ind$cohort <- scale( ind$baseline_age ) # mean 0 sd 1 quantile( ind$cohort ) cqs <- quantile( ind$cohort, (1:5)/6) # five knots some(gsp( ind$cohort, cqs, c(2,2,2,2,2,2), c(1,1,1,1,1))) spcohort <- function( x ) gsp( x, cqs, 2, 1)[,-1] # remove 1st column fitspc <- glmmPQL( resp ~ sp(age.y) + time + Season + spcohort( cohort) , ind , random = ~ 1 | id, family = binomial ) summary( fitspc ) wald( fitspc, "cohort" ) # there is no striking evidence of a cohort effect suggesting that # the earlier analysis is reasonable # Recall that our effective sample size is somewhat small so this # result needs to be interpreted cautiously # # Alternatives to glmmPQL # # Advantages of glmmPQL: # # -- Based on lme and can model G and R side models like lme, # -- Can be used for smoothing splines with mixed models # -- Its limitations are somewhat understood. # -- Genearally fast # -- Wide set of methods for lme can be applied # # Disadvantages: # # -- Questionable estimates of components of variance # especially with binary # -- not based on likelihood - can't do LRTs # # glmmPQL has flaws in theory but is still a popular workhorse # because it generally converges and is easy to use # # Alternatives for generalized linear mixed models: # # glmmML in glmmML: # likelihood based, random intercept only. # # lmer in package lme4: # likelihood based, handles crossed effects well (so does lme # but not quite as well), very advanced numerical underpinnings, # no R side, more limited options for G side # hence no smoothing splines. # Many methods have not been developed yet. # lme4 has a function to perform a MCMC analysis on a # fitted model but there are bugs with the method # that need to be resolved. # # glmm in GLMMGibbs: beta release, limited choices: binomial and poisson # Use Bayesian inference. Can use glmmPQL for starting values. # # # glmmML # library(glmmML) fit.ML <- glmmML( resp ~ sp(age.y) * Sex + time + Season + height_for_age, data = ind, cluster = id, family = binomial ) summary( fit.ML ) coef( fit.ML ) str( fit.ML ) str(fitsp) sqrt( diag( fitsp$varFix)) # comparison with glmmPQL comparison <- data.frame( "PQL_coef" = pc <- fixef( fitsp ) , "PQL_sd" = ps <- sqrt( diag( fitsp$varFix)), "PQL_z" = pc/ps, "ML_coef" = mc <- coef(fit.ML), "ML_sd" = ms <- fit.ML$coef.sd, "ML_z" = mc/ms) comparison # ML is more conservative. I would expect them to be more similar # with a larger effective sample size and more observations per cluster. # # Some disadvantage: # Limited G side, no R side, only two families, limited links, # Few methods: (e.g. predict, resid, vcov do no exist yet) but it would # probably be easy to write them # # # lmer # # recall: fitsp <- glmmPQL( resp ~ sp(age.y) * Sex + time + Season + height_for_age, ind, random = ~ 1 | id, family = binomial ) library( lme4 ) # this used to conflict with nlme and one would have # to close and restart R before loading the other # I have experienced some problems but the two # should work much better now fitlmer <- lmer( resp ~ sp(age.y) * Sex + time + Season + height_for_age + (1|id), data = ind, family = binomial) ss <- summary(fitlmer) fixef( fitlmer) # coefficient sqrt( diag( vcov ( fitlmer ))) # their sds # EXERCISE: # # Add lmer to the comparison matrix. # Does lmer look closer to glmmPQL or to glmmML? # Plot the profile of coefs, sds and zs superposed for each method # What patterns emerge? # # # # Future projects # # 1) Add robust covariance option to 'wald' 2) Add influence diagnostics to nlme 3) Wald test for non-linear functions of parameters 4) Bootstrapping for mixed models 5) Weighted bootstrapping for mixed models 6) Clean up documentation for 'spida' and 'p3d'