Nothing
###{{{ twinlm
##' Fits a classical twin model for quantitative traits.
##'
##' @title Classic twin model for quantitative traits
##' @return Returns an object of class \code{twinlm}.
##' @author Klaus K. Holst
##' @seealso \code{\link{bptwin}}, \code{\link{twinlm.time}}, \code{\link{twinlm.strata}}, \code{\link{twinsim}}
##' @aliases twinlm twinlm.strata
##' @export
##' @examples
##' ## Simulate data
##' set.seed(1)
##' d <- twinsim(1000,b1=c(1,-1),b2=c(),acde=c(1,1,0,1))
##' ## E(y|z1,z2) = z1 - z2. var(A) = var(C) = var(E) = 1
##'
##' ## E.g to fit the data to an ACE-model without any confounders we simply write
##' ace <- twinlm(y ~ 1, data=d, DZ="DZ", zyg="zyg", id="id")
##' ace
##' ## An AE-model could be fitted as
##' ae <- twinlm(y ~ 1, data=d, DZ="DZ", zyg="zyg", id="id", type="ae")
##' ## LRT:
##' lava::compare(ae,ace)
##' ## AIC
##' AIC(ae)-AIC(ace)
##' ## To adjust for the covariates we simply alter the formula statement
##' ace2 <- twinlm(y ~ x1+x2, data=d, DZ="DZ", zyg="zyg", id="id", type="ace")
##' ## Summary/GOF
##' summary(ace2)
##' \donttest{ ## Reduce Ex.Timings
##' ## An interaction could be analyzed as:
##' ace3 <- twinlm(y ~ x1+x2 + x1:I(x2<0), data=d, DZ="DZ", zyg="zyg", id="id", type="ace")
##' ace3
##' ## Categorical variables are also supported
##' d2 <- transform(d,x2cat=cut(x2,3,labels=c("Low","Med","High")))
##' ace4 <- twinlm(y ~ x1+x2cat, data=d2, DZ="DZ", zyg="zyg", id="id", type="ace")
##' }
##' @keywords models
##' @keywords regression
##' @param formula Formula specifying effects of covariates on the response
##' @param data \code{data.frame} with one observation pr row. In
##' addition a column with the zygosity (DZ or MZ given as a factor) of
##' each individual much be
##' specified as well as a twin id variable giving a unique pair of
##' numbers/factors to each twin pair
##' @param id The name of the column in the dataset containing the twin-id variable.
##' @param zyg The name of the column in the dataset containing the
##' zygosity variable
##' @param DZ Character defining the level in the zyg variable
##' corresponding to the dyzogitic twins. If this argument is missing,
##' the reference level (i.e. the first level) will be interpreted as
##' the dyzogitic twins
##' @param group Optional. Variable name defining group for interaction analysis (e.g., gender)
##' @param group.equal If TRUE marginals of groups are asummed to be the same
##' @param strata Strata variable name
##' @param weights Weights matrix if needed by the chosen estimator. For use
##' with Inverse Probability Weights
##' @param type Character defining the type of analysis to be
##' performed. Can be a subset of "aced" (additive genetic factors, common
##' environmental factors, unique environmental factors, dominant
##' genetic factors). Other choices are:
##' \itemize{
##' \item {"0" (or "sat"): Saturated model where twin 1 and twin 2 within each twin
##' pair may have a different marginal distribution.}
##' \item {"1" (or "flex","zyg"): Within twin pairs the marginal distribution is
##' the same, but the marginal distribution may differ between MZ and DZ
##' twins. A free correlation structure within MZ and DZ twins. }
##' \item {"2" (or "u", "eqmarg"): All individuals have the same marginals
##' but a free correlation structure within MZ and DZ twins.}
##' }
##' The default value is an additive polygenic model \code{type="ace"}.
##' @param twinnum The name of the column in the dataset numbering the
##' twins (1,2). If it does not exist in \code{data} it will
##' automatically be created.
##' @param binary If \code{TRUE} a liability model is fitted. Note that if the right-hand-side of the formula is a factor, character vector, og logical variable, then the liability model is automatically chosen (wrapper of the \code{bptwin} function).
##' @param ordinal If non-zero (number of bins) a liability model is fitted.
##' @param keep Vector of variables from \code{data} that are not
##' specified in \code{formula}, to be added to data.frame of the SEM
##' @param estimator Choice of estimator/model
##' @param constrain Development argument
##' @param control Control argument parsed on to the optimization routine
##' @param messages Control amount of messages shown
##' @param ... Additional arguments parsed on to lower-level functions
twinlm <- function(formula, data, id, zyg, DZ, group=NULL,
group.equal=FALSE, strata=NULL, weights=NULL, type=c("ace"),
twinnum="twinnum",
binary=FALSE,ordinal=0,
keep=weights,estimator=NULL,
constrain=TRUE,control=list(),messages=1,...)
{
cl <- match.call(expand.dots=TRUE)
opt <- options(na.action="na.pass")
mf <- model.frame(formula,data)
mt <- attr(mf, "terms")
y <- model.response(mf, "any")
## formula <- update(formula, ~ . + 1)
yvar <- getoutcome(formula)
if (missing(zyg)) stop("Zygosity variable not specified")
if (!(zyg%in%colnames(data))) stop("'zyg' not found in data")
if (!(id%in%colnames(data))) stop("'id' not found in data")
if (missing(id)) stop("Twin-pair variable not specified")
if (binary || ordinal || is.factor(mf[,yvar]) || is.character(mf[,yvar]) || is.logical(mf[,yvar])) {
if (length(unique(mf[,yvar]))>2 || ordinal) {
if (is.character(mf[,yvar])) {
y <- as.factor(y)
}
if (ordinal<2) ordinal <- length(unique(mf[,yvar]))
y <- as.numeric(y)-1
} else {
args <- as.list(cl)
names(args)[which(names(args)=="formula")] <- "x"
args[[1]] <- NULL
return(do.call("bptwin",args,envir=parent.frame()))
}
}
cl$ordinal <- ordinal
formulaId <- unlist(Specials(formula,"cluster"))
formulaStrata <- unlist(Specials(formula,"strata"))
formulaSt <- paste("~.-cluster(",formulaId,")-strata(",paste(formulaStrata,collapse="+"),")")
formula <- update(formula,formulaSt)
if (!is.null(formulaId)) {
id <- formulaId
cl$id <- id
}
if (!is.null(formulaStrata)) strata <- formulaStrata
cl$formula <- formula
if (!is.null(strata)) {
dd <- split(data,interaction(data[,strata]))
nn <- unlist(lapply(dd,nrow))
dd[which(nn==0)] <- NULL
if (length(dd)>1) {
fit <- lapply(seq(length(dd)),function(i) {
if (messages>0) message("Strata '",names(dd)[i],"'")
cl$data <- dd[[i]]
eval(cl)
})
res <- list(model=fit)
res$strata <- names(res$model) <- names(dd)
class(res) <- c("twinlm.strata","twinlm")
res$coef <- unlist(lapply(res$model,coef))
res$vcov <- blockdiag(lapply(res$model,vcov))
res$N <- length(dd)
res$idx <- seq(length(coef(res$model[[1]])))
rownames(res$vcov) <- colnames(res$vcov) <- names(res$coef)
return(res)
}
}
type <- tolower(type)
if (type%in%c("u", "samemarg", "eqmarg", "2"))
type <- c("u")
if (type%in%c("flex", "zyg", "1"))
type <- c("flex")
if (type%in%c("sat", "0")) {
type <- c("sat")
}
varnames <- all.vars(formula)
latentnames <- c("a1","a2","c1","c2","d1","d2","e1","e2")
if (any(latentnames%in%varnames))
stop(paste(paste(latentnames,collapse=",")," reserved for names of latent variables.",sep=""))
M <- model.matrix(formula,mf)
options(opt)
covars <- colnames(M)
hasIntercept <- FALSE
if (attr(terms(formula),"intercept")==1) {
hasIntercept <- TRUE
covars <- covars[-1]
}
if(length(covars)<1) covars <- NULL
zygstat <- data[,zyg]
if(!is.factor(zygstat)) {
zygstat <- as.factor(zygstat)
}
zyglev <- levels(zygstat)
if (length(zyglev)>2) stop("More than two zygosity levels found. For opposite sex (OS) analysis use the 'group' argument (and regroup OS group as DZ, e.g. DZ=c('OS','DZ'))")
if (tolower(type)=="cor") type <- "u"
if (!is.null(group) && type%in%c("u","flex","sat")) stop("Only polygenic models are allowed with 'group' ('type' subset of 'acde'). See also the 'strata' argument.")
## To wide format:
num <- NULL; if (twinnum%in%colnames(data)) num <- twinnum
if (!is.null(group)) data[,group] <- as.factor(data[,group])
y <- data.frame(y); names(y) <- yvar
#y <- cbind(as.vector(y)); colnames(y) <- yvar
data <- cbind(y,data[,c(keep,num,zyg,id,group)],M)
ddd <- fast.reshape(data,id=c(id),varying=c(yvar,keep,covars,group),keep=zyg,num=num,sep=".",labelnum=TRUE)
groups <- paste(group,".",1:2,sep="")
outcomes <- paste(yvar,".",1:2,sep="")
if (missing(DZ)) {
warning("Using first level, `",zyglev[1],"', in status variable as indicator for 'dizygotic'", sep="")
DZ <- zyglev[1]
}
MZ <- setdiff(zyglev,DZ)
wide1 <- ddd[which(ddd[,zyg]==MZ),,drop=FALSE]
wide2 <- ddd[which(ddd[,zyg]%in%DZ),,drop=FALSE]
mm <- nn <- c()
dd <- list()
levgrp <- NULL
if (!is.null(group)) {
levgrp <- levels(data[,group])
for (i1 in levgrp) {
for (i2 in levgrp) {
idxMZ <- which(wide1[,groups[1]]==i1 & wide1[,groups[2]]==i2)
dMZ <- wide1[idxMZ,,drop=FALSE]
idxDZ <- which(wide2[,groups[1]]==i1 & wide2[,groups[2]]==i2)
dDZ <- wide2[idxDZ,,drop=FALSE]
m0 <- twinsem1(outcomes,c(i1,i2),
levels=levgrp,covars=covars,type=type,
data=list(dMZ,dDZ),constrain=constrain,
equal.marg=group.equal,intercept=hasIntercept,
ordinal=ordinal)$model
if (length(idxMZ)>0) {
nn <- c(nn,paste("MZ:",i1,sep=""))
dd <- c(dd,list(dMZ))
mm <- c(mm,list(m0[[1]]))
}
if (length(idxDZ)>0) {
nn <- c(nn,paste("DZ:",i1," ",i2,sep=""))
dd <- c(dd,list(dDZ))
mm <- c(mm,list(m0[[2]]))
}
}
}; names(mm) <- nn; names(dd) <- nn
} else {
mm <- twinsem1(outcomes,NULL,
levels=NULL,covars=covars,type=type,
data=list(wide1,wide2),constrain=constrain,
intercept=hasIntercept,ordinal=ordinal)$model
dd <- list(MZ=wide1,DZ=wide2)
}
newkeep <- unlist(sapply(keep, function(x) paste(x, 1:2, sep = ".")))
if (!is.null(estimator) && (inherits(estimator,c("numeric","logical")) && !estimator)) return(multigroup(mm, dd, missing=TRUE,fix=FALSE,keep=newkeep,type=2))
optim <- list()
if (is.null(control$start)) {
optim <- list(start=rep(0.1,length(coef(mm[[1]]))*length(mm)))
## optim <- list(method="nlminb2",refit=FALSE,gamma=1,start=rep(0.1,length(coef(mm[[1]]))*length(mm)))
optim$start <- twinlmStart(formula,na.omit(mf),type,hasIntercept,
surv=inherits(data[,yvar],"Surv"),ordinal=ordinal,
model=mm, group=levgrp, group.equal=group.equal)
}
if (length(control)>0) {
optim[names(control)] <- control
}
if (inherits(data[,yvar],"Surv")) {
##if (!requireNamespace("lava.tobit",quietly=TRUE)) stop("lava.tobit required")
if (!is.null(estimator) && estimator%in%c("gaussian","tobit")) optim$method <- "nlminb1"
suppressWarnings(e <- estimate(mm,dd,control=optim,estimator=estimator,...))
} else {
suppressWarnings(e <- estimate(mm,dd,weights=weights,estimator=estimator,fix=FALSE,control=optim,...))
}
if (!is.null(optim$refit) && optim$refit) {
optim$method <- "NR"
optim$start <- pars(e)
if (inherits(mf[,yvar],"Surv")) {
suppressWarnings(e <- estimate(mm,dd,estimator=estimator,fix=FALSE,control=optim,...))
} else {
suppressWarnings(e <- estimate(mm,dd,weights=weights,estimator=estimator,fix=FALSE,control=optim,...))
}
}
e$vcov <- Inverse(information(e,type="hessian"))
counts <- function(dd) {
dd0 <- apply(dd,2,function(x) !is.na(x))
pairs <- sum(dd0[,1]*dd0[,2])
singletons <- sum((!dd0[,1])*dd0[,2] + (!dd0[,2])*dd0[,1])
return(c(pairs,singletons))
}
counts <- lapply(dd, function(x) counts(x[,outcomes]))
## mz <- counts(object$data.mz[,object$outcomes])
## dz <- counts(object$data.dz[,object$outcomes])
if (!e$model$missing) {
zygtab <- c("MZ-pairs"=counts[[1]][1],"DZ-pairs"=counts[[2]][1])
} else {
zygtab <- c(paste(counts[[1]],collapse="/"),paste(counts[[2]],collapse="/"))
names(zygtab) <- c("MZ-pairs/singletons","DZ-pairs/singletons")
}
res <- list(coefficients=e$opt$estimate, vcov=e$vcov,
estimate=e, model=mm, call=cl, data=data, zyg=zyg,
id=id, twinnum=twinnum, type=type, group=group,
constrain=constrain, outcomes=outcomes, zygtab=zygtab,
nam=nn, groups=levgrp, group.equal=group.equal,
counts=counts,ordinal=ordinal)
class(res) <- "twinlm"
return(res)
}
###}}} twinlm
###{{{ twinsem1 (create lava model)
##outcomes <- c("y1","y2"); groups <- c("M","F"); covars <- NULL; type <- "ace"
##twinsem1(c("y1","y2"),c("M","F"))
twinsem1 <- function(outcomes,groups=NULL,levels=NULL,covars=NULL,type="ace",data,constrain=TRUE,equal.marg=FALSE,intercept=TRUE,ordinal=0,...) {
isA <- length(grep("a",type))>0 & type!="sat"
isC <- length(grep("c",type))>0
isD <- length(grep("d",type))>0
isE <- length(grep("e",type))>0 | type=="sat" | type=="u"
lambdas <- c("lambda[a]","lambda[c]","lambda[d]","lambda[e]")
varidx <- which(c(isA,isC,isD,isE))
vMZ1 <- c("a1","c1","d1","e1")
vMZ2 <- c("a1","c1","d1","e2")
vDZ2 <- c("a2","c1","d2","e2")
rhoA <- rhoD <- zA <- zD <- NULL
##if (ordinal>0) {
## mm <- lapply(mm, function(x) {
## x
## })
## }
if (is.list(outcomes)) {
if (!is.null(groups) & is.null(levels)) stop("missing levels")
if (is.null(groups)) groups <- c("","")
grp <- paste(sort(groups),collapse=" ")
sameGroup <- groups[1]==groups[2]
model1 <- outcomes[[1]]
model2 <- outcomes[[2]]
if (!is.null(levels)) {
pars <- c()
for (i in seq(length(levels)-1)) for (j in seq(i+1,length(levels)))
pars <- c(pars,paste(sort(levels)[c(i,j)],collapse=" "))
if (isA) {
parameter(model1) <- paste("z(A):",pars,sep="")
parameter(model2) <- paste("z(A):",pars,sep="")
}
if (isD) {
parameter(model1) <- paste("z(D):",pars,sep="")
parameter(model2) <- paste("z(D):",pars,sep="")
}
}
outcomes <- endogenous(model1)
if (!(type%in%c("u","flex","sat"))) {
if (equal.marg) {
regression(model1,to=outcomes[1],from=vMZ1[varidx]) <-
lambdas[varidx]
regression(model1,to=outcomes[2],from=vMZ2[varidx]) <-
lambdas[varidx]
regression(model2,to=outcomes[1],from=vMZ1[varidx]) <-
lambdas[varidx]
regression(model2,to=outcomes[2],from=vDZ2[varidx]) <-
lambdas[varidx]
} else {
regression(model1,to=outcomes[1],from=vMZ1[varidx]) <-
paste(lambdas,groups[1],sep="")[varidx]
regression(model1,to=outcomes[2],from=vMZ2[varidx]) <-
paste(lambdas,groups[2],sep="")[varidx]
regression(model2,to=outcomes[1],from=vMZ1[varidx]) <-
paste(lambdas,groups[1],sep="")[varidx]
regression(model2,to=outcomes[2],from=vDZ2[varidx]) <-
paste(lambdas,groups[2],sep="")[varidx]
}
}
if (sameGroup) {
if (isA) covariance(model2,a1~a2) <- 0.5
if (isD) covariance(model2,d1~d2) <- 0.25
} else {
if (isA) {
rhoA <- paste("Kinship[A]:",grp,sep="")
zA <- paste("z(A):",grp,sep="")
covariance(model2, a1~a2) <- rhoA
constrain(model2, rhoA,zA) <- function(x) tanh(x)
}
if (isD) {
rhoD <- paste("Kinship[D]:",grp,sep="")
zD <- paste("z(D):",grp,sep="")
covariance(model2, d1~d2) <- rhoD
constrain(model2, rhoD,zD) <- function(x) tanh(x)
}
}
if (ordinal>0) {
for (i in c("model1","model2")) {
model <- get(i)
yy <- lava::endogenous(model)
if (type=="sat") {
lava::ordinal(model,K=ordinal) <- yy
} else {
pname <- paste0("t",seq_len(ordinal-1))
if (type=="flex") pname <- paste0(pname,gsub("model","",i))
lava::ordinal(model,K=ordinal,pname) <- yy
}
if (!(type%in%c("u","flex","sat"))) {
lava::covariance(model,yy) <- 0
lava::regression(model,from="e1",to=yy[1]) <- 1
lava::regression(model,from="e2",to=yy[2]) <- 1
}
assign(i,model)
}
}
return(list(model=list(MZ=model1,DZ=model2),
zA=zA, rhoA=rhoA, zD=zD, rhoD=rhoD))
}
### Build model from scratch....
model1<- lvm()
if (!(type%in%c("u","flex","sat"))) {
model1 <- regression(model1,to=outcomes[1],from=vMZ1[varidx])
model1 <- regression(model1,to=outcomes[2],from=vMZ2[varidx])
} else {
addvar(model1) <- outcomes
covariance(model1) <- outcomes
}
latent(model1) <- union(vMZ1[varidx],vMZ2[varidx])
intercept(model1,latent(model1)) <- 0
if (!is.null(covars))
for (i in 1:length(covars)) {
regression(model1, from=paste(covars[i],".1",sep=""), to=outcomes[1],silent=TRUE) <- paste("beta[",i,"]",sep="")
regression(model1, from=paste(covars[i],".2",sep=""), to=outcomes[2],silent=TRUE) <- paste("beta[",i,"]",sep="")
}
covariance(model1,outcomes) <- 0
covariance(model1, latent(model1)) <- 1
if (!type%in%c("sat","flex")) {
intercept(model1,outcomes) <- "mu"
}
if (type%in%c("u","flex","sat")) {
kill(model1) <- ~e1+e2
covariance(model1,outcomes) <- "v1"
}
model2 <- model1
if (!(type%in%c("u","flex","sat"))) {
model2 <- cancel(model2,c(outcomes[2],vMZ2[varidx]))
model2 <- regression(model2,to=outcomes[2],from=vDZ2[varidx])
latent(model2) <- vDZ2[varidx]
intercept(model2, latent(model2)) <- 0
covariance(model2, latent(model2)) <- 1
}
if (type=="flex") {
varMZ <- paste("var(MZ)",groups,sep="")
varDZ <- paste("var(DZ)",groups,sep="")
intercept(model1,outcomes) <- "mu1"
intercept(model2,outcomes) <- "mu2"
covariance(model1,outcomes) <- varMZ
covariance(model2,outcomes) <- varDZ
}
if (type=="sat") {
varMZ <- c(paste("var(MZ)1",groups[1],""),
paste("var(MZ)2",groups[2],""))
varDZ <- c(paste("var(DZ)1",groups[1],""),
paste("var(DZ)2",groups[2],""))
covariance(model1,outcomes) <- varMZ
covariance(model2,outcomes) <- varDZ
}
if (type%in%c("u","flex","sat")) {
if (constrain) {
if (ordinal) {
covariance(model1,outcomes,pairwise=TRUE) <- "covMZ"
covariance(model2,outcomes,pairwise=TRUE) <- "covDZ"
constrain(model1,"atanh(rhoMZ)","covMZ") <- tanh
constrain(model2,"atanh(rhoDZ)","covDZ") <- tanh
}
if (type=="sat") {
if (!ordinal) {
model1 <- covariance(model1,outcomes,constrain=TRUE,rname="atanh(rhoMZ)",cname="covMZ",lname="log(var(MZ)).1",l2name="log(var(MZ)).2")
model2 <- covariance(model2,outcomes,constrain=TRUE,rname="atanh(rhoDZ)",cname="covDZ",lname="log(var(DZ)).1",l2name="log(var(DZ)).2")
}
} else {
if (type=="flex") {
if (!ordinal) {
model1 <- covariance(model1,outcomes,constrain=TRUE,rname="atanh(rhoMZ)",cname="covMZ",lname="log(var(MZ))")
model2 <- covariance(model2,outcomes,constrain=TRUE,rname="atanh(rhoDZ)",cname="covDZ",lname="log(var(DZ))")
}
} else {
if (!ordinal) {
model1 <- covariance(model1,outcomes,constrain=TRUE,rname="atanh(rhoMZ)",cname="covMZ",lname="log(var)")
model2 <- covariance(model2,outcomes,constrain=TRUE,rname="atanh(rhoDZ)",cname="covDZ",lname="log(var)")
}
}
}
} else {
covariance(model1,outcomes[1],outcomes[2]) <- "covMZ"
covariance(model2,outcomes[1],outcomes[2]) <- "covDZ"
}
}
if (!is.null(covars) & type%in%c("flex","sat")) {
sta <- ""
if (type=="sat") sta <- "b"
for (i in 1:length(covars)) {
regression(model1, from=paste(covars[i],".1",sep=""), to=outcomes[1],silent=TRUE) <- paste("beta1[",i,"]",sep="")
regression(model1, from=paste(covars[i],".2",sep=""), to=outcomes[2],silent=TRUE) <- paste("beta1",sta,"[",i,"]",sep="")
regression(model2, from=paste(covars[i],".1",sep=""), to=outcomes[1],silent=TRUE) <- paste("beta2[",i,"]",sep="")
regression(model2, from=paste(covars[i],".2",sep=""), to=outcomes[2],silent=TRUE) <- paste("beta2",sta,"[",i,"]",sep="")
}
}
if (!intercept) {
intercept(model1,outcomes) <- 0
intercept(model2,outcomes) <- 0
}
## Full rank covariate/design matrix?
if (!missing(data))
for (i in covars) {
myvars <- paste(i,c(1,2),sep=".")
dif <- data[[1]][,myvars[1]]-data[[1]][,myvars[2]]
mykeep <- myvars
if (all(na.omit(dif)==00)) {
mykeep <- mykeep[-2]
}
trash <- setdiff(myvars,mykeep)
if (length(mykeep)==1) {
regression(model1, to=outcomes[2], from=mykeep) <- lava::regfix(model1)$label[trash,outcomes[2]]
kill(model1) <- trash
}
dif <- data[[2]][,myvars[1]]-data[[2]][,myvars[2]]
mykeep <- myvars
if (all(na.omit(dif)==00)) {
mykeep <- mykeep[-2]
}
trash <- setdiff(myvars,mykeep)
if (length(mykeep)==1) {
regression(model2, to=outcomes[2], from=mykeep) <- lava::regfix(model2)$label[trash,outcomes[2]]
kill(model2) <- trash
}
}
twinsem1(list(MZ=model1,DZ=model2),groups=groups,type=type,levels=levels,constrain=constrain,equal.marg=equal.marg,ordinal=ordinal)
}
###}}} twinsem1
###{{{ twinlmStart (starting values)
twinlmStart <- function(formula,mf,type,hasIntercept,surv=FALSE,ordinal=0,model,group=NULL,group.equal,...) {
if (surv) {
l <- survival::survreg(formula,data=mf,dist="gaussian")
beta <- coef(l)
sigma <- l$scale
} else {
if (ordinal) {
l <- lava::ordreg(formula,mf)
beta <- coef(l)
sigma <- 1
} else {
l <- lm.fit(model.matrix(formula,mf),mf[,1])
beta <- l$coefficients
sigma <- sd(l$residuals)
}
}
intidx <- 1
if (ordinal) intidx <- seq_len(ordinal-1)
start <- rep(sigma/sqrt(nchar(type)),nchar(type))
if (ordinal) start <- start[-length(start)]
if (hasIntercept) {
start <- c(beta[intidx],start)
start <- c(start,beta[-intidx])
} else start <- c(start,beta)
varp <- 0.5
if (type=="sat") {
if (!ordinal) varp <- c(rep(log(sigma^2),2),varp)
start <- c()
if (hasIntercept) {
start <- rep(beta[intidx],4)
beta <- beta[-intidx]
}
start <- c(start,rep(c(rep(beta,2),varp),2))
}
if (type=="flex") {
if (!ordinal) varp <- c(log(sigma^2),varp)
start <- c()
if (hasIntercept) {
start <- c(beta[intidx],beta[intidx])
beta <- beta[-intidx]
}
start <- c(start,rep(c(beta,varp),2))
}
if (type=="u") {
start <- c(varp,varp)
if (!ordinal) start <- c(log(sigma^2),start)
start <- c(beta,start)
}
names(start) <- NULL
if (!is.null(group)) {
cc <- coef(model[[1]])
iA <- grep(c("z(A):"),cc,fixed=TRUE)
iD <- grep(c("z(D):"),cc,fixed=TRUE)
nplus <- length(iA) + length(iD) + max(length(iA),length(iD))
start <- c(start,rep(.2,nplus))
ii <- sort(c(iA,iD))
start <- c(start[seq(ii[1]-1)],rep(0.3,length(ii)),start[seq(ii[1]+1,length(start))])
}
return(start)
}
###}}} twinlmStart
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.