# Generate multivariate normal responses with equal correlations (.7)
# within subjects and no correlation between subjects
# Simulate realizations from a piecewise linear population time-response
# profile with large subject effects, and fit using a 6-knot spline
# Estimate the correlation structure from the residuals, as a function
# of the absolute time difference
# Function to generate n p-variate normal variates with mean vector u and
# covariance matrix S
# Slight modification of function written by Bill Venables
# See also the built-in function rmvnorm
mvrnorm <- function(n, p = 1, u = rep(0, p), S = diag(p)) {
Z <- matrix(rnorm(n * p), p, n)
t(u + t(chol(S)) %*% Z)
}
n <- 20 # Number of subjects
sub <- .5*(1:n) # Subject effects
# Specify functional form for time trend and compute non-stochastic component
times <- seq(0, 1, by=.1)
g <- function(times) 5*pmax(abs(times-.5),.3)
ey <- g(times)
# Generate multivariate normal errors for 20 subjects at 11 times
# Assume equal correlations of rho=.7, independent subjects
nt <- length(times)
rho <- .7
set.seed(19)
errors <- mvrnorm(n, p=nt, S=diag(rep(1-rho,nt))+rho)
# Note: first random number seed used gave rise to mean(errors)=0.24!
# Add E[Y], error components, and subject effects
y <- matrix(rep(ey,n), ncol=nt, byrow=TRUE) + errors +
matrix(rep(sub,nt), ncol=nt)
# String out data into long vectors for times, responses, and subject ID
y <- as.vector(t(y))
times <- rep(times, n)
id <- sort(rep(1:n, nt))
# Show lowess estimates of time profiles for individual subjects
f <- rm.boot(times, y, id, plot.individual=TRUE, B=25, cor.pattern='estimate',
smoother=lowess, bootstrap.type='x fixed', nk=6)
# In practice use B=400 or 500
# This will compute a dependent-structure log-likelihood in addition
# to one assuming independence. By default, the dep. structure
# objective will be used by the plot method (could have specified rho=.7)
# NOTE: Estimating the correlation pattern from the residual does not
# work in cases such as this one where there are large subject effects
# Plot fits for a random sample of 10 of the 25 bootstrap fits
plot(f, individual.boot=TRUE, ncurves=10, ylim=c(6,8.5))
# Plot pointwise and simultaneous confidence regions
plot(f, pointwise.band=TRUE, col.pointwise=1, ylim=c(6,8.5))
# Plot population response curve at average subject effect
ts <- seq(0, 1, length=100)
lines(ts, g(ts)+mean(sub), lwd=3)
## Not run:
# #
# # Handle a 2-sample problem in which curves are fitted
# # separately for males and females and we wish to estimate the
# # difference in the time-response curves for the two sexes.
# # The objective criterion will be taken by plot.rm.boot as the
# # total of the two sums of squared errors for the two models
# #
# knots <- rcspline.eval(c(time.f,time.m), nk=6, knots.only=TRUE)
# # Use same knots for both sexes, and use a times vector that
# # uses a range of times that is included in the measurement
# # times for both sexes
# #
# tm <- seq(max(min(time.f),min(time.m)),
# min(max(time.f),max(time.m)),length=100)
#
#
# f.female <- rm.boot(time.f, bp.f, id.f, knots=knots, times=tm)
# f.male <- rm.boot(time.m, bp.m, id.m, knots=knots, times=tm)
# plot(f.female)
# plot(f.male)
# # The following plots female minus male response, with
# # a sequence of shaded confidence band for the difference
# plot(f.female,f.male,multi=TRUE)
#
#
# # Do 1000 simulated analyses to check simultaneous coverage
# # probability. Use a null regression model with Gaussian errors
#
#
# n.per.pt <- 30
# n.pt <- 10
#
#
# null.in.region <- 0
#
#
# for(i in 1:1000) {
# y <- rnorm(n.pt*n.per.pt)
# time <- rep(1:n.per.pt, n.pt)
# # Add the following line and add ,id=id to rm.boot to use clustering
# # id <- sort(rep(1:n.pt, n.per.pt))
# # Because we are ignoring patient id, this simulation is effectively
# # using 1 point from each of 300 patients, with times 1,2,3,,,30
#
#
# f <- rm.boot(time, y, B=500, nk=5, bootstrap.type='x fixed')
# g <- plot(f, ylim=c(-1,1), pointwise=FALSE)
# null.in.region <- null.in.region + all(g$lower<=0 & g$upper>=0)
# prn(c(i=i,null.in.region=null.in.region))
# }
#
#
# # Simulation Results: 905/1000 simultaneous confidence bands
# # fully contained the horizontal line at zero
# ## End(Not run)
Run the code above in your browser using DataLab