source("MAMMA.R") ################################################ # Percentage Bodyfat Data # ################################################ # Processing data # Read data fat.NA <- read.table(file = "percentage_body_fat.csv", header = TRUE, sep = ",", na.strings = ".") # > str(fat.NA) # 'data.frame': 1718 obs. of 5 variables: # $ subject : int 101 101 101 101 101 101 101 101 101 101 ... # $ visitno : int 1 1 2 2 3 3 4 4 5 5 ... # $ age : num 12.3 12.3 12.9 12.9 13.3 ... # $ instrument: Factor w/ 2 levels "caliper","DEXA": 1 2 1 2 1 2 1 2 1 2 ... # $ bodyfat : num 20.8 NA 21.7 17.5 23.2 ... # > # > colnames(fat.NA) # [1] "subject" "visitno" "age" "instrument" "bodyfat" # > colnames(fat.NA) <- c("subject", "rep", "time", "method", "meas") fat.NA$subject <- as.factor(fat.NA$subject) fat.NA$rep <- as.factor(fat.NA$rep) # order the data by subject, then by method, and then by rep # (this is ordering is assumed for creating design and covariance matrices) fat.NA <- with(fat.NA, fat.NA[order(subject, method, rep), ]) # Omit the missing observations fat <- na.omit(fat.NA) # > str(fat) # 'data.frame': 1515 obs. of 5 variables: # $ subject: Factor w/ 112 levels "101","102","103",..: 1 1 1 1 1 1 1 1 1 1 ... # $ rep : Factor w/ 9 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 2 ... # $ time : num 12.3 12.9 13.3 13.8 14.3 ... # $ method : Factor w/ 2 levels "caliper","DEXA": 1 1 1 1 1 1 1 1 1 2 ... # $ meas : num 20.8 21.7 23.2 26.3 23.9 ... # - attr(*, "na.action")=Class 'omit' Named int [1:203] 10 28 46 51 64 66 82 87 90 100 ... # .. ..- attr(*, "names")= chr [1:203] "2" "20" "38" "48" ... # > # Note: removal of missing observations breaks the pairing (i.e., one instrument # may be present, the other may be absent) # > table(fat$method) # caliper DEXA # 858 657 # > ############################# # Reshape the data frame so that caliper and DEXA measurements appear side-by-side caliper <- fat.NA[fat.NA$method == "caliper", ] DEXA <- fat.NA[fat.NA$method == "DEXA", ] # > summary(caliper[, 1:3] == DEXA[, 1:3]) # subject rep time # Mode:logical Mode:logical Mode:logical # TRUE:859 TRUE:858 TRUE:858 # NA's:0 NA's:1 NA's:1 # > fat.wide <- data.frame(subject = as.factor(caliper[, 1]), rep = as.factor(caliper[, 2]), time = caliper[, 3], caliper = caliper[, 5], DEXA = DEXA[, 5], diff = DEXA[, 5] - caliper[, 5]) # > str(fat.wide) # 'data.frame': 859 obs. of 6 variables: # $ subject: Factor w/ 112 levels "101","102","103",..: 1 1 1 1 1 1 1 1 1 2 ... # $ rep : Factor w/ 9 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 1 ... # $ time : num 12.3 12.9 13.3 13.8 14.3 ... # $ caliper: num 20.8 21.7 23.2 26.3 23.9 ... # $ DEXA : num NA 17.5 18.8 19.9 18.9 ... # $ diff : num NA -4.21 -4.4 -6.4 -4.97 ... # > fat.wide <- na.omit(fat.wide) # > str(fat.wide) # 'data.frame': 657 obs. of 6 variables: # $ subject: Factor w/ 112 levels "101","102","103",..: 1 1 1 1 1 1 1 1 2 2 ... # $ rep : Factor w/ 9 levels "1","2","3","4",..: 2 3 4 5 6 7 8 9 2 3 ... # $ time : num 12.9 13.3 13.8 14.3 14.8 ... # $ caliper: num 21.7 23.2 26.3 23.9 27.6 ... # $ DEXA : num 17.5 18.8 19.9 18.9 22.4 ... # $ diff : num -4.21 -4.4 -6.4 -4.97 -5.26 ... # - attr(*, "na.action")=Class 'omit' Named int [1:202] 1 10 19 24 28 30 37 42 45 46 ... # .. ..- attr(*, "names")= chr [1:202] "1" "10" "19" "24" ... # > # > length(unique(fat.wide$subject)) # [1] 91 # > # Note: The difference data has 91 subjects, whereas the original data has 112. # This happens because 21 subjects have DEXA measurements missing for all visits. # They can be examined using the following code: # > temp <- fat.NA[!(fat.NA$subject %in% fat.wide$subject), ] # > unique(temp$subject) # [1] 108 111 116 117 118 128 135 209 210 211 305 306 308 318 321 329 330 401 # [19] 404 419 422 # > ################################################ # Exploratory data analysis # ################################################ # Plot of individual trajectories xyplot(meas ~ time | method, group = subject, data = fat, type = c("g", "o"), col = "black", panel = function(...) { panel.superpose(..., lwd = 0.5, pch = 1, cex = 0.5, lty = 1) }, xlab = "age (years)", ylab = "percentage body fat") # Plot of trajectories of DEXA - caliper differences xyplot(diff ~ time, group = subject, data = fat.wide, type = c("g", "o"), col = "black", panel = function(...) { panel.superpose(..., lwd = 0.5, pch = 1, cex = 0.5, lty = 1) }, xlab = "age (years)", ylab = "difference (DEXA $-$ caliper)") # Scatterplot of caliper vs DEXA by visitno xyplot(DEXA ~ caliper | rep, groups = subject, data = fat.wide, type = c("g", "p"), col = "black", xlab = "percentage body fat, caliper", ylab = "percentage body fat, DEXA", prepanel = function(...) { list(xlim = c(12, 38), ylim = c(12, 38)) }, panel = function(...) { panel.xyplot(..., pch = 1, col.symbol = "black") panel.abline(a = 0, b = 1) }) # Bland-Altman plot by visitno xyplot(diff ~ I((DEXA + caliper)/2) | rep, groups = subject, data = fat.wide, type = c("g", "p"), xlab = "average percentage body fat", ylab = "difference (DEXA $-$ caliper)", panel = function(...) { panel.xyplot(..., pch = 1, col.symbol = "black") panel.abline(h = 0, col.line = "black") }) # Boxplot of difference by visitno bwplot(I(DEXA - caliper) ~ rep, data = fat.wide, xlab = "visit number", ylab = "difference (DEXA $-$ caliper)", pch = "|", panel = function(...) { panel.bwplot(...) panel.abline(h = 0, col.line = "grey") }) ######################################################################### # Modeling of paired data # ######################################################################### subid <- unique(fat$subject) nsub <- length(subid) nobs <- nrow(fat) # > nsub # [1] 112 # > nobs # [1] 1515 # > method.names <- levels(droplevels(fat$method)) rep.names <- levels(droplevels(fat$rep)) # > method.names # [1] "caliper" "DEXA" # > rep.names # [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" # > # the fixed-effect design matrix X <- model.matrix(~-1 + method + method:time + method:I(time^2) + method:I(time^3), data = fat) # create an alias for "method" for use in corCAR1 below instrument <- fat$method # model without subject-method interactions and continuous AR(1) errors fit.noint <- lme(meas ~ -1 + X, data = fat, method = "ML", random = list(subject = pdBlocked(list(~1, pdIdent(~-1 + rep)))), weights = varIdent(form = ~1 | instrument), correlation = corCAR1(form = ~time | subject/instrument)) # model with subject-method interactions and continuous AR(1) errors fit.int <- lme(meas ~ -1 + X, data = fat, method = "ML", random = list(subject = pdBlocked(list(~1, pdIdent(~-1 + method), pdIdent(~-1 + rep)))), weights = varIdent(form = ~1 | method), correlation = corCAR1(form = ~time | subject/instrument)) # compare models # > anova(fit.noint, fit.int) # Model df AIC BIC logLik Test L.Ratio p-value # fit.noint 1 13 6265.683 6334.884 -3119.842 # fit.int 2 14 6265.752 6340.276 -3118.876 1 vs 2 1.93149 0.1646 # > # We will use fit.noint (model without interaction) for all inferences from now on. ############################### # Plot fitted mean functions # ############################### fit <- fit.noint # get the fitted values # population level estimates: fitted(fit, level = 0) # individual level estimates: fitted(fit, level = 1) # prepare a data frame for plotting fat.new <- data.frame(fat, mean.hat = fitted(fit, level = 0)) fat.new <- fat.new[order(fat.new$time), ] # Plot of individual trajectories with fitted mean functions superimposed p1 <- xyplot(meas ~ time | method, group = subject, data = fat, type = c("g", "o"), col = "black", panel = function(...) { panel.superpose(..., lwd = 0.5, pch = 1, cex = 0.5, lty = 1) }, xlab = "age (years)", ylab = "percentage body fat") p2 <- xyplot(mean.hat ~ time | method, data = fat.new, type = c("g", "l"), col = "dark grey", lwd = 4) p1 + as.layer(p2) ############################### # Perform model diagnostics # ############################### # Note: fit represents model without interaction plot(fit, resid(., type = "n") ~ fitted(.), type = c("p", "smooth")) plot(fit, resid(., type = "n") ~ fitted(.) | method, type = c("p", "smooth")) plot(fit, resid(., type = "n") ~ time | method, type = c("p", "smooth")) plot(fit, abs(resid(., type = "n")) ~ fitted(.) | method, type = c("p", "smooth")) qqnorm(fit, ~resid(., type = "n") | method) # Semivariogram estimates of normalized residuals # model without interaction and with CAR errors vgram.car <- variogram.resid(fat, resid(fit, type = "n")) df.vgram.car <- data.frame(rbind(vgram.car[[1]], vgram.car[[2]]), method = rep(names(vgram.car), each = 20)) xyplot(variog ~ dist | method, data = df.vgram.car, type = c("g", "p", "smooth"), xlab = "distance", ylab = "semivariogram") # ACF of normalized residuals of model without interaction and with CAR errors acf.car <- acf.resid(fat, resid(fit, type = "n"), 8) # > acf.car # $caliper # lag n.used acf cval accept.null # 1 0 858 1.00000000 NA NA # 2 1 745 -0.04404446 0.07180752 TRUE # 3 2 641 -0.04100931 0.07741392 TRUE # 4 3 542 -0.05102849 0.08418766 TRUE # 5 4 448 -0.09403735 0.09259959 FALSE # 6 5 357 -0.02042377 0.10373229 TRUE # 7 6 267 0.07711813 0.11994785 TRUE # 8 7 177 -0.05785589 0.14731991 TRUE # 9 8 88 -0.08350459 0.20893286 TRUE # $DEXA # lag n.used acf cval accept.null # 1 0 657 1.000000000 NA NA # 2 1 513 -0.041698516 0.08653452 TRUE # 3 2 435 0.036444824 0.09397308 TRUE # 4 3 355 0.076437355 0.10402409 TRUE # 5 4 294 -0.066793005 0.11430742 TRUE # 6 5 239 0.052363320 0.12677953 TRUE # 7 6 159 -0.007431597 0.15543525 TRUE # 8 7 79 -0.048659177 0.22051318 TRUE # 9 8 0 NaN Inf NA # > # plot the acf of model without interaction and with CAR errors acf.car <- data.frame(rbind(acf.car[[1]][-1, ], acf.car[[2]][-1, ]), method = rep(names(acf.car), each = 8)) p5 <- xyplot(acf ~ lag | method, data = acf.car, type = c("g", "h"), ylim = c(-0.23, 0.23), col.line = "black", lwd = 2, xlab = "lag", ylab = "autocorrelation", panel = function(...) { panel.xyplot(...) panel.abline(h = 0, col.line = "dark gray") }) p6 <- xyplot(cval ~ lag | method, data = acf.car, type = "l", ylim = c(-0.23, 0.23), lty = 2, col.line = "black") p7 <- xyplot(-cval ~ lag | method, data = acf.car, type = "l", ylim = c(-0.23, 0.23), lty = 2, col.line = "black") p5 + as.layer(p6) + as.layer(p7) # model with interaction and independent errors fit.indep <- lme(meas ~ -1 + X, data = fat, method = "ML", random = list(subject = pdBlocked(list(~1, pdIdent(~-1 + method), pdIdent(~-1 + rep)))), weights = varIdent(form = ~1 | method)) vgram.indep <- variogram.resid(fat, resid(fit.indep, type = "pear")) df.vgram.indep <- data.frame(rbind(vgram.indep[[1]], vgram.indep[[2]]), method = rep(names(vgram.indep), each = 20)) xyplot(variog ~ dist | method, data = df.vgram.indep, type = c("g", "p", "smooth"), xlab = "distance (years)", ylab = "semivariogram") # ACF of standardized residuals of model with interaction and independent errors acf.indep <- acf.resid(fat, resid(fit.indep, type = "n"), 8) # > acf.indep # $caliper # lag n.used acf cval accept.null # 1 0 858 1.00000000 NA NA # 2 1 745 0.22558298 0.07180752 FALSE # 3 2 641 -0.09488933 0.07741392 FALSE # 4 3 542 -0.20672852 0.08418766 FALSE # 5 4 448 -0.29978446 0.09259959 FALSE # 6 5 357 -0.26476776 0.10373229 FALSE # 7 6 267 -0.18529736 0.11994785 FALSE # 8 7 177 -0.43435552 0.14731991 FALSE # 9 8 88 -0.49838080 0.20893286 FALSE # # $DEXA # lag n.used acf cval accept.null # 1 0 657 1.0000000 NA NA # 2 1 513 0.2181688 0.08653452 FALSE # 3 2 435 -0.1116702 0.09397308 FALSE # 4 3 355 -0.1996512 0.10402409 FALSE # 5 4 294 -0.3877363 0.11430742 FALSE # 6 5 239 -0.3471122 0.12677953 FALSE # 7 6 159 -0.3789045 0.15543525 FALSE # 8 7 79 -0.4049424 0.22051318 FALSE # 9 8 0 NaN Inf NA # > ########################################### # Inference based on fitted model # ########################################### # set up a grid ngrid <- 25 x.g <- with(fat, seq(f = min(time), t = max(time), l = ngrid)) X.g <- cbind(rep(1, ngrid), x.g, x.g^2, x.g^3) # ngrid x 3 matrix colnames(X.g) <- c("intercept", "x", "x^2", "x^3") # estimates prob <- 0.9 estimates.noint <- get.param.estimates(fit.noint, interaction = F, method.names, rep.names, X.g, colnames(X), prob, err.cor = T) estimates.int <- get.param.estimates(fit.int, interaction = T, method.names, rep.names, X.g, colnames(X), prob, err.cor = T) # Verify computation of likelihood function # random effects design matrix # model without subject x method interaction # nobs x nsub*(1 + # levels of rep) block diagonal matrix Z.noint <- model.matrix(~1 + rep, data = fat, contrasts = list(rep = contrasts(fat$rep, F))) Z.noint <- lapply(split(data.frame(Z.noint), fat$subject), as.matrix) Z.noint <- as.matrix(bdiag(Z.noint)) # Note: After the splitting, the lists are ordered according to levels of # subject factor. In other words, the i-th list in Z.noint corresponds to level i # of the subject factor. When the lists are converted into a block diagonal # matrix, the block elements may not be in the right order, which is the order # in rows in Z.noint appear. But this works here, because # the original dat (fat) is by default ordered according to subject levels. # > nsub*(length(rep.names) + 1) # [1] 1120 # > dim(Z.noint) # [1] 1515 1120 # > # > loglik.fun(estimates.noint$param, fat, err.cor = T, Zmat = Z.noint, interaction = F) # [1] -3119.842 # > # > logLik(fit.noint) # 'log Lik.' -3119.842 (df=13) # > # > loglik.fun(estimates.noint$param, fat, err.cor = T, Zmat = Z.noint, interaction = F) - logLik(fit.noint) # 'log Lik.' -1.364242e-12 (df=13) # > # model with subject x method interaction # nobs x nsub*(3 + # levels of rep) block diagonal matrix Z.int <- model.matrix(~1 + method + rep, data = fat, contrasts = list(method = contrasts(fat$method, F), rep = contrasts(fat$rep, F))) Z.int <- lapply(split(data.frame(Z.int), fat$subject), as.matrix) Z.int <- as.matrix(bdiag(Z.int)) # > nsub*(length(rep.names) + 3) # [1] 1344 # > # > dim(Z.int) # [1] 1515 1344 # > # > loglik.fun(estimates.int$param, fat, err.cor = T, Zmat = Z.int, interaction = T) # [1] -3118.876 # > # > logLik(fit.int) # 'log Lik.' -3118.876 (df=14) # > # Get Hessian matrix # without interaction # hess.noint <- -hessian(loglik.fun, estimates.noint$param, dat = fat, Xmat = X, # Zmat = Z.noint, err.cor = T, interaction = F) # # Takes about 4 hours to run on my laptop # > hess.noint # V1 V2 V3 V4 V5 # 1 2.307108e+01 -1.768597e+01 3.187123e+02 -2.505591e+02 4.466813e+03 # 2 -1.768597e+01 2.234032e+01 -2.465195e+02 3.162409e+02 -3.484842e+03 # 3 3.187123e+02 -2.465195e+02 4.553616e+03 -3.503179e+03 6.589254e+04 # 4 -2.505591e+02 3.162409e+02 -3.503179e+03 4.616707e+03 -4.967874e+04 # 5 4.466813e+03 -3.484842e+03 6.589254e+04 -4.967874e+04 9.815519e+05 # 6 -3.589036e+03 4.525824e+03 -5.033593e+04 6.803357e+04 -7.161221e+05 # 7 6.348714e+04 -4.993260e+04 9.648916e+05 -7.141245e+05 1.475223e+07 # 8 -5.196269e+04 6.546261e+04 -7.310529e+05 1.011444e+06 -1.043452e+07 # 9 -4.579723e-02 2.487320e-02 -9.526456e-01 -3.421144e-01 -1.647767e+01 # 10 7.368869e-01 -7.578061e-01 1.001661e+01 -5.885880e+00 1.454654e+02 # 11 -8.039309e-01 8.133387e-01 -1.232620e+01 4.237218e+00 -1.980805e+02 # 12 1.128412e-01 -8.040577e-02 3.262239e+00 1.990776e+00 6.909284e+01 # 13 6.023052e-01 -6.883020e-01 6.047141e+00 -4.905460e+00 5.007481e+01 # V6 V7 V8 V9 V10 # 1 -3.589036e+03 6.348714e+04 -5.196269e+04 -0.04579723 0.7368869 # 2 4.525824e+03 -4.993260e+04 6.546261e+04 0.02487320 -0.7578061 # 3 -5.033593e+04 9.648916e+05 -7.310529e+05 -0.95264555 10.0166107 # 4 6.803357e+04 -7.141245e+05 1.011444e+06 -0.34211437 -5.8858800 # 5 -7.161221e+05 1.475223e+07 -1.043452e+07 -16.47766683 145.4653711 # 6 1.029617e+06 -1.032816e+07 1.567942e+07 -13.87768454 -21.9019063 # 7 -1.032816e+07 2.267887e+08 -1.509941e+08 -264.63873928 2236.6246278 # 8 1.567942e+07 -1.509941e+08 2.438414e+08 -317.25658118 496.1239888 # 9 -1.387768e+01 -2.646387e+02 -3.172566e+02 38.94543110 0.4446325 # 10 -2.190191e+01 2.236625e+03 4.961240e+02 0.44463250 52.2544698 # 11 -3.474042e+01 -3.276857e+03 -1.765092e+03 -1.19237285 41.6655508 # 12 7.052001e+01 1.304871e+03 1.586225e+03 9.05723697 29.0055514 # 13 -6.790625e+00 1.721337e+02 7.432798e+02 0.74258617 -43.3277252 # V11 V12 V13 # 1 -0.8039309 0.11284124 0.6023052 # 2 0.8133387 -0.08040577 -0.6883020 # 3 -12.3262038 3.26223866 6.0471410 # 4 4.2372179 1.99077649 -4.9054600 # 5 -198.0805423 69.09283805 50.0748080 # 6 -34.7404215 70.52001231 -6.7906254 # 7 -3276.8568182 1304.87092977 172.1336508 # 8 -1765.0924983 1586.22509079 743.2798101 # 9 -1.1923729 9.05723697 0.7425862 # 10 41.6655508 29.00555136 -43.3277252 # 11 272.7899537 24.93281641 -147.0492705 # 12 24.9328164 185.68331497 -96.3282237 # 13 -147.0492705 -96.32822367 174.3445383 # > # # > sqrt(diag(hess.noint)) # [1] 4.803236 4.726555 67.480485 67.946357 990.733001 # [6] 1014.700389 15059.503919 15615.420407 6.240627 7.228725 # [11] 16.516354 13.626567 13.203959 # > # get the **individual** confidence intervals and bands ci.noint <- conf.measures.longitudinal(estimates.noint$param, hess.noint, interaction = F, prob = prob) # confidence intervals for model parameters # > round(ci.noint$param, 2) # estimate se lcl ucl # methodcaliper -319.10 62.98 -442.54 -195.66 # methodDEXA 473.89 82.87 311.47 636.32 # methodcaliper:time 74.81 13.75 47.87 101.75 # methodDEXA:time -95.50 17.67 -130.14 -60.86 # methodcaliper:I(time^2) -5.40 0.99 -7.34 -3.45 # methodDEXA:I(time^2) 6.64 1.25 4.19 9.09 # methodcaliper:I(time^3) 0.13 0.02 0.08 0.18 # methodDEXA:I(time^3) -0.15 0.03 -0.21 -0.09 # lsigmasq.b 2.24 0.16 1.92 2.56 # lsigmasq.bstar -0.75 0.16 -1.06 -0.44 # lsigmasq.ecaliper 1.70 0.09 1.52 1.88 # lsigmasq.eDEXA 1.44 0.10 1.25 1.63 # logitphi 0.17 0.14 -0.09 0.44 # > # confidence bands for parameters that are functions of time ci.df <- data.frame(x.g = x.g, mean1.est = ci.noint$mean.1$result[, "estimate"], mean1.lcl = ci.noint$mean.1$result[, "lcl"], mean1.ucl = ci.noint$mean.1$result[, "ucl"], mean2.est = ci.noint$mean.2$result[, "estimate"], mean2.lcl = ci.noint$mean.2$result[, "lcl"], mean2.ucl = ci.noint$mean.2$result[, "ucl"], meanD.est = ci.noint$mean.D$result[, "estimate"], meanD.lcl = ci.noint$mean.D$result[, "lcl"], meanD.ucl = ci.noint$mean.D$result[, "ucl"], ccc.est = ci.noint$ccc$result[, "estimate"], ccc.lcl = ci.noint$ccc$result[, "lcl"], tdi.est = ci.noint$tdi$result[, "estimate"], tdi.ucl = ci.noint$tdi$result[, "ucl"], lloa = estimates.noint$loa[, "lower"], uloa = estimates.noint$loa[, "upper"]) # plot the mean difference bands and 95% LOA plot(ci.df$x.g, ci.df$meanD.est, ylim = c(-11, 11), type = "l", panel.first = grid(lty = "solid"), xlab = "age (years)", ylab = "difference (DEXA $-$ caliper)") polygon(c(ci.df$x.g, rev(ci.df$x.g)), c(ci.df$meanD.lcl, rev(ci.df$meanD.ucl)), col = "grey80", border = F) lines(ci.df$x.g, ci.df$meanD.est) lines(ci.df$x.g, ci.df$lloa, lty = 2) lines(ci.df$x.g, ci.df$uloa, lty = 2) # plot CCC and TDI bands and their 95% individual confidence bands p8 <- xyplot(ccc.est ~ x.g | rep("CCC", length(x.g)), data = ci.df, ylim = c(0.14, 0.7), ylab = NULL, xlab = "age (years)", panel = function() { panel.xyplot(ci.df$x.g, ci.df$ccc.est, type = c("l", "g"), col = "black") panel.polygon(c(ci.df$x.g, rev(ci.df$x.g)), c(ci.df$ccc.est, rev(ci.df$ccc.lcl)), border = F, col = "grey80") }) p9 <- xyplot(tdi.est ~ x.g | rep("TDI", length(x.g)), data = ci.df, ylim = c(4.5, 13), ylab = NULL, xlab = "age (years)", panel = function() { panel.xyplot(ci.df$x.g, ci.df$tdi.est, type = c("l", "g"), col = "black") panel.polygon(c(ci.df$x.g, rev(ci.df$x.g)), c(ci.df$tdi.est, rev(ci.df$tdi.ucl)), border = F, col = "grey80") }) print(p8, split = c(1, 1, 2, 1), more = T) print(p9, split = c(2, 1, 2, 1)) #######################################################################################