Nothing
#' Estimate average hazard ratios from k independent samples
#'
#' This function is a simple wrapper for \code{\link{ahrWKM}} and \code{\link{ahrAJ}}.
#'
#' @title avgHR
#' @param L time-limit specifying time-interval [0,L] on which average hazard ratios will be calculated
#' @param data data frame (see \code{data} argument to \code{\link{ahrWKM}} (if \code{method == "wkm" || "km"}) or \code{\link{ahrAJ}} (if \code{method == "aj"}))
#' @param method method used for estimating survival functions (default: Kaplan-Meier estimator)
#' @param ... additional arguments passed to \code{\link{ahrWKM}} or \code{\link{ahrAJ}}
#' @return An object of class '"ahr"'
#' @seealso \code{\link{ahrWKM}}, \code{\link{ahrAJ}}
#' @references J.~D. Kalbfleisch and R.~L. Prentice. Estimation of the average hazard ratio. \emph{Biometrika}, 68(1):105--112, Apr. 1981.
#' @export
#' @examples
#' T <- c(rexp(100, 1), rexp(100, 2))
#' C <- c(rexp(100, 1), rexp(100, 2))
#' Y <- pmin(T, C)
#' D <- T <= C
#' Z <- rep(c(0,1), c(100, 100))
#' fit <- avgHR(2, data.frame(Y=Y, D=D, Z=Z), formula=Surv(Y, D) ~ Z)
avgHR <- function(L, data, method="km", ...) {
if(method %in% c("km", "wkm")) ahrWKM(L=L, data=data, ...)
else if(method == "aj") ahrAJ(L=L, data=data, ...)
else if(method == "user") ahrUser(L=L, data=data, ...)
else stop(paste("Unknown method:", method))
}
#' Estimate average hazard ratios from k independent samples based on the weighted Kaplan-Meier (WKM) estimator
#'
#' @title ahrWKM
#' @param L time-limit specifying time-interval [0,L] over which average hazard ratios will be calculated
#' @param formula an object of class '"formula"' specifying the conditional survival model
#' @param data data frame containing the variables in formula
#' @param null.theta vector specifying the null hypothesis for the average hazard ratios (H_0: theta = null.theta)
#' @param contrast vector of contrasts to test H_0: contrast * (theta - null.theta) = 0
#' @param multi.test calculate multivariate test statistic if TRUE
#' @param cov if TRUE calculate covariance matrix estimator (direct)
#' @param bootstrap if > 0 then use bootstrap to estimate covariance matrix (ignore if cov is TRUE)
#' @param alpha exponent of the weight function
#' @param left.limit if TRUE use left-continuous interpolation of WKM estimates instead of right-continuous interpolation
#' @param rr.subset logical vector defining subset of observations to use for response rate estimation (default: use all observations)
#' @return An object of class '"ahr"'
#' @references J.~D. Kalbfleisch and R.~L. Prentice. Estimation of the average hazard ratio. \emph{Biometrika}, 68(1):105--112, Apr. 1981.
#' @export
#' @seealso \code{\link{wkm}}
#' @examples
#' T <- c(rexp(100, 1), rexp(100, 2))
#' C <- c(rexp(100, 1), rexp(100, 2))
#' Y <- pmin(T, C)
#' D <- T <= C
#' Z <- rep(c(0,1), c(100, 100)) # treatment indicator
#' fit <- ahrWKM(2, Surv(Y, D) ~ Z, data.frame(Y=Y, D=D, Z=Z))
#' fit
#'
#' ## the same as above, but estimate covariance matrix using bootstrap
#' \dontrun{fitBS <- ahrWKM(2, Surv(Y, D) ~ Z, data.frame(Y=Y, D=D, Z=Z), cov=FALSE,
#' bootstrap=1000)}
#'
ahrWKM <- function(L, formula, data, null.theta=NULL, contrast=NULL, multi.test=FALSE, cov=TRUE, bootstrap=0, alpha=1, left.limit=FALSE, rr.subset=rep(TRUE, nrow(data))) {
if(!is.null(formula)) data <- parseFormula(formula, data)
wkm.param <- list(alpha=alpha, var=cov, cov=cov, left.limit=left.limit, rr.subset=rr.subset)
fit <- ahrSurv(L, data, null.theta, contrast, multi.test, cov, bootstrap, wkm, wkm.param)
fit <- c(fit, logHR(fit))
class(fit) <- "ahr"
fit
}
#' Estimate average hazard ratios from k independent samples based on the Kaplan-Meier estimator
#'
#' @title ahrKM
#' @param L time-limit specifying time-interval [0,L] over which average hazard ratios will be calculated
#' @param formula an object of class '"formula"' specifying the conditional survival model
#' @param data data frame containing the variables in formula
#' @param null.theta vector specifying the null hypothesis for the average hazard ratios (H_0: theta = null.theta)
#' @param contrast vector of contrasts to test H_0: contrast * (theta - null.theta) = 0
#' @param multi.test calculate multivariate test statistic if TRUE
#' @param cov if TRUE calculate covariance matrix estimator (direct)
#' @param bootstrap if > 0 then use bootstrap to estimate covariance matrix (ignore if cov is TRUE)
#' @param left.limit if TRUE use left-continuous interpolation of WKM estimates
#' @return An object of class '"ahr"'
#' @references J.~D. Kalbfleisch and R.~L. Prentice. Estimation of the average hazard ratio. \emph{Biometrika}, 68(1):105--112, Apr. 1981.
#' @export
#' @seealso \code{\link{survfit}}
#' @examples
#' T <- c(rexp(100, 1), rexp(100, 2))
#' C <- c(rexp(100, 1), rexp(100, 2))
#' Y <- pmin(T, C)
#' D <- T <= C
#' Z <- rep(c(0,1), c(100, 100)) # treatment indicator
#' fit <- ahrKM(2, Surv(Y, D) ~ Z, data.frame(Y=Y, D=D, Z=Z))
#' fit
#'
#' ## the same as above, but estimate covariance matrix using bootstrap
#' \dontrun{fitBS <- ahrKM(2, Surv(Y, D) ~ Z, data.frame(Y=Y, D=D, Z=Z), cov=FALSE,
#' bootstrap=1000)}
ahrKM <- function(L, formula, data, null.theta=NULL, contrast=NULL, multi.test=FALSE, cov=TRUE, bootstrap=0, left.limit=FALSE) {
if(!is.null(formula)) data <- parseFormula(formula, data)
wkm.param <- list(alpha=1, var=cov, cov=FALSE, left.limit=left.limit, rr.subset=rep(TRUE, nrow(data)))
fit <- ahrSurv(L, data, null.theta, contrast, multi.test, cov, bootstrap, wkm, wkm.param, TRUE)
fit <- c(fit, logHR(fit))
class(fit) <- "ahr"
fit
}
#' Estimate average hazard ratios from k independent samples based on the Aalen-Johansen estimator of the empirical transition probabilities (NOTE: variance estimation not yet implemented)
#'
#' @title ahrAJ
#' @param L time-limit specifying time-interval [0,L] over which average hazard ratios will be calculated
#' @param target string specifying the target transition, for which the Aalen-Johansen estimator is to be calculated
#' @param states list of state names
#' @param transitions matrix of possible transitions
#' @param censoring name of censoring 'state'
#' @param data data frame containing variables id, time, from, to (see \code{\link{etm}}) and Trt (factor giving treatment groups)
#' @param null.theta vector specifying the null hypothesis for the average hazard ratios
#' @param contrast vector of contrasts to test H_0: contrast * (theta - null.theta) = 0
#' @param multi.test calculate multivariate test statistic if TRUE
#' @param cov if TRUE calculate covariance matrix estimator (direct)
#' @param bootstrap number of bootstrap samples to draw for variance estimation (default: 0 = no bootstrap, direct variance estimation). This parameter is ignored if cov=TRUE
#' @return An object of class '"ahr"'
#' @seealso \code{\link{aj}}
#' @references J.~D. Kalbfleisch and R.~L. Prentice. Estimation of the average hazard ratio. \emph{Biometrika}, 68(1):105--112, Apr. 1981.
#' @export
#' @examples
#' ## competing risks
#' Trt <- factor(rep(c(0,1), c(100, 100)))
#' T <- c(rexp(100, 1), rexp(100, 2))
#' C <- c(rexp(100, 1), rexp(100, 2))
#' r <- c(rbinom(100, 2, 0.5), rbinom(100, 2, 0.4))
#' r[(r == 0) | (T > C)] <- "cens"
#' data <- data.frame(id=1:200, time=pmin(T,C), from=rep(0, 200), to=r, Trt=Trt)
#' tra <- matrix(FALSE, nrow=3, ncol=3)
#' tra[1, 2:3] <- TRUE
#' # estimate average subdistribution hazard ratio up to L=2 for event type 1
#' fit <- ahrAJ(2, target="0 1", states=c("0", "1", "2"), transitions=tra, censoring="cens",
#' data=data, cov=TRUE)
#' fit
ahrAJ <- function(L, target, states, transitions, censoring, data, null.theta=NULL, contrast=NULL, multi.test=FALSE, cov=FALSE, bootstrap=0) {
msm <- list(target=target, states=states, transitions=transitions, censoring=censoring, s=0, t=L)
if(cov) {
if(isCmprsk(msm$transitions)) msm$covariance <- TRUE
else stop("Variance estimation only supported for competing risks models")
if(bootstrap > 0) warning("Bootstrap parameter ignored since 'cov' is TRUE")
} else msm$covariance <- FALSE
fit <- ahrSurv(L, data, null.theta, contrast, multi.test, cov, bootstrap, aj, msm, FALSE)
fit <- c(fit, logHR(fit))
class(fit) <- "ahr"
fit
}
#' Estimate average hazard ratios from k independent samples based on user-supplied survival function estimator
#'
#' @title ahrUser
#' @param L time-limit specifying time-interval [0,L] over which average hazard ratios will be calculated
#' @param formula an object of class '"formula"' specifying the conditional survival model
#' @param data data frame containing the variables in formula
#' @param null.theta vector specifying the null hypothesis for the average hazard ratios (H_0: theta = null.theta)
#' @param contrast vector of contrasts to test H_0: contrast * (theta - null.theta) = 0
#' @param multi.test calculate multivariate test statistic if TRUE
#' @param cov if TRUE calculate covariance matrix estimator (direct)
#' @param bootstrap if > 0 then use bootstrap to estimate covariance matrix (ignore if cov is TRUE)
#' @param user.survfit user defined function taking vector of times, data.frame and list of parameters returning survival function estimate
#' @param user.param list of parameters passed to function user.survfit
#' @return An object of class '"ahr"'
#' @references J.~D. Kalbfleisch and R.~L. Prentice. Estimation of the average hazard ratio. \emph{Biometrika}, 68(1):105--112, Apr. 1981.
#' @export
#' @seealso \code{\link{wkm}}
#' @details user.survfit must return logV (if user supplied survival function estimator hast independent increments property) or logCOV
#' @examples
#' ## User supplied survival function estimator (should be exactly the same as
#' T <- c(rexp(100, 1), rexp(100, 2))
#' C <- c(rexp(100, 1), rexp(100, 2))
#' time <- pmin(T, C)
#' status <- T <= C
#' trt <- rep(c(0,1), c(100, 100)) # treatment indicator
#'
#' sfit <- function(times, data, param) {
#' fit <- survfit(Surv(Y, D) ~ 1, data=data)
#' f <- approxfun(fit$time, fit$surv, method="constant", f=0, yleft=1, rule=2)
#' fv <- approxfun(fit$time, fit$std.err^2, method="constant", f=0, yleft=1, rule=2)
#'
#' S <- f(times)
#' logV <- fv(times) * nrow(data)
#' V <- S^2 * logV
#'
#' list(times=times, S=S, V=V, logV=logV)
#' }
#' fit1 <- ahrUser(2, Surv(time, status) ~ trt,
#' data.frame(time=time, status=status, trt=trt), user.survfit=sfit, user.param=list())
#' fit1
#' fit2 <- ahrKM(2, Surv(time, status) ~ trt, data.frame(time=time, status=status, trt=trt), cov=FALSE)
#' fit2
#'
ahrUser <- function(L, formula, data, null.theta=NULL, contrast=NULL, multi.test=FALSE, cov=FALSE, bootstrap=0, user.survfit=wkm, user.param=list(alpha=1, var=FALSE, cov=FALSE, left.limit=FALSE, rr.subset=rep(TRUE, nrow(data)))) {
if(!is.null(formula)) data <- parseFormula(formula, data)
fit <- ahrSurv(L, data, null.theta, contrast, multi.test, cov, bootstrap, user.survfit, user.param, TRUE)
fit <- c(fit, logHR(fit))
class(fit) <- "ahr"
fit
}
#' Print ahr object
#'
#' @title print.ahr
#' @param x an object of class '"ahr"'.
#' @param digits minimal number of significant digits.
#' @param ... further arguments passed to or from other methods.
#' @method print ahr
#' @export
print.ahr <- function(x, digits=3, ...) {
cat("\n")
cat("\t Average Group to Total Hazard Ratios ", paste0("(L=", x$L, ")"), "\n")
cat("\n")
cat(paste(x$k, "treatment groups:"), levels(x$groups), "\n")
cat("\n")
if(is.null(x$cov.theta)) tmpcov <- rep(NA, x$k) ## covariance matrix missing
else if(is.null(dim(x$cov.theta))) tmpcov <- x$cov.theta ## scalar value
else tmpcov <- diag(x$cov.theta)
if(is.null(x$Z.theta)) Ztheta <- NA
else Ztheta <- round(x$Z.theta, digits)
tmp <- cbind(round(x$theta, digits), round(sqrt(tmpcov), digits), Ztheta)
dimnames(tmp) <- list(levels(x$groups), c("Theta", "Std.dev", "Z"))
print(tmp)
cat("\n\n")
cat("\t Average Hazard Ratios ", paste0("(L=", x$L, ")"), "\n")
cat("\n")
if(is.null(x$cov.hr)) tmpcov <- rep(NA, x$k-1)
else if(is.null(dim(x$cov.hr))) tmpcov <- x$cov.hr
else tmpcov <- diag(x$cov.hr)
cat("\n")
if(is.null(x$Z.hr)) Zhr <- NA
else Zhr <- round(x$Z.hr, digits)
tmp <- cbind(round(x$hr, digits), round(sqrt(tmpcov), digits), Zhr)
dimnames(tmp) <- list(levels(x$groups)[-1], c("HR", "Std.dev", "Z"))
print(tmp)
cat("\n")
if(!is.null(x$contrast)) cat
invisible(x)
}
## INTERNAL
##
## title ahrSurv
## param data data frame containing the variables
## param L time-limit specifying time-interval [0,L] over which average hazard ratios will be calculated
## param null.theta vector specifying the null hypothesis for the average hazard ratios
## param contrast vector of contrasts to test H_0: contrast * (theta - null.theta) = 0
## param multi.test calculate multivariate test statistic if TRUE
## param cov if TRUE calculate asymptotic covariance matrix estimator (direct estimation)
## param bootstrap if > 0 then use bootstrap to estimate covariance matrix (ignored if cov is TRUE)
## param surv.fit.fun function which calculates survival function from data
## param surv.fit.param list of parameters passed to surv.fit.fun
## param log.iis asymptotic covariance function of logarithm of survival function has independent increments structure (TRUE for (weighted) Kaplan-Meier estimator, FALSE for Aalen-Johansen estimator)
## return an object of class '"ahr"'
ahrSurv <- function(L, data, null.theta=NULL, contrast=NULL, multi.test=FALSE, cov=TRUE, bootstrap=0, surv.fit.fun=wkm, surv.fit.param=list(alpha=1, cov=TRUE), log.iis=FALSE) {
ahr.obj <- ahrFit(data, L, cov, surv.fit.fun, surv.fit.param, log.iis)
if(!cov && bootstrap > 0) ahr.obj$cov.theta <- ahrBootstrap(ahr.obj, data, bootstrap, surv.fit.fun, surv.fit.param)$cov.theta
if(cov && bootstrap > 0) warning("Bootstrap parameter ignored since 'cov' is TRUE")
k <- ahr.obj$k
n <- ahr.obj$n
theta <- ahr.obj$theta
cov.theta <- ahr.obj$cov.theta
if(is.null(null.theta)) null.theta <- rep(1/k, k)
ahr.obj$null.theta <- null.theta
if(!is.null(cov.theta)) {
## Wald-type test
Z.theta <- (theta - null.theta) / sqrt(diag(cov.theta))
ahr.obj$cov.theta <- cov.theta
ahr.obj$Z.theta <- Z.theta
## linear contrast test
if(!is.null(contrast)) {
ahr.obj$contrast <- contrast
ahr.obj$Z.contrast <- contrast %*% (theta - null.theta) / sqrt(contrast %*% cov.theta %*% contrast)
}
if(multi.test) {
ahr.obj$multi.test <- multi.test
ahr.obj$Z.multi <- (theta - null.theta) %*% MASS::ginv(cov.theta) %*% (theta - null.theta) ## has chi^2 distribution with k-1 degrees of freedom
}
}
ahr.obj
}
## INTERNAL
## bootstrap estimation of covariance matrix
ahrBootstrap <- function(ahr.obj, data, B, surv.fit.fun, surv.fit.param) {
k <- ahr.obj$k
lev <- levels(ahr.obj$groups)
L <- ahr.obj$L
sel <- list()
n <- numeric(k)
if(is.null(data$id)) data$id <- 1:nrow(data)
## sel <- lapply(1:k, function(i) unique(data$id[data$Trt == lev[i]]))
## n <- do.call(c, lapply(sel, function(x) length(x)))
for(i in 1:k) {
sel[[i]] <- unique(data$id[data$Trt == lev[i]])
n[i] <- length(sel[[i]])
}
surv.fit.param$var <- FALSE
surv.fit.param$cov <- FALSE
f <- function() {
## stratified resampling
ids <- do.call(c, lapply(1:k, function(i) sample(sel[[i]], n[i], replace=TRUE)))
## draw bootstrap sample
Bdata <- data[data$id %in% ids,]
Bdata$id <- 1:nrow(Bdata) ## ids are not unique after sampling WITH replacement
## calculate ahr from bootstrap sample
ahrFit(Bdata, L, FALSE, surv.fit.fun, surv.fit.param)$theta
}
res <- replicate(B, f())
list(theta=rowMeans(res), cov.theta=cov(t(res)))
}
## INTERNAL
ahrFit <- function(data, L, cov, surv.fit.fun, surv.fit.param, log.iis=FALSE) {
grps <- levels(data$Trt)
k <- length(grps)
if(k <= 1) stop("Need at least two groups!")
times <- getTimes(L, data)
n.times <- length(times)
n <- NULL
fit <- list()
trt.sub <- data$Trt[surv.fit.param$rr.subset]
## estimate transition probabilities in each treatment group
for(i in grps) {
cur.data <- data[data$Trt == i,]
par <- surv.fit.param
par$rr.subset <- surv.fit.param$rr.subset[trt.sub == i]
n <- c(n, nrow(cur.data))
fit[[length(fit)+1]] <- surv.fit.fun(times, cur.data, par)
}
p <- n / sum(n)
tmp <- lapply(1:k, function(i) fit[[i]]$S)
## dummy to ensure length(tmp[-l]) > 1
tmp[[k+1]] <- rep.int(1, n.times)
xi <- sapply(1:k, function(l) do.call(function(...) mapply(prod, ...), tmp[-l]))
G <- xi[,k] * fit[[k]]$S
x <- sapply(1:(k-1), function(i) stepIntegrate(xi[,i], fit[[i]]$S))
theta <- -x / (1 - G[n.times])
theta[k] <- 1 - sum(theta)
ahr.obj <- list(L=L, surv.fit.param=surv.fit.param, k=k, n=n, p=p,
surv.fit=fit, times=times, n.times=n.times, groups=data$Trt, strata=data$W, log.iis=log.iis)
class(ahr.obj) <- "ahr"
ahr.obj$theta <- theta
if(cov) ahr.obj$cov.theta <- ahrVar(theta, G, xi, times, p, fit, log.iis) / sum(n)
ahr.obj
}
## INTERNAL
ahrVar <- function(theta, G, xi, times, p, fit, log.iis=FALSE) {
k <- length(fit)
GL <- G[length(G)]
n.times <- length(times)
snt <- 1:n.times
## logCOV = COV / S
## COV == 0 <=> S = 0
get.A <- function(i, j) {
if(log.iis) ft <- xi[,j] * fit[[i]]$logV
else ft <- xi[,j] * fit[[i]]$logCOV[ , ncol(fit[[i]]$logCOV)]
GL * stepIntegrate(ft, fit[[j]]$S) / p[i]
}
get.B <- function(i, j, kk) {
Sj <- fit[[j]]$S
if(log.iis) {
logV <- fit[[i]]$logV
tmp <- simplify2array(lapply(snt, function(l) stepIntegrate(xi[,j] * logV[pmin(snt, l)], Sj)))
} else {
logCOV <- fit[[i]]$logCOV
tmp <- simplify2array(lapply(snt, function(l) stepIntegrate(xi[,j] * logCOV[, l], Sj)))
}
stepIntegrate(xi[,kk] * tmp, fit[[kk]]$S) / p[i]
}
get.C <- function(i) {
if(log.iis) GL^2 * fit[[i]]$logV[n.times] / p[i]
else GL^2 * fit[[i]]$logCOV[n.times, n.times] / p[i]
}
A <- matrix(0, nrow=k, ncol=k)
B <- array(0, dim=rep.int(k, 3))
C <- numeric(k)
sel <- 1:k
for(i in sel) {
C[i] <- get.C(i)
for(j in sel) {
if(i != j) A[i,j] <- get.A(i,j)
for(l in j:k) {
B[i,j,l] <- get.B(i,j,l)
B[i,l,j] <- B[i,j,l]
}
}
}
VG <- sum(C)
VxG <- numeric(k)
Sigma <- matrix(0, nrow=k, ncol=k)
for(i in sel) {
si <- sel[-i]
Vxii <- sum(B[i, si, si]) + sum(B[si, i, i]) + C[i] - 2*sum(A[i, si])
VxG[i] <- sum(A[i, si] - A[si, i]) - C[i]
Sigma[i,i] <- Vxii + 2*theta[i]*VxG[i] + theta[i]^2 * VG
}
for(i in sel[-k]) {
si <- sel[-i]
for(j in sel[(i+1):k]) {
Vxij <- -sum(B[i, si, j]) - sum(B[j, sel[-j], i]) + sum(B[sel[-c(i,j)], i, j]) + A[i,j] + A[j,i]
Sigma[i,j] <- Vxij + theta[j]*VxG[i] + theta[i]*VxG[j] + theta[i]*theta[j]*VG
Sigma[j,i] <- Sigma[i,j]
}
}
## must be zero:
##print(rowSums(Sigma))
##print(colSums(Sigma))
Sigma / (1 - GL)^2
}
## INTERNAL
logHR <- function(ahr.obj, null.hr=1, contrast=NULL) {
k <- ahr.obj$k
theta <- ahr.obj$theta
cov.theta <- ahr.obj$cov.theta
null.loghr <- log(null.hr)
hr <- theta[-1] / theta[1]
loghr <- log(hr)
res <- list(hr=hr, loghr=loghr, null.hr=null.hr, null.loghr=null.loghr)
if(!is.null(cov.theta)) {
## covariance matrix cov.beta has dimension k-1
if(k == 2) {
cov.loghr <- cov.theta[1,1] / (theta[1] * theta[2])^2
Z.loghr <- (loghr - null.loghr) / sqrt(cov.loghr)
cov.hr <- cov.loghr * hr^2
Z.hr <- (hr - null.hr) / sqrt(cov.hr)
} else {
J <- matrix(0, nrow=k-1, ncol=k)
J[,1] <- -1 / theta[1]
J[,-1] <- diag(1 / theta[-1])
cov.loghr <- J %*% cov.theta %*% t(J)
Z.loghr <- sum(solve(cov.loghr, loghr - null.loghr) * (loghr - null.loghr)) ## asympt. chi^2 distribution with k-1 degrees of freedom
J2 <- diag(loghr)
cov.hr <- J2 %*% cov.loghr %*% J2
Z.hr <- sum(solve(cov.hr, hr - null.hr) * (hr - null.hr)) ## asympt. chi^2 distribution with k-1 degrees of freedom
}
res$Z.loghr <- Z.loghr
res$cov.loghr <- cov.loghr
res$Z.hr <- Z.hr
res$cov.hr <- cov.hr
if(!is.null(contrast)) { ## asympt. standard normal distribution
res$contrast <- contrast
res$Z.contrast.loghr <- sum(contrast * (loghr - null.loghr)) / sqrt(contrast %*% cov.loghr %*% contrast)
res$Z.contrast.hr <- sum(contrast * (hr - null.hr)) / sqrt(contrast %*% cov.hr %*% contrast)
}
}
res
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.