# # Lab3 # Non-linear models # library( spida2 ) library(nlme) library(lattice) library(latticeExtra) library(car) # # Polynomials and interepretation of parameters # x = seq(-3, 2, .01) y = 1 + 2*x - .5 * x^2 - .4 * x^3 plot( x, y , type = 'l', lwd = 2, asp = 1) abline ( h =0) abline ( v =0) lines( x , 1 + 2*x, col = 'blue', lwd = 2) lines( x , 1 + 2*x + .5 * x^2, col = 'red', lwd= 2) # # IQ recovery data # iq <- spida2::iq ?iq data ( iq ) ls ( all = T ) head( iq ) names(iq) <- tolower( names( iq )) # make sure no conflicts head( iq ) xqplot( iq ) xqplot( up( iq, ~ id )) xyplot( viq + piq ~ dayspc, iq, groups = id, type = 'b', outer = T) xyplot( viq + piq ~ dayspc, iq, groups = id, type = 'b', outer = T, layout = c(1,2), subset = dayspc < 3000) # # small simulation for a single subject # plot(0,0, xlim = c(0,800), ylim = c(85,105)) iqsim.ex <- locator( 10 , type = 'p') iqsim.ex iqsim.ex <- as.data.frame( iqsim.ex ) iqsim.ex names( iqsim.ex ) <- c('days','iq') iqsim.ex <- iqsim.ex[ order(iqsim.ex$days),] iqsim.ex # However we will use a precooked version: data( iqsim ) # in spida2 iqsim plot( iq ~ days , iqsim, pch = 16, xlim = c(0,800), ylim = c(85,105) ) fit.lin <- lm ( iq ~ days, iqsim ) summary( fit.lin ) # We would like to show predicted values as a function of days # over the whole range in the plot pred <- expand.grid( days = seq( -20, 850, 1)) some(pred) pred$iq.lin <- predict( fit.lin, pred ) some( pred ) lines( iq.lin ~ days, pred, col = 'blue', lwd = 2) # Not very realistic. Need something better. Try a quadratic: fit.quad <- lm( iq ~ days + I(days ^2), iqsim) summary( fit.quad ) pred$iq.quad <- predict( fit.quad, pred ) lines( iq.quad ~ days , pred, col = 'red', lwd = 2) # Not very realistic either fit.cube <- lm( iq ~ days + I(days ^2) + I( days^ 3), iqsim) summary( fit.cube ) pred$iq.cube <- predict( fit.cube, pred ) lines( iq.cube ~ days , pred, col = 'purple', lwd = 2) p8 <- function( x ) poly( x, 8, raw = TRUE) fit.high <- lm( iq ~ p8( days ), iqsim ) summary(fit.high) # look at R-Squared!! pred$iq.high <- predict( fit.high, pred) lines( iq.high ~ days , pred, col = 'magenta', lwd = 2) # Beyond a point, the better you fit the data, the worse you fit the population or process # Try a non-linear (in parameters) asymptotic model fit.nl <-nls ( iq ~ b0 + b1*exp( -a*days ), iqsim, start= list(b0 = 100, b1 = -30, a = .01)) # where do we get the starting values?? # Draw a sketch and understand the parameters! # We'll see how soon. summary( fit.nl ) pred$iq.nl <- predict( fit.nl, pred ) plot( iq ~ days , iqsim, pch = 16, xlim = c(0,800), ylim = c(85,105) ) lines( iq.nl ~ days , pred, col = 'blue', lwd = 2) coef( fit.nl ) abline( h = coef(fit.nl)[1], col = 'gray', lwd = 2) # asymptotic level # # Turning a non-linear model into a linear model by transformation # ttime <- function( x ) exp( -0.0058 * x) fit.lin <- lm( iq ~ ttime( days ), iqsim) summary(fit.lin) # # Fitting non-linear model # library(nlme) fit.nlme <- nlme( piq ~ b0 + b1*exp(-a*dayspc), data = iq, fixed = list( b0 ~ 1 + sqrt(dcoma) , b1 ~ 1 + sqrt(dcoma) , a ~ 1), random = list( id = list( b0 ~ 1, b1~ 1 )), start = list( fixed = c(100, -10, -10, 0,.05)), control = list( maxIter = 100, msMaxIter = 100, returnObject = T), verbose = TRUE) summary( fit.nlme ) fit.nlme <- nlme( piq ~ b0 + b1*exp(-a*(dayspc - 30)), data = iq, fixed = list( b0 ~ 1 + sqrt(dcoma) , b1 ~ 1 + sqrt(dcoma) , a ~ 1), random = list( id = list( b0 ~ 1, b1~ 1 )), start = list( fixed = c(100, -10, -10, 0,.05)), control = list( maxIter = 100, returnObject = T), verbose = TRUE) summary( fit.nlme ) margSD <- sqrt( 13.769323^2 + 6.736146^2) # this is 15 for IQ tests margSD TRreliab <- 13.769293^2 / margSD^2 # this is thought to be 0.80 for IQ tests TRreliab plot( fit.nlme, resid(. , type = 'p') ~ fitted(.), id = .05) plot( fit.nlme, sqrt( abs( resid(. ,type='p'))) ~fitted(.), id = .05) plot( ranef( fit.nlme )) pairs( ranef( fit.nlme )) # # The VIQ story # fit.nlme.viq <- nlme( viq ~ b0 + b1*exp(-a*dayspc), data = iq, fixed = list( b0 ~ 1 + sqrt(dcoma) , b1 ~ 1 + sqrt(dcoma) , a ~ 1), random = list( id = list( b0 ~ 1, b1~ 1 )), start = list( fixed = c(100, -.3, -10, -5,.01)), control = list( maxIter = 100, returnObject = T), verbose = T) # Unset buffered output to see what's happening fit.nlme.viq.1 <- nlme( viq ~ b0 + b1*exp(-a*dayspc), data = iq, fixed = list( b0 ~ 1 + sqrt(dcoma) , b1 ~ 1 + sqrt(dcoma) , a ~ 1), random = list( id = list( b0 ~ 1, b1~ 1 )), start = list( fixed = c(100, -.3, -10, -5,.1)), control = list( maxIter = 400, returnObject = T), verbose = T) # Unset buffered output to see what's happening plot(resid(fit.nlme.viq), resid(fit.nlme.viq.1)) fit.nlme.viq2 <-nlme( viq ~ b0 + b1*exp(-a*(dayspc-30)), data = iq, fixed = list( b0 ~ 1 + sqrt(dcoma) , b1 ~ 1 + sqrt(dcoma) , a ~ 1), random = list( id = list( b0 ~ 1, b1~ 1 )), start = list( fixed = c(100, -.3, -10, 0,.3)), control = list( maxIter = 100, returnObject = T), verbose = T, subset = dcoma < 100) summary( fit.nlme.viq2 ) fit.nlme.piq2 <-nlme( piq ~ b0 + b1*exp(-a*(dayspc-30)), data = iq, fixed = list( b0 ~ 1 + sqrt(dcoma) , b1 ~ 1 + sqrt(dcoma) , a ~ 1), random = list( id = list( b0 ~ 1, b1~ 1 )), start = list( fixed = c(100, -.3, -10, 0,.1)), control = list( maxIter = 100, returnObject = T), verbose = T, subset = dcoma < 100) summary( fit.nlme.piq2 ) summary( fit.nlme.viq2 ) plot( ranef( fit.nlme.viq2)) pairs( ranef( fit.nlme.viq2)) plot( fit.nlme.viq2 ) summary( fit.nlme.piq2 ) plot( ranef( fit.nlme.piq2)) pairs( ranef( fit.nlme.piq2)) plot( fit.nlme.piq2 ) # # Visualizing the comparison # windows( height = 7, width = 8.5) pred <- expand.grid( dcoma = c(0,1,7,16,25,100), dayspc = seq(30,365*2,5)) pred$piq <- predict( fit.nlme.piq2, pred, level = 0 ) pred$viq <- predict( fit.nlme.viq2, pred, level = 0 ) zz <- factor( paste( 'dcoma =', pred$dcoma)) pred$dcoma.lab <- reorder( zz, pred$dcoma) td( col = c('blue','red'), lwd = 2) xyplot( viq + piq ~ dayspc | dcoma.lab, pred, type = 'l', ylim = c(60,102), lwd = 2, auto.key = list(columns = 2, points = F, lines = T)) pred2 <- expand.grid( dcoma = 0:100, dayspc = c(30,60,90.180,360,720)) pred2$piq <- predict( fit.nlme.piq2, pred2, level = 0 ) pred2$viq <- predict( fit.nlme.viq2, pred2, level = 0 ) zz <- factor( paste( 'dayspc =', pred2$dayspc)) pred2$dayspc.lab <- reorder( zz, pred2$dayspc) xyplot( viq + piq ~ dcoma | dayspc.lab, pred2, type = 'l', ylim = c(60,102), lwd = 2, auto.key = list(columns = 2, points = F, lines = T)) predviq <- expand.grid( dcoma = seq(0,100,10), dayspc = seq(30,720,30)) predpiq <- predviq predviq $ iq <- predict( fit.nlme.viq2, predviq, level = 0) predpiq $ iq <- predict( fit.nlme.piq2 , predpiq, level = 0) predpiq$type <- factor( "PIQ" ) predviq$type <- factor( "VIQ" ) wireframe( iq ~ dayspc + dcoma | type, Rbind( predpiq, predviq)) wireframe( iq ~ dcoma + dayspc | type, Rbind( predpiq, predviq), scales = list( arrows = F), col = 'blue') wireframe( iq ~ dcoma + dayspc | type, Rbind( predpiq, predviq), scales = list( arrows = F), col = 'blue') wireframe( iq ~ dcoma + dayspc | type, Rbind(predpiq, predviq), scales = list( arrows = F), col = 'blue', xlab = 'Coma', ylab = 'Days post coma', screen = list( z = -65, x = -75 )) library( p3d ) Init3d() Plot3d( piq ~ dayspc + dcoma , subset(iq, dayspc < 800 & dcoma < 100)) wireframe( iq ~ dcoma + dayspc , Rbind(predpiq, predviq), groups = type, scales = list( arrows = F), col = 'blue', xlab = 'Coma', ylab = 'Days post coma', screen = list( z = -65, x = -75 ), auto.key = list(columns=2, lines = T, points = F), alpha = .5) #################################################### #################################################### # # Multivariate test of equivalence of IQ types # head(iq) iq <- within(iq, { iq__PIQ <- piq iq__VIQ <- viq occ <- 1:nrow(iq) # occasion ID }) head(iq) # # stacked multivariate model # iql <- tolong(iq, timevar = 'type', idvar = 'occ', sep = '__') iql$type <- factor(iql$type) # necessary for 'nlme' head(iql) fit.mult <- nlme( iq ~ b0 + b1*exp(-a*(dayspc-30)), data = iql, fixed = list( b0 ~ type/(1 + sqrt(dcoma)) - 1, b1 ~ type - 1 , a ~ type - 1), random = list( id = list( b0 + b1 ~ 1 ), occ = list(b0 ~ 1)), weights = varIdent( form = ~ 1 | type ), control = list( maxIter = 200, returnObject = T, msVerbose = T), start = list( fixed = c(100,-1,100,-1, -10,-10, .01,.01)) ) summary( fit.mult) L <- rbind( asymp = c(1,-1, 0,0, 0,0, 0,0), dcoma = c(0,0, 1,-1, 0,0,0,0), b1 = c(0,0, 0,0, 1,-1, 0,0), curvature = c(0,0, 0,0, 0,0, 1,-1)) wald(fit.mult, list('PIQ - VIQ' = L))