Macro name: %spmm ~~~~~~~~ SEMI-PARAMETRIC MIXED MODELS ~~~~~~~~ Version 1.0 Function: This macro fits the following model to longitudinal Gaussian data: Yij = Xij*beta + f(tij) + Zij*bi + Ui(tij) + esp(ij), where beta is parametric fixed effects, f(.) is a smooth function, bi is random effects, Ui(.) is a Gaussian process, esp(ij) is the measurement error. Maximum penalized likelihood was used to estimate the beta and f(.), while smoohting parameter and the parameters in the variance matrix are estimated by REML method, which treats f(.) as an integrated Wiener process. Macro developer: Daowen Zhang Department of Biostatistics The University of Michigan Ann Arbor, MI 48109 E-mail: dzhang@umich.edu (C) Daowen Zhang Date: April 17, 1996 Call statement: Default %spmm(data=, dep=, fixed=, smthvar=, smooth=, method=, lintest=, random=, process=, time=, gt=, g1t=, g2t=, id=, msrerr=, {Y} maxiter=, {50} conv=, {0.00001} print= {Y} outbeta= outbstd= outvar= outvstd= outran_F= outran_B= outprs= outsmth= outband= outscore= symsize= {1022976} worksize= {1022976}) where data = name of the data set to be used; dep = dependent variable (Y); fixed = fixed effects covariates (X); NO INTERCEPT ! /s -- print the solution; smthvar = covariate that needs smoothed (t); if missing, a linear mixed model is fitted. smooth = 1/smoothing parameter, can be missing; lintest = Y/N -- request the performance of linearity test based on the score test statistics for f(t) or the first covariate in X if smthvar is missing; The default is N; method = ML/REML -- Method to be used to fit the model; Default is REML; ML is only effective if smthvar is missing; random = random effects covariates (Z); The random effects covariance matrix is assumed to be unstructured; /s -- print the solution; process = name of the process: AR --- AR(1); OU --- Ornstein-Uhlenbeck process, NOU1 --- Nonhomogeneous OU process with log-variance function: log(v(t)) = a0 + a*g(t), where g(t) is specified in gt, NOU2 --- Nonhomogeneous OU process with log-variance function: log(v(t)) = a0 + a1*exp(alpha*t), NOU3 --- Nonhomogeneous OU process with log-variance function: log(v(t)) = a0 + a1*t + a2*t^2, NOU4 --- Nonhomogeneous OU process with log-variance function: log(v(t)) = a0 + a1*g1(t) + a2*g2(t), NOU5 --- Nonhomogeneous OU process with variance function: v(t) = a0 + a1*t + a2*t^2, NOU6 --- Nonhomogeneous OU process with variance function: v(t) = a0 + a1*g(t), NOU7 --- Nonhomogeneous OU process with variance function: v(t) = a0 + a1*g1(t) + a2*g2(t), where g1(t), g2(t) must be two pieces of curves, IOU --- Integrated OU process, W --- Wiener process (Brownian motion), IW --- Integrated Wiener process; time = time variable representing the order of obs; not required for AR, but the data should be sorted by time; TIME MUST BE NONNEGATIVE IF precess = IOU, W OR IW; gt = g(t), the function in NOU1, defined by users. The default is g(t) = t; g1t = g1(t) in NOU4 and NOU7; g2t = g2(t) in NOU4 and NOU7; id = identification variable for subjects; if missing, then the data are assumed to be independent; msrerr = whether or not include measurement error (Y/N); maxiter = maximum number of iterations (maximum=50); conv = convergence criterion using relative difference of the likelihood; print = flag requesting output for model fitting information and variance components be printed (Y/N); outbeta = output data set storing the estimates of beta; outbstd = output data set of std of beta; outvar = output data set for smoothing parameter and variance components; outvstd = output data set for std of variance components; outran_F = output data set storing the estimates of random effects and their frequentist variance; outran_B = output data set storing the estimates of random effects and their Bayesian variance; outprs = output data set storing the estimated Gaussian process; outsmth = output data set storing the est of f(.) in 203 equally spaced points, with two more boundary points added; outband = output data set storing the est std and band for f(.) at the knots. outscore = output data set storing the score test related info. symsize = size of symbolspace. worksize = size of workspace. Examples: Here we used the progesterone data used in the paper cited in the Reference and fitted different models. All the models had random intercept, Gaussian stochastic process and measurement error. The differeces were the choice of the process. There are usually three pages of output. The first page is about the model fitting information; the second page is about the variance components; the third page is about linear fixed effects. Here we only gave all output for AR(1) process, and gave only second page putput for other process. The estimated nonparametric function should be output using outsmth= and its bands as well as SEs should be output using outband=. They can be plotted using S-plus or Proc GPLOT in SAS. When there are continuous fixed effects covariates, for computational stability, it is recomemded that these covariates be centered. 1. AR(1) process. a. program: %include 'spmm.mac'; %spmm(data=control, dep=logpdg, smthvar=day, fixed= age bmi / s, random=int, process=ar, msrerr = y, time=, id=bmid, lintest=y); b. output: The SAS System 1 10:45 Friday, October 25, 1996 Macro: SPMM -- Semi-Parametric Mixed Models Iteration History ITERNUM LIKELIHD OBJECTIV 0 -643.1215 . 1 -425.5509 0.5112682 2 -425.0478 0.0011835 3 -425.0461 4.0446E-6 Model Fitting Information for logpdg DESCRIPT VALUE REML Log Likelihood -425.0461 Number of Subjects 34 Number of Observations 492 Number of Knots 98 Covariance Matrix for Random Effects D INT INT 0.1155396 Score test of linearity for covariate day VERSION CHISQ DF P_VALUE Bias-uncorrected 18.371578 1.6222555 0.0001 Bias-corrected 19.208323 1.6185349 0.0000 SPMM Macro: dep=logpdg, id=bmid, smoothed variable=day fixed effects=age bmi random effects=int Gaussian process=AR(1) The SAS System 2 Estimates of the variance components 10:45 Friday, October 25, 1996 COV_PARM ESTIMATE SE CHISQ P_VALUE SMOOTH 8.0202 4.7525 2.8480 0.0915 D(1,1) 0.1155 0.0843 1.8769 0.1707 AR 0.7680 0.0866 78.6704 0.0000 Diag 0.3408 0.0730 21.8100 0.0000 Residual 0.1145 0.0279 16.8717 0.0000 SPMM Macro: dep=logpdg, id=bmid, smoothed variable=day fixed effects=age bmi random effects=int Gaussian process=AR(1) The SAS System 3 Estimates of the fixed effects 10:45 Friday, October 25, 1996 XNAMES BETA F_SE F_CHISQ F_P B_SE B_CHISQ B_P AGE 1.8464 1.8016 1.0504 0.3054 1.8017 1.0503 0.3054 BMI -2.1479 2.1962 0.9565 0.3281 2.1964 0.9563 0.3281 c. interpretaion of the output: The first page of the putput is about the iteration history and the model fitting information, where the REML Log Likelihood is the REML log likelihood of the smoothing parameter and variance components, as well as the linearity test, if requested. The second page is the estimates of 1/(smoothing parameter), called SMOOTH, and the parameters in the random effects, stochastic process and measurement error. The following is the meaning of the parameters: SMOOTH --- estimate of one over smoothing parameter D(1,1) --- the (1,1)th element of D matrix of random effects AR --- the correlation between two observations in AR(1) process Diag --- marginal variance of the AR(1) process Residual --- the variance of the measurement error The third page contains the information about estimation and inference about beta in the linear fixed effects. F_SE, F_CHISQ, F_P are the standard error, chi-square, p-value from Frequentist inference, while B_SE, B_CHISQ, B_P are the standard error, chi-square, p-value from Bayesian inference. At the end of each page, you can find the information about dependent variable, the variable that needs smoothed, linear fixed effects covariates, random effects covariates, the name of the Gaussian process and the corresponding variance function. 2. Ornstein-Uhlenbeck (OU) process a. program %spmm(data=control, dep=logpdg, smthvar=day, fixed= age bmi / s, random=int, process=ou, msrerr = y, time=day, id=bmid); b. output. (The first page and the third page of the output are similar. We only presented the second page about the variance components information.) 5 Estimates of the variance components 10:45 Friday, October 25, 1996 COV_PARM ESTIMATE SE CHISQ P_VALUE SMOOTH 8.2500 5.0348 2.6850 0.1013 D(1,1) 0.1714 0.0731 5.4945 0.0191 Rho 0.1427 0.1133 1.5874 0.2077 Diag 0.2662 0.0505 27.7506 0.0000 Residual 0.1267 0.0327 15.0153 0.0001 SPMM Macro: dep=logpdg, id=bmid, smoothed variable=day fixed effects=age bmi random effects=int Gaussian process=Ornstein-Uhlenbeck process c. interpretation: SMOOTH --- estimate of one over smoothing parameter D(1,1) --- the (1,1)th element of D matrix of random effects Rho --- the rho in cov(U(t), U(s)) = rho^|t-s| Diag --- marginal variance of the OU process Residual --- the variance of the measurement error 3. NOU1: Nonhomogeneous OU process with log-variance function: log(V(t)) = a0 + a1*g(t) a. program: %spmm(data=control, dep=logpdg, smthvar=day, fixed= age bmi / s, random=int, process=nou1, msrerr = y, time=day, id=bmid, lintest); NOTE: Here gt is missing, so gt=t is assumed. b. output: 8 Estimates of the variance components 10:45 Friday, October 25, 1996 COV_PARM ESTIMATE SE CHISQ P_VALUE SMOOTH 7.7999 4.7710 2.6728 0.1021 D(1,1) 0.2364 0.0723 10.6849 0.0011 Rho 0.0690 0.0544 1.6086 0.2047 A0 -1.7849 0.2572 48.1735 0.0000 A1 0.9792 0.2497 15.3814 0.0001 Residual 0.1166 0.0264 19.4554 0.0000 SPMM Macro: dep=logpdg, id=bmid, smoothed variable=day fixed effects=age bmi random effects=int Gaussian process=Nonhomogeneous OU process with log-variance function: log(V(t)) = a0 + a1*g(t) c. interpretation: SMOOTH --- estimate of one over smoothing parameter D(1,1) --- the (1,1)th element of D matrix of random effects Rho --- the rho in cov(U(t), U(s)) = rho^|t-s| A0 --- ao in log(V(t)) = a0 + a1*g(t) A1 --- a1 in log(V(t)) = a0 + a1*g(t) Residual --- the variance of the measurement error 4. NOU2: Nonhomogeneous OU process with log-variance function: log(V(t)) = a0 + a1*exp(alpha*t) a. program: %spmm(data=control, dep=logpdg, smthvar=day, fixed= age bmi / s, random=int, process=nou2, msrerr = y, time=day, id=bmid, lintest=y); b. output: 11 Estimates of the variance components 10:45 Friday, October 25, 1996 COV_PARM ESTIMATE SE CHISQ P_VALUE SMOOTH 22.7023 15.0126 2.2868 0.1305 D(1,1) 0.0000 0.9805 0.0000 1.0000 Rho 0.7942 0.0523 230.742 0.0000 Alpha -1.1543 0.2761 17.4776 0.0000 A0 0.8766 0.4835 3.2876 0.0698 A1 0.7076 0.3243 4.7623 0.0291 Residual 0.6047 0.0799 57.3022 0.0000 SPMM Macro: dep=logpdg, id=bmid, smoothed variable=day fixed effects=age bmi random effects=int Gaussian process=Nonhomogeneous OU process with log-variance function: log(V(t)) = a0 + a1*exp(alpha*t) c. interpretation: SMOOTH --- estimate of one over smoothing parameter D(1,1) --- the (1,1)th element of D matrix of random effects Rho --- the rho in cov(U(t), U(s)) = rho^|t-s| Alpha --- the alpha in log(V(t)) = a0 + a1*exp(alpha*t) A0 --- the a0 in log(V(t)) = a0 + a1*exp(alpha*t) A1 --- the A1 in log(V(t)) = a0 + a1*exp(alpha*t) Residual --- the variance of the measurement error 5. NOU3: Nonhomogeneous OU process with log-variance function: log(V(t)) = a0 + a1*t + a2*t^2 a. program: %spmm(data=control, dep=logpdg, smthvar=day, fixed= age bmi / s, random=int, process=nou3, msrerr = y, time=day, id=bmid, lintest=y); b. output: Estimates of the variance components 10:45 Friday, October 25, 1996 COV_PARM ESTIMATE SE CHISQ P_VALUE SMOOTH 7.3159 4.5124 2.6286 0.1050 D(1,1) 0.2598 0.0719 13.0513 0.0003 Rho 0.0968 0.0613 2.4955 0.1142 A0 -2.1034 0.4206 25.0069 0.0000 A1 3.5210 1.1524 9.3351 0.0022 A2 -2.0745 0.7708 7.2433 0.0071 Residual 0.1360 0.0166 67.4678 0.0000 SPMM Macro: dep=logpdg, id=bmid, smoothed variable=day fixed effects=age bmi random effects=int Gaussian process=Nonhomogeneous OU process with log-variance function: log(V(t)) = a0 + a1*t + a2*t^2 c. interpretation: SMOOTH --- estimate of one over smoothing parameter D(1,1) --- the (1,1)th element of D matrix of random effects Rho --- the rho in cov(U(t), U(s)) = rho^|t-s| A0 --- the a0 in log(V(t)) = a0 + a1*t + a2*t^2 A1 --- the a1 in log(V(t)) = a0 + a1*t + a2*t^2 A2 --- the a2 in log(V(t)) = a0 + a1*t + a2*t^2 Residual --- the variance of the measurement error 6. IOU: Integrated OU process a. program: %spmm(data=control, dep=logpdg, smthvar=day, fixed= age bmi / s, random=int, process=iou, msrerr = y, time=day, id=bmid, lintest=y); NOTE: the variable in time= should be nonnegative b. output: 17 Estimates of the variance components 10:45 Friday, October 25, 1996 COV_PARM ESTIMATE SE CHISQ P_VALUE SMOOTH 9.0315 5.3381 2.8625 0.0907 D(1,1) 0.0917 0.0967 0.8990 0.3430 Alpha 26.9190 133.544 0.0406 0.8402 Sigma^2 6.1340 30.1805 0.0413 0.8389 Residual 0.1652 0.0712 5.3788 0.0204 SPMM Macro: dep=logpdg, id=bmid, smoothed variable=day fixed effects=age bmi random effects=int Gaussian process=Integrated OU process c. interpretation: SMOOTH --- estimate of one over smoothing parameter D(1,1) --- the (1,1)th element of D matrix of random effects Alpha --- -log(rho), where rho is the rho in the OU process Sigma^2 --- the variance of the OU process Residual --- the variance of the measurement error 7. W: Wiener process (Brownian motion) a. program: %spmm(data=control, dep=logpdg, smthvar=day, fixed= age bmi / s, random=int, process=w, msrerr = y, time=day, id=bmid, lintest=y); NOTE: the variable in time= should be nonnegative b. output: 20 Estimates of the variance components 10:45 Friday, October 25, 1996 COV_PARM ESTIMATE SE CHISQ P_VALUE SMOOTH 8.7015 5.2265 2.7719 0.0959 D(1,1) 0.1165 0.0755 2.3780 0.1231 Sigma^2 0.4334 0.0825 27.5880 0.0000 Residual 0.1735 0.0179 93.4297 0.0000 SPMM Macro: dep=logpdg, id=bmid, smoothed variable=day fixed effects=age bmi random effects=int Gaussian process=Wiener process (Brownian motion) c. interpretation: SMOOTH --- estimate of one over smoothing parameter D(1,1) --- the (1,1)th element of D matrix of random effects Sigma^2 --- multiple of the Brownian motion Residual --- the variance of the measurement error 8. IW: Integrated Wiener process a. program: %spmm(data=control, dep=logpdg, smthvar=day, fixed= age bmi / s, random=int, process=iw, msrerr = y, time=day, id=bmid); NOTE: the variable in time= should be nonnegative b. output: 23 Estimates of the variance components 10:45 Friday, October 25, 1996 COV_PARM ESTIMATE SE CHISQ P_VALUE SMOOTH 8.8050 5.4335 2.6260 0.1051 D(1,1) 0.2806 0.0844 11.0554 0.0009 Sigma^2 0.0939 0.0295 10.1456 0.0014 Residual 0.2755 0.0194 202.667 0.0000 SPMM Macro: dep=logpdg, id=bmid, smoothed variable=day fixed effects=age bmi random effects=int Gaussian process=Integrated Wiener process c. interpretation: SMOOTH --- estimate of one over smoothing parameter D(1,1) --- the (1,1)th element of D matrix of random effects Sigma^2 --- Sigma^2 the Wiener process Residual --- the variance of the measurement error Reference: Zhang, D., Lin, X., Raz, J. and M.F. Sowers, (1998), Semiparametric Stochastic Mixed Models for Longitudinal Data, JASA, 93, 710-719. CAUTION: All effects should be defined except the fixed effect of intercept. There must be at least 3 distinct tij's in order to fit f(.). Otherwise the program would crash. Problems --- Report problems to Daowen Zhang at dzhang@umich.edu