Nothing
predict.class <- function(results) {
### results should be an object of class inarmix
a <- apply(results$post.probs,1,function(x) which(x==max(x))[1])
return(a)
}
EnforceConstraints <- function(params,Xlist,XX,subj.factor) {
beta <- as.vector(params$beta)
ff <- function(X) { return(X %*% beta)}
gg <- function(x) {
tmp <- abs(diff(x))
ans <- ifelse(length(tmp) > 0,max(tmp),0)
return(ans)
}
lin.stack <- XX%*%beta
split.lin <- split(lin.stack,subj.factor)
bounds <- unlist(lapply(split.lin,gg))
bound <- max(bounds)
alpha.bound <- exp(-bound/2) + 1e-4
newalpha <- min(params$alpha,alpha.bound)
alpha <- max(newalpha,1e-4)
alpha <- min(alpha,1 - 1e-4)
gamma <- max(params$gamma,1e-4)
return(list(beta=beta,alpha=alpha,gamma=gamma))
}
MyPoisson_Reg <- function(yy,XX,numiter=5) {
### This function computes regression coefficients for
### a Poisson regression under a log-link function.
d <- ncol(XX)
beta <- rep(1,d)
### Perform 5 iterations of the Fisher
### scoring algorithm
for(k in 1:numiter) {
mu <- as.vector(exp(XX%*%beta))
W <- diag(mu)
A <- crossprod(W%*%XX,XX)
### Check that the condition number of the
### above matrix is not extremely large
if(rcond(A) < sqrt(.Machine$double.eps)) {
A <- A + diag(rep(2e-8,d))
}
### update beta
beta <- beta + solve(A,crossprod(XX,yy-mu))
}
return(beta)
}
InitializePars = function(formula,data,subj.col,nclasses,nvec,yy,XX,ylist,Xlist) {
beta <- matrix(0,nrow=nclasses,ncol=ncol(XX))
alpha <- gamma <- mix.prop <- rep(1/nclasses,nclasses)
deviance <- matrix(0,nrow=length(yy),ncol=nclasses)
m <- length(nvec)
fitted <- numeric(length(yy))
subject.ind <- rep(1:m,nvec)
#### Sample nclasses number of people
subject.init <- sample(1:m,size=nclasses)
for(k in 1:nclasses) {
#tmpdata <- data.frame(y=ylist[[subject.init[k]]],Xlist[[subject.init[k]]])
current_y <- ylist[[subject.init[k]]]
current_X <- Xlist[[subject.init[k]]]
### run glm
### still need to worry if there are fewer observations
### than predictors.
if(length(current_y) < ncol(XX)) {
ss <- c(1:m)
orig.subj <- subject.init[[k]]
add.subj <- sample(ss[-orig.subj],size=1)
ty <- unlist(ylist[c(orig.subj,add.subj)])
tx <- rbind(Xlist[[orig.subj]],Xlist[[add.subj]])
#tmpdata <- data.frame(y=ty,tx)
if(length(ty) >= ncol(XX)) {
#suppressWarnings(bet <- glm(y ~ . - 1,data=tmpdata,family=poisson(),maxit=5)$coef)
bet <- MyPoisson_Reg(ty,tx)
}
else {
bet <- rep(1,ncol(XX))
}
}
else {
# suppressWarnings(bet <- glm(y ~ . - 1,data=tmpdata,family=poisson(),maxit=5)$coef)
bet <- MyPoisson_Reg(current_y,current_X)
}
fitted <- exp(XX%*%bet)
deviance[,k] <- 2*fitted - 2*yy*log(fitted)
}
tau <- 3
for(i in 1:tau) {
### make assignments based on minimum deviance.
RR <- aggregate(deviance,by=list(subject.ind),FUN="sum")
assign <- apply(RR[,-1],1,function(x) which.min(x))
assign.expand <- rep(assign,nvec)
tt <- data.frame(subject=subject.ind,class=assign.expand)
for(k in 1:nclasses) {
assign.ind <- which(tt$class==k)
if(length(assign.ind)==0) {
assign.ind <- sample(1:m,size=floor(m/nclasses))
}
tmp.yy <- as.vector(yy[assign.ind])
tmp.XX <- as.matrix(XX[assign.ind,])
if(i==tau) {
tmp.nvec <- as.vector(nvec[assign==k])
tmp.m <- length(tmp.nvec)
myfit <- OneclassEsts(yy=tmp.yy,XX=tmp.XX,nvec=tmp.nvec,m=tmp.m)
beta[k,] <- as.numeric(myfit$beta)
alpha[k] <- max(myfit$alpha,.05)
gamma[k] <- max(myfit$gamma,.1)
}
else {
if(length(assign.ind) == 0) {
bet <- bet + rnorm(length(bet))
}
else {
bet <- MyPoisson_Reg(tmp.yy,tmp.XX)
fitted <- exp(XX%*%bet)
deviance[,k] <- 2*fitted - 2*yy*log(fitted)
}
}
}
}
### In case any NA's are present
beta[which(is.na(beta))] <- 1
alpha[which(is.na(alpha))] <- .2
gamma[which(is.na(gamma))] <- 1
tmp.mix <- table(assign)/m
if(length(tmp.mix) < nclasses) {
a <- rep(1,nclasses-length(tmp.mix))
mix.prop <- c(tmp.mix,a)/(sum(tmp.mix) + sum(a))
}
else {
mix.prop <- tmp.mix
}
return(list(beta=beta,alpha=alpha,gamma=gamma,mix.prop=mix.prop))
}
MultiStart <- function(parlist,nclasses,ylist,Xlist) {
m <- length(ylist)
alpha <- gamma <- mix <- numeric(nclasses)
beta.mat <- matrix(0,nrow=nclasses,ncol=ncol(Xlist[[1]]))
aa <- rnorm(nclasses)
subject.ind <- sample(1:m,size=nclasses)
for(k in 1:nclasses) {
#tmpdata <- data.frame(y=ylist[[subject.ind[k]]],Xlist[[subject.ind[k]]])
#suppressWarnings(beta.vals <- glm(y ~ . - 1,data=tmpdata,family="poisson"))
beta.vals <- MyPoisson_Reg(yy=ylist[[subject.ind[k]]],XX=Xlist[[subject.ind[k]]])
resids <- ylist[[subject.ind[k]]] - exp(Xlist[[subject.ind[k]]]%*%beta.vals)
beta.mat[k,] <- beta.vals
gamma[k] <- max(mean(resids^2) - 1,.1)
alpha[k] <- runif(1)
mix[k] <- exp(aa[k])/sum(exp(aa))
}
### Take a convex combination of previous values and new random values
### for the parameters
beta.mat <- .75*parlist$coef + .25*beta.mat
alpha <- .75*parlist$auto.corr + .25*alpha
gamma <- .75*(parlist$scale - 1) + .25*gamma
mix <- .75*parlist$mix.prop + .25*mix
inpars <- list()
inpars$beta <- beta.mat
inpars$alpha <- alpha
inpars$gamma <- gamma
inpars$mix.prop <- mix
return(inpars)
}
PsiFn <- function(par.stack,len.beta,nclasses,yy,XX,ylist,Xlist,nvec,E1,E2,E3,E4,group.remove=NULL) {
beta.vec <- par.stack[1:(len.beta*nclasses)]
beta <- matrix(beta.vec,nrow=nclasses,ncol=len.beta)
alpha.start <- len.beta*nclasses + 1
alpha <- par.stack[alpha.start:(alpha.start + nclasses - 1)]
gam.start <- alpha.start + nclasses
gamma <- par.stack[gam.start:(gam.start + nclasses - 1)]
mix.start <- gam.start + nclasses
tmp.mix <- par.stack[mix.start:(mix.start + nclasses - 2)]
mix.prop <- c(tmp.mix,1 - sum(tmp.mix))
est.eq <- rep(0,length(par.stack))
m <- length(nvec)
betalist <- lapply(1:nrow(beta), function(i) beta[i,])
tmp.pp <- PostProbs(betalist,alpha,gamma,mix.prop,ylist,Xlist,m,nclasses)
post.probs <- tmp.pp$postprobs
N <- sum(nvec)
m <- length(nvec)
A1.store <- Diagonal(N,x=rep(1,N))
A2.store <- Diagonal(N,x=rep(1,N))
##### initialize R.alpha, and A.store
ind1 <- 0
ind2 <- 0
for(k in 1:nclasses) {
ind1 <- ind2 + 1
ind2 <- ind2 + len.beta
weights <- post.probs[k,]
expanded.weights <- rep(weights,nvec)
mu.stack <- as.vector(exp(XX%*%beta[k,]))
diag(A1.store) <- sqrt(expanded.weights*mu.stack)
diag(A2.store) <- (1/sqrt(mu.stack))
R.alpha <- (alpha[k]*E1 + (1 + alpha[k]^2)*E2 + E3)/(1 - alpha[k]^2) + E4
weighted.dat <- sqrt(expanded.weights)*(yy - mu.stack)
est.eq[ind1:ind2] <- as.vector(crossprod(A1.store%*%XX,R.alpha%*%A2.store%*%weighted.dat))
}
alpha.ind <- len.beta*nclasses + 1
gamma.ind <- alpha.ind + nclasses
mix.ind <- seq(gamma.ind + nclasses,gamma.ind + 2*nclasses-2,by=1)
for(k in 1:nclasses) {
weights <- post.probs[k,]
expanded.weights <- rep(weights,nvec)
mu.stack <- as.vector(exp(XX%*%beta[k,]))
diag(A1.store) <- (1/sqrt(mu.stack))
weighted.dat <- sqrt(expanded.weights)*(yy - mu.stack)
#### alpha est.eq
#### First update R.alpha as R.sandwich
N.weight <- sum(weights*(nvec - 1))
R.alpha <- -((1 + alpha[k]^2)*E1 + 4*alpha[k]*E2 + 2*alpha[k]*E3)/((1 - alpha[k]^2)^2)
LHT <- as.numeric(crossprod(A1.store%*%weighted.dat,R.alpha%*%A1.store%*%weighted.dat))
RHT <- (2*alpha[k]*N.weight*(1+gamma[k]))/(1 - alpha[k]^2)
est.eq[alpha.ind] <- LHT + RHT
#### gamma est.eq
#### update R.alpha
N.weight <- sum(weights*nvec)
R.alpha <- (alpha[k]*E1 + (1 + alpha[k]^2)*E2 + E3)/(1 - alpha[k]^2) + E4
comp2 <- (1/(1 + gamma[k]))*as.numeric(crossprod(A1.store%*%weighted.dat,R.alpha%*%A1.store%*%weighted.dat))
est.eq[gamma.ind] <- comp2 - N.weight
gamma.ind <- gamma.ind + 1
alpha.ind <- alpha.ind + 1
}
est.eq[mix.ind] <- rowSums(post.probs)[1:(nclasses-1)] - m*tmp.mix
return(est.eq/m)
}
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.