Nothing
R/SpeciesMix.R
In SpeciesMix: Fit Mixtures of Archetype Species
# This is package SpeciesMix
"additive.logistic" <-
function (x,inv=FALSE)
{
if(inv){
x <- log(x/x[length(x)])
return(x)
}
x.t <- exp(x)
x.t <- x.t/(1+sum(x.t))
x.t[length(x.t)+1] <- 1-sum(x.t)
return(x.t)
}
"apply.glm" <-
function (i,form,datsp,tau,n)
{
dat.tau <- rep(tau[,i],each=n)
x <- model.matrix(as.formula(form),data=datsp)
y <- datsp$obs
##dat.tau <- get("dat.tau")
##datsp <- get("datsp")
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
f.mix <- glm.fit(x,y,dat.tau,family=binomial())
return(list(coef=f.mix$coef))
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="binomial",x=T,y=T)
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
}
"apply.glm.gaussian" <-
function (i,form,datsp,tau,n)
{
dat.tau <- rep(tau[,i],each=n)
x <- model.matrix(as.formula(form),data=datsp)
y <- datsp$obs
environment(form) <- environment()
##dat.tau <- get("dat.tau")
##datsp <- get("datsp")
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
##f.mix <- glm.nb(as.formula(form),data=datsp,weights=dat.tau)
f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="gaussian")
##f.mix <- glm.nbinom(as.formula(form),data=datsp,weights=dat.tau)
sp.int <- rep(f.mix$coef[1],dim(tau)[1])
##return(list(coef=f.mix$coef[-1],theta=f.mix$theta))
return(list(coef=f.mix$coef[-1],theta=sqrt(f.mix$deviance/f.mix$df.residual),sp.intercept=sp.int))
##return(list(coef=f.mix$coef,theta=f.mix$theta))
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="binomial",x=T,y=T)
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
}
"apply.glm.nbinom" <-
function (i,form,datsp,tau,n)
{
dat.tau <- rep(tau[,i],each=n)
x <- model.matrix(as.formula(form),data=datsp)
y <- datsp$obs
environment(form) <- environment()
##dat.tau <- get("dat.tau")
##datsp <- get("datsp")
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
##f.mix <- glm.nb(as.formula(form),data=datsp,weights=dat.tau)
#f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="poisson")
f.mix <- glm.nbinom(as.formula(form),data=datsp,weights=dat.tau)
sp.int <- rep(f.mix$coef[1],dim(tau)[1])
##return(list(coef=f.mix$coef[-1],theta=f.mix$theta))
return(list(coef=f.mix$coef[-1],theta=f.mix$theta,sp.intercept=sp.int))
##return(list(coef=f.mix$coef,theta=f.mix$theta))
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="binomial",x=T,y=T)
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
}
"apply.glm.tweedie" <-
function (i,form,datsp,tau,n)
{
dat.tau <- rep(tau[,i],each=n)
x <- model.matrix(as.formula(form),data=datsp)
y <- datsp$obs
environment(form) <- environment()
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
f.mix <- tglm(as.formula(form),data=datsp,wts=dat.tau,vcov=FALSE,residuals=FALSE,trace=0)
sp.int <- rep(f.mix$coef[1],dim(tau)[1])
return(list(coef=f.mix$coef[c(-1,-(length(f.mix$coef)-1),-(length(f.mix$coef)))],phi=f.mix$coef["phi"],p=f.mix$coef["p"],sp.intercept=sp.int))
}
"artificial.data" <-
function (formula, data, theta, S, dist = "bernoulli")
{
X <- model.matrix(formula, data)
out <- matrix(0, dim(X)[1], S)
k <- dim(theta)[1]
sp.int <- rep(0, S)
group <- rep(0, S)
for (s in 1:S) {
g <- ceiling(runif(1) * k)
if (dist == "bernoulli") {
lgtp <- X %*% theta[g, ]
p <- exp(lgtp)/(1 + exp(lgtp))
out[, s] <- rbinom(dim(X)[1], 1, p)
}
if (dist == "negbin") {
tmp <- rep(1e+05, dim(X)[1])
while (max(tmp, na.rm = T) > 5000 | sum(tmp) < 100) {
theta[g, 1] <- runif(1, -15, 5)
sp.int[s] <- theta[g, 1]
lgtp <- X %*% theta[g, ]
p <- exp(lgtp)
tmp <- rnbinom(dim(X)[1], mu = p, size = 1)
}
out[, s] <- tmp
}
if (dist == "tweedie") {
tmp <- rep(6e+05, dim(X)[1])
while (max(tmp, na.rm = T) > 5e+05 | sum(tmp) < 100) {
theta[g, 1] <- runif(1, -15, 5)
theta[g, 1] <- runif(1, 1, 5)
sp.int[s] <- theta[g, 1]
lgtp <- X %*% theta[g, ]
p <- exp(lgtp)
tmp <- rTweedie(dim(X)[1], mu = p, phi = 2, p = 1.6)
}
out[, s] <- tmp
}
if (dist == "gaussian") {
tmp <- rep(1e+05, dim(X)[1])
while (max(tmp, na.rm = T) > 50000 | sum(tmp) < 100) {
theta[g, 1] <- runif(1, 100, 500)
lgtp <- X %*% theta[g, ]
p <- (lgtp)
tmp <- rnorm(dim(X)[1], mean = p, sd = 1)
}
out[, s] <- tmp
}
group[s] <- g
}
pi <- tapply(group, group, length)/S
if (dist == "negbin")
return(list(pa = out, group = group, pi = pi, sp.int = sp.int))
if (dist == "tweedie")
return(list(pa = out, group = group, pi = pi, sp.int = sp.int))
list(pa = out, group = group, pi = pi)
}
"clusterSelect" <-
function (sp.form,sp.data,covar.data,G=1:10,em.prefit=TRUE, em.steps=4 ,em.refit=3,est.var=FALSE,trace=TRUE)
{
my.fun <- function(g,form,sp.data,covar.data){
cat("Fitting group",g,"\n")
try(SpeciesMix(form,sp.data,covar.data,g,em.prefit=em.prefit,em.steps=em.steps,em.refit=em.refit,est.var=est.var,trace=trace))
}
# if(mc){ out <- mclapply(G,my.fun,form,dat,mc.preschedule = FALSE, mc.set.seed = TRUE, mc.silent = FALSE, mc.cores = set.cores)} else
{ out <- lapply(G,my.fun,sp.form,sp.data,covar.data)}
aic <- rep(0,length(G))
bic <- rep(0,length(G))
fm <- list()
for(i in 1:length(G))
if(!is.atomic(out[[i]])){
aic[i] <- out[[i]]$aic
bic[i] <- out[[i]]$bic
fm[[i]] <- list(logl=out[[i]]$logl,coef=out[[i]]$coef,tau=out[[i]]$tau,pi=out[[i]]$pi,covar=out[[i]]$covar)
}
return(list(aic=aic,bic=bic,fm=fm))
}
"create.starting.values" <-
function (S,G,n,form,datsp)
{
##apply.glm <- function(i,form,datsp,tau,n,dat.tau){
## dat.tau <- rep(tau[,i],each=n)
## dat.tau <- get("dat.tau")
## datsp <- get("datsp")
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
## f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="binomial",x=T,y=T)
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
## list(coef=f.mix$coef)
##}
environment(form) <- environment()
tau <- matrix(runif(S*G),S,G)
tau <- (tau/rowSums(tau))
fmM <- list()
for(i in 1:G){
##dat.tau[,i] <- rep(tau[,i],each=n)
pi[i] <- sum(tau[,i])/S
}
##dat.tau <- rep(0,dim(datsp)[1])
fmM <- lapply(1:G,apply.glm,form,datsp,tau,n)
first.fit <- list(x=model.matrix(as.formula(form),data=datsp),y=datsp$obs,formula=form)
##return(list(pi=pi,fmM=fmM,tau=tau,dat.tau=dat.tau,first.fit=first.fit))
return(list(pi=pi,fmM=fmM,tau=tau,first.fit=first.fit))
}
"create.starting.values.gaussian" <-
function (S,G,n,form,datsp,mc=FALSE,set.cores=2)
{
##apply.glm <- function(i,form,datsp,tau,n,dat.tau){
## dat.tau <- rep(tau[,i],each=n)
## dat.tau <- get("dat.tau")
## datsp <- get("datsp")
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
## f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="binomial",x=T,y=T)
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
## list(coef=f.mix$coef)
##}
environment(form) <- environment()
tau <- matrix(runif(S*G),S,G)
tau <- (tau/rowSums(tau))
fmM <- list()
for(i in 1:G){
##dat.tau[,i] <- rep(tau[,i],each=n)
pi[i] <- sum(tau[,i])/S
}
##dat.tau <- rep(0,dim(datsp)[1])
fmM <- lapply(1:G,apply.glm.gaussian,form,datsp,tau,n)
first.fit <- list(x=model.matrix(as.formula(form),data=datsp)[,-1],y=datsp$obs,formula=form)
##return(list(pi=pi,fmM=fmM,tau=tau,dat.tau=dat.tau,first.fit=first.fit))
return(list(pi=pi,fmM=fmM,tau=tau,first.fit=first.fit))
}
"create.starting.values.nbinom" <-
function (S,G,n,form,datsp,mc=FALSE,set.cores=2)
{
##apply.glm <- function(i,form,datsp,tau,n,dat.tau){
## dat.tau <- rep(tau[,i],each=n)
## dat.tau <- get("dat.tau")
## datsp <- get("datsp")
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
## f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="binomial",x=T,y=T)
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
## list(coef=f.mix$coef)
##}
environment(form) <- environment()
tau <- matrix(runif(S*G),S,G)
tau <- (tau/rowSums(tau))
fmM <- list()
for(i in 1:G){
##dat.tau[,i] <- rep(tau[,i],each=n)
pi[i] <- sum(tau[,i])/S
}
##dat.tau <- rep(0,dim(datsp)[1])
fmM <- lapply(1:G,apply.glm.nbinom,form,datsp,tau,n)
offset <- model.frame(as.formula(form),data=datsp)
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,length(datsp$obs))
first.fit <- list(x=model.matrix(as.formula(form),data=datsp)[,-1],y=datsp$obs,offset=offset,formula=form)
##return(list(pi=pi,fmM=fmM,tau=tau,dat.tau=dat.tau,first.fit=first.fit))
return(list(pi=pi,fmM=fmM,tau=tau,first.fit=first.fit))
}
"create.starting.values.nbinom.kmeans" <-
function (S, G, n, form, datsp, tol = 0.1)
{
MM <- model.matrix(form, datsp)
offset <- model.frame(form,datsp)
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,nrow(MM))
first.fit <- list(x = model.matrix(as.formula(form), data = datsp)[,
-1], y = datsp$obs, offset=offset, formula = form)
if(class(first.fit$x)=="numeric") {first.fit$x <- matrix(first.fit$x,length(first.fit$x),1) }#deal with model matrix returning a vector for covar == 1
if (tol < 0 || tol >= 1)
stop("Minimum Prevalence % must be between 0 and 1")
sp.name <- 1:S
sp <- rep(sp.name, each = n)
starting.fitem <- list(intercepts = rep(0, S), alpha = rep(0,
S))
all.betas <- matrix(0, nrow = S, ncol = ncol(MM))
colnames(all.betas) <- colnames(MM)
for (j in 1:S) {
fit <- glm.nbinom(form, datsp[((j - 1) * n):(j * n -
1), ], est.var = FALSE)
starting.fitem$sp.intercepts[j] <- fit$coef[1]
all.betas[j, ] <- fit$coef
starting.fitem$theta[j] = fit$theta
}
all.betas[, 1] <- 0
cat("Clustering...\n")
fmmvnorm <- kmeans(x = all.betas, centers = G, iter.max = 100,
nstart = 50)
starting.fitem$coef <- fmmvnorm$centers
fmM <- list()
for (i in 1:G) {
B <- matrix(rep(fmmvnorm$centers[i, ], nrow(datsp)),
nrow(datsp), ncol(fmmvnorm$centers), byrow = T)
B[, 1] <- rep(starting.fitem$sp.intercepts, each = n)
fitted <- exp(rowSums(MM * B)+offset)
fmM[[i]] <- list(coef = fmmvnorm$centers[i, 2:ncol(fmmvnorm$centers)],
theta = mean(starting.fitem$theta[fmmvnorm$cluster ==
i]), sp.intercept = starting.fitem$sp.intercepts,
fitted = fitted)
}
tau <- matrix(0, S, G)
pi <- rep(1/G, G)
pi <- runif(G, 0.2, 0.8)
pi <- pi/sum(pi)
est.tau <- lapply(1:S, estimate.pi.nbinom, sp, sp.name, datsp,
fmM, pi, G, first.fit)
max.newTau <- 0.8
alpha <- (1 - max.newTau * G)/(max.newTau * (2 - G) - 1)
for (j in 1:S) {
newTau <- (2 * alpha * est.tau[[j]]$tau - alpha + 1)/(2 *
alpha - alpha * G + G)
tau[j, ] <- newTau
}
for (i in 1:G) {
pi[i] <- sum(tau[, i])/S
}
return(list(pi = pi, fmM = fmM, tau = tau, first.fit = first.fit))
}
"create.starting.values.tweedie" <-
function (S,G,n,form,datsp)
{
##apply.glm <- function(i,form,datsp,tau,n,dat.tau){
## dat.tau <- rep(tau[,i],each=n)
## dat.tau <- get("dat.tau")
## datsp <- get("datsp")
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
## f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="binomial",x=T,y=T)
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
## list(coef=f.mix$coef)
##}
environment(form) <- environment()
tau <- matrix(runif(S*G),S,G)
tau <- (tau/rowSums(tau))
fmM <- list()
for(i in 1:G){
##dat.tau[,i] <- rep(tau[,i],each=n)
pi[i] <- sum(tau[,i])/S
}
##dat.tau <- rep(0,dim(datsp)[1])
offset <- model.frame(form,datsp)
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,length(datsp$obs))
fmM <- lapply(1:G,apply.glm.tweedie,form,datsp,tau,n)
first.fit <- list(x=model.matrix(as.formula(form),data=datsp)[,-1],y=datsp$obs,formula=form)
##return(list(pi=pi,fmM=fmM,tau=tau,dat.tau=dat.tau,first.fit=first.fit))
return(list(pi=pi,fmM=fmM,tau=tau,first.fit=first.fit))
}
"create.starting.values.tweedie.kmeans" <-
function (S, G, n, form, datsp, tol = 0.1)
{
MM <- model.matrix(form, datsp)
offset <- model.frame(form,datsp)
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,nrow(MM))
first.fit <- list(x = model.matrix(as.formula(form), data = datsp)[,
-1], y = datsp$obs, offset=offset,formula = form)
if(class(first.fit$x)=="numeric") {first.fit$x <- matrix(first.fit$x,length(first.fit$x),1) }#deal with model matrix returning a vector for covar == 1
if (tol < 0 || tol >= 1)
stop("Minimum Prevalence % must be between 0 and 1")
sp.name <- 1:S
sp <- rep(sp.name, each = n)
starting.fitem <- list(intercepts = rep(0, S), alpha = rep(0,
S))
all.betas <- matrix(0, nrow = S, ncol = ncol(MM))
colnames(all.betas) <- colnames(MM)
for (j in 1:S) {
fit <- tglm(form, datsp[((j - 1) * n):(j * n - 1), ],
vcov = FALSE, residuals = FALSE, trace = 0, p = 1.6)
starting.fitem$sp.intercepts[j] <- fit$coef[1]
all.betas[j, ] <- fit$coef[-length(fit$coef)]
starting.fitem$phi[j] = fit$coef["phi"]
}
all.betas[, 1] <- 0
cat("Clustering...\n")
fmmvnorm <- kmeans(x = all.betas, centers = G, iter.max = 100,
nstart = 50)
starting.fitem$coef <- fmmvnorm$centers
fmM <- list()
for (i in 1:G) {
B <- matrix(rep(fmmvnorm$centers[i, ], nrow(datsp)),
nrow(datsp), ncol(fmmvnorm$centers), byrow = T)
B[, 1] <- rep(starting.fitem$sp.intercepts, each = n)
fitted <- exp(rowSums(MM * B)+offset)
fmM[[i]] <- list(coef = fmmvnorm$centers[i, 2:ncol(fmmvnorm$centers)],
phi = mean(starting.fitem$phi[fmmvnorm$cluster ==
i]), p = 1.6, sp.intercept = starting.fitem$sp.intercepts,
fitted = fitted)
}
tau <- matrix(0, S, G)
pi <- rep(1/G, G)
pi <- runif(G, 0.2, 0.8)
pi <- pi/sum(pi)
est.tau <- lapply(1:S, estimate.pi.tweedie, sp, sp.name,
datsp, fmM, pi, G, first.fit)
max.newTau <- 0.8
alpha <- (1 - max.newTau * G)/(max.newTau * (2 - G) - 1)
for (j in 1:S) {
newTau <- (2 * alpha * est.tau[[j]]$tau - alpha + 1)/(2 *
alpha - alpha * G + G)
tau[j, ] <- newTau
}
for (i in 1:G) {
pi[i] <- sum(tau[, i])/S
}
return(list(pi = pi, fmM = fmM, tau = tau, first.fit = first.fit))
}
"distr.binom" <-
function( p){
nobs <- length( p)
new.dist <- old.dist <- rep( 0, nobs+1)
old.dist[1] <- 1-p[1]
old.dist[2] <- p[1]
for( ii in 2:nobs){
new.dist[1] <- old.dist[1]*(1-p[ii])
for( jj in 2:ii)
new.dist[jj] <- old.dist[jj-1]*p[ii] + old.dist[jj]*(1-p[ii])
new.dist[ii+1] <- old.dist[ii]*p[ii]
old.dist <- new.dist
}
return( new.dist)
}
"dPoisGam" <-
function ( y, lambda, mu.Z, alpha, LOG=TRUE)
{
#function to calculate Random sum (Tweedie) densities.
#y is the value of the r.v. Can be a vector
#mu.N is the mean of the Poisson summing r.v. Can be a vector of length(y)
#mu.Z is the mean of the Gamma rv Can be a vector of length(y)
#alpha is the `other' parameter of the gamma distribution s.t. var = ( mu.Z^2)/alpha Can be a vector of length(y)
#If mu.N, mu.Z or alpha are scalare but y isn't then they will be used for all y. If lengths mis-match then error
#LOG=TRUE gives the density on the log scale
#do.checks=TRUE checks the input vectors for compatability and gives errors / changes them as appropriate.
#do.checks=FALSE doesn't check and relies on the user to have things right. If not right then catastrophic failure may occur.
# if( any( is.null( c( y, mu.N, mu.Z, alpha)))){
# print( "Error: null input values -- please check. Null values are:")
# tmp <- double( is.null( c( y, mu.N, mu.Z, alpha)))
# names( tmp) <- c( "y", "mu.N","mu.Z","alpha")
# print( tmp)
# print( "Exitting")
# return()
# }
mu.N <- lambda
if( !all( is.element( c( length( mu.N), length( mu.Z), length( alpha)), c( length( y), 1)))){
print( "Error: length of parameter vectors does not match length of random variable vector")
print( "Exitting")
return()
}
if( length( mu.N) != length( y))
mu.N <- rep( mu.N, length( y))
if( length( mu.Z) != length( y))
mu.Z <- rep( mu.Z, length( y))
if( length( alpha) != length( y))
alpha <- rep( alpha, length( y))
res <- .Call( "dTweedie", as.numeric( y), as.numeric( mu.N), as.numeric( mu.Z), as.numeric( alpha), as.integer( LOG),PACKAGE="SpeciesMix")
return( res)
}
"dPoisGamDerivs" <-
function ( y=NULL, lambda=NULL, mu.Z=NULL, alpha=NULL, do.checks=TRUE)
{
#function to calculate Random sum (Tweedie) densities.
#y is the value of the r.v. Can be a vector
#mu.N is the mean of the Poisson summing r.v. Can be a vector of length(y)
#mu.Z is the mean of the Gamma rv Can be a vector of length(y)
#alpha is the `other' parameter of the gamma distribution s.t. var = ( mu.Z^2)/alpha Can be a vector of length(y)
#If mu.N, mu.Z or alpha are scalare but y isn't then they will be used for all y. If lengths mis-match then error
#LOG=TRUE gives the density on the log scale
#do.checks=TRUE checks the input vectors for compatability and gives errors / changes them as appropriate.
#do.checks=FALSE doesn't check and relies on the user to have things right. If not right then catastrophic failure may occur.
mu.N <- lambda
if( do.checks){
if( any( is.null( c( y, mu.N, mu.Z, alpha)))){
print( "Error: null input values -- please check. Null values are:")
tmp <- double( is.null( c( y, mu.N, mu.Z, alpha)))
names( tmp) <- c( "y", "mu.N","mu.Z","alpha")
print( tmp)
print( "Exitting")
return()
}
if( !all( is.element( c( length( mu.N), length( mu.Z), length( alpha)), c( length( y), 1)))){
print( "Error: length of parameter vectors does not match length of random variable vector")
print( "Exitting")
}
if( length( mu.N) != length( y))
mu.N <- rep( mu.N, length( y))
if( length( mu.Z) != length( y))
mu.Z <- rep( mu.Z, length( y))
if( length( alpha) != length( y))
alpha <- rep( alpha, length( y))
}
res <- .Call( "dTweedieDeriv", as.numeric( y), as.numeric( mu.N), as.numeric( mu.Z), as.numeric( alpha),PACKAGE="SpeciesMix")
colnames( res) <- c("lambda","mu.Z","alpha")
return( res)
}
"dTweedie" <-
function ( y, mu, phi, p, LOG=TRUE)
{
lambda <- ( mu^( 2-p)) / ( phi*(2-p))
alpha <- ( 2-p) / ( p-1)
tau <- phi*(p-1)*mu^(p-1)
mu.Z <- alpha * tau
dens <- dPoisGam( y, lambda, mu.Z, alpha, LOG)
return( dens)
}
"estimate.pi" <-
function (j,sp,spname,datsp,fmM,pi,G,first.fit)
{
tmp.like <- rep(0,G)
tau <- rep(0,G)
sel.sp <- which(sp==spname[j])
link.fun <- make.link("logit")
for(i in 1:G) {
if(length(fmM[[i]]$coef)==1){lpre <- link.fun$linkinv(first.fit$x[sel.sp,]*fmM[[i]]$coef)
}else{ lpre <- link.fun$linkinv(first.fit$x[sel.sp,]%*%fmM[[i]]$coef)}
##lpre <- link.fun$linkinv(first.fit$x[sel.sp,]%*%fmM[[i]]$coef)
## tmp.like[i] <- sum(dbinom(datsp$obs[sel.sp],1,fmM[[i]]$fitted[sel.sp],log=T))
obs <- datsp$obs[sel.sp]
lpre[obs==0]<- 1- lpre[obs==0]
tmp.like[i] <- sum(log(lpre))
}
eps <- max(tmp.like)
sum.like <- (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
for(i in 1:G) {
tau[i] <- exp((log(pi[i]) + tmp.like[i]) - sum.like )
}
return(list(tau=tau,sum.like=sum.like))
}
"estimate.pi.gaussian" <-
function (j,sp,spname,datsp,fmM,pi,G,first.fit)
{
tmp.like <- rep(0,G)
tau <- rep(0,G)
sel.sp <- which(sp==spname[j])
for(i in 1:G) {
##lpre <- dnorm(first.fit$y[sel.sp],mean=fmM[[i]]$fitted[sel.sp],sd=fmM[[i]]$theta,log=TRUE)
lpre <- dnorm(first.fit$y[sel.sp],mean=fmM[[i]]$fitted[sel.sp],sd=1,log=TRUE)
##lpre <- dnbinom(first.fit$y[sel.sp],mu=exp(cbind(1,first.fit$x[sel.sp,])%*%c(fmM[[i]]$sp.intercept[j],fmM[[i]]$coef)),size=fmM[[i]]$theta,log=TRUE)
##lpre <- dpois(first.fit$y[sel.sp],exp(cbind(1,first.fit$x[sel.sp,])%*%c(fmM[[i]]$sp.intercept[j],fmM[[i]]$coef)),log=TRUE)
#lpre <- dpois(first.fit$y[sel.sp],exp(first.fit$x[sel.sp,]%*%fmM[[i]]$coef),log=TRUE)
tmp.like[i] <- sum(lpre)
}
eps <- max(tmp.like)
sum.like <- (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
for(i in 1:G) {
tau[i] <- exp((log(pi[i]) + tmp.like[i]) - sum.like )
}
return(list(tau=tau,sum.like=sum.like))
}
"estimate.pi.nbinom" <-
function (j,sp,spname,datsp,fmM,pi,G,first.fit)
{
tmp.like <- rep(0,G)
tau <- rep(0,G)
sel.sp <- which(sp==spname[j])
for(i in 1:G) {
lpre <- dnbinom(first.fit$y[sel.sp],mu=fmM[[i]]$fitted[sel.sp],size=fmM[[i]]$theta,log=TRUE)
##lpre <- dnbinom(first.fit$y[sel.sp],mu=exp(cbind(1,first.fit$x[sel.sp,])%*%c(fmM[[i]]$sp.intercept[j],fmM[[i]]$coef)),size=fmM[[i]]$theta,log=TRUE)
##lpre <- dpois(first.fit$y[sel.sp],exp(cbind(1,first.fit$x[sel.sp,])%*%c(fmM[[i]]$sp.intercept[j],fmM[[i]]$coef)),log=TRUE)
#lpre <- dpois(first.fit$y[sel.sp],exp(first.fit$x[sel.sp,]%*%fmM[[i]]$coef),log=TRUE)
tmp.like[i] <- sum(lpre)
}
eps <- max(tmp.like)
sum.like <- (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
for(i in 1:G) {
tau[i] <- exp((log(pi[i]) + tmp.like[i]) - sum.like )
}
return(list(tau=tau,sum.like=sum.like))
}
"estimate.pi.tweedie" <-
function (j, sp, spname, datsp, fmM, pi, G, first.fit)
{
tmp.like <- rep(0, G)
tau <- rep(0, G)
sel.sp <- which(sp == spname[j])
for (i in 1:G) {
lpre <- dTweedie(first.fit$y[sel.sp], mu = fmM[[i]]$fitted[sel.sp],
phi = fmM[[i]]$phi, p = fmM[[i]]$p, LOG = TRUE)
tmp.like[i] <- sum(lpre)
}
eps <- max(tmp.like)
sum.like <- (log(sum(pi * exp((tmp.like) - (eps)))) + (eps))
for (i in 1:G) {
tau[i] <- exp((log(pi[i]) + tmp.like[i]) - sum.like)
}
return(list(tau = tau, sum.like = sum.like))
}
"fitMix" <-
function (form,datsp,sp,G=2,ite.max=500,trace=TRUE,full.model=FALSE,r1=FALSE)
{
## dat2 has colums obs,sp
##
temp.warn <- getOption( "warn")
options( warn=-1)
sp.name <- unique(sp)
S <- length(unique(sp))
n <- length(which(sp==sp.name[1]))
cat("Fitting Group",G,"\n")
if(trace) cat("Iteration | LogL \n")
##dat.tau <- data.frame(matrix(0,dim(datsp)[1],G))
dat.tau <- 0
##dat <- data.frame(datsp,dat.tau)
pi <- rep(0,G)
ite <- 1
logL <- -99999999
old.logL <- -88888888
## set up initial GLM
t1 <- create.starting.values(S,G,n,form,datsp)
pi <- t1$pi;fmM <- t1$fmM;tau <- t1$tau;dat.tau <- t1$dat.tau;first.fit <- t1$first.fit
while(abs(logL-old.logL) > 0.0001 & ite<=ite.max){
old.logL <- logL
for(i in 1:G){
##dat.tau[,i] <- rep(tau[,i],each=n)
pi[i] <- sum(tau[,i])/S
}
if(any(pi==0)) { ## occasionally with complicated models the random starting values result in a pi[i]==0; so restart with new random starts
cat("pi has gone to zero - restarting fitting \n")
t1 <- create.starting.values(S,G,n,form,datsp)
pi <- t1$pi;fmM <- t1$fmM;tau <- t1$tau;dat.tau <- t1$dat.tau;first.fit <- t1$first.fit
ite <- 1
}
fmM <- lapply(1:G,weighted.glm,first.fit,tau,n,fmM,sp)
logL <- 0
tmp.like <- matrix(0,S,G)
est.tau <- lapply(1:S,estimate.pi,sp,sp.name,datsp,fmM,pi,G,first.fit)
for(j in 1:S){
if(is.atomic(est.tau[[j]])){ print (est.tau[[j]])} else
{
tau[j,] <- est.tau[[j]]$tau
logL <- logL+est.tau[[j]]$sum.like
}
}
if(trace) cat(ite," | ",logL,"\n")
ite <- ite+1
}
fm.out <- data.frame(matrix(0,G,length(fmM[[1]]$coef)))
names(fm.out) <- names(fmM[[1]]$coef)
tau <- data.frame(tau)
names(tau) <- paste("grp.",1:G,sep="")
EN <- -sum(unlist(tau)*log(unlist(tau)))
d <- length(unlist(fm.out)) + length(tau)-1
for(i in 1:G) {
fm.out[i,] <- fmM[[i]]$coef
##dat.tau[,i] <- rep(tau[,i],each=n)
}
names(pi) <- paste("G",1:G,sep=".")
t.pi <- additive.logistic(pi,TRUE)
parms <- c(t.pi[1:(G-1)],unlist(fm.out))
logL.full <- logL
logL <- logLmix(parms,first.fit,G,S,sp,sp.name)
options(warn=temp.warn)
if(full.model) return(list(logl=logL,aic = -2*logL + 2*d,tau=round(tau,4),pi=pi,bic=-2*logL + log(S)*d,ICL= -2*logL + log(S)*d +2*EN,coef=fm.out,fmM=fmM,model.tau=dat.tau,covar=0,aic.full= -2*logL.full + 2*d,bic.full= -2*logL.full + log(S)*d))
return(list(logl=logL,aic = -2*logL + 2*d,tau=round(tau,4),pi=pi,bic=-2*logL + log(S)*d,ICL= -2*logL + log(S)*d +2*EN,coef=fm.out,covar=0,aic.full= -2*logL.full + 2*d,bic.full= -2*logL.full + log(S)*d))
}
"fitmix.cpp" <-
function (form,datsp,sp,G=2,pars=NA,trace=TRUE,calc.hes=FALSE,r1)
{
if(!is.numeric(sp)){
sp <- as.integer(factor(sp))
}
sp.name <- unique(sp)
S <- length(unique(sp))
n <- length(which(sp==sp.name[1]))
X <- model.matrix(form, data = datsp[sp==sp.name[1],])
##X <- model.frame(form, data = datsp)
# X <- X[sp==sp.name[1],]
#y <- model.response(form, data = datsp)
y <- datsp$obs
if(is.na(pars[1])) {
##pars <- rep(0.01,G-1+(ncol(X)*G))
pars <- runif(G-1+(ncol(X)*G),-2,2)
pars[1:(G-1)] <- runif(G-1)
}
offset<-rep(0,n)
## hes <- rep(0,length(pars)^2)
gradient <- rep(0,length(pars))
tau <- matrix(0,S,G) ##must leave this in as defines S & G
model.type<- as.integer(1)
loglike <- try(.Call("SpeciesMix",pars,y,X,sp,tau,gradient,offset,model.type,PACKAGE="SpeciesMix"))
calc.deriv <- function(p){
gradient <- rep(0,length(pars))
ll <- .Call("Calculate_Gradient",p,y,X,sp,tau,gradient,offset,as.integer(1),PACKAGE="SpeciesMix")
return(gradient)
}
hes <- 0
if(calc.hes){
hes <- nd2(pars,calc.deriv)
dim(hes) <- rep(length(pars),2)
dim(hes) <- rep(length(pars),2)
rownames(hes) <- colnames(hes) <- c(paste("G.",1:(G-1),sep=""),paste("G",1:G,rep(colnames(X),each=G),sep="."))
}
if(!is.numeric(loglike)) loglike <- 0
pi <- pars[1:(G-1)]
coef <- pars[ (G):length(pars)]
r.logl <- logLmix(pars,list(y=y,x=model.matrix(form, data = datsp)),G,S,sp,sp.name,out.tau=TRUE)
pi <- additive.logistic(pi)
names(pi) <- paste("G.",1:G,sep="")
coef <- matrix(coef,G,ncol(X))
rownames(coef) <- paste("G.",1:G,sep="")
colnames(coef) <- colnames(X)
AIC <- 2*loglike + 2*length(pars)
BIC <- 2*loglike + log(S)*length(pars)
##list(logl=loglike,r.logl=r.logl$logl,pi=pi,coef=coef,tau=round(exp(r.logl$tau),4),aic=AIC,bic=BIC,hessian=hes,gradient=gradient)
list(logl=loglike,pi=pi,coef=coef,tau=round(exp(r.logl$tau),4),aic=AIC,bic=BIC,hessian=hes,gradient=gradient)
}
"fitmix.gaussian" <-
function (form,datsp,sp,G=2,ite.max=500,trace=TRUE,full.model=FALSE)
{
## fitting abundance data with a negative binomial distribution
## dat2 has colums obs,sp
##
temp.warn <- getOption( "warn")
options( warn=-1)
sp.name <- unique(sp)
S <- length(unique(sp))
n <- length(which(sp==sp.name[1]))
cat("Fitting Group",G,"\n")
if(trace) cat("Iteration | LogL \n")
##dat.tau <- data.frame(matrix(0,dim(datsp)[1],G))
dat.tau <- 0
##dat <- data.frame(datsp,dat.tau)
pi <- rep(0,G)
ite <- 1
logL <- -99999999
old.logL <- -88888888
## set up initial GLM
t1 <- create.starting.values.gaussian(S,G,n,form,datsp)
pi <- t1$pi;fmM <- t1$fmM;tau <- t1$tau;dat.tau <- t1$dat.tau;first.fit <- t1$first.fit
while(abs(logL-old.logL) > 0.0001 & ite<=ite.max){
old.logL <- logL
for(i in 1:G){
##dat.tau[,i] <- rep(tau[,i],each=n)
pi[i] <- sum(tau[,i])/S
}
if(any(pi==0)) { ## occasionally with complicated models the random starting values result in a pi[i]==0; so restart with new random starts
cat("pi has gone to zero - restarting fitting \n")
t1 <- create.starting.values.gaussian(S,G,n,form,datsp)
pi <- t1$pi;fmM <- t1$fmM;tau <- t1$tau;dat.tau <- t1$dat.tau;first.fit <- t1$first.fit
ite <- 1
}
fmM <- lapply(1:G,weighted.glm.gaussian,first.fit,tau,n,fmM,sp)
logL <- 0
tmp.like <- matrix(0,S,G)
est.tau <- lapply(1:S,estimate.pi.gaussian,sp,sp.name,datsp,fmM,pi,G,first.fit)
for(j in 1:S){
if(is.atomic(est.tau[[j]])){ print (est.tau[[j]])} else
{
tau[j,] <- est.tau[[j]]$tau
logL <- logL+est.tau[[j]]$sum.like
}
}
if(trace) cat(ite," | ",logL,"\n")
ite <- ite+1
}
fm.out <- data.frame(matrix(0,G,length(fmM[[1]]$coef)))
int.out <- data.frame(matrix(0,S,G))
names(int.out) <- paste("grp.",1:G,sep="")
names(fm.out) <- names(fmM[[1]]$coef)
tau <- data.frame(tau)
names(tau) <- paste("grp.",1:G,sep="")
EN <- -sum(unlist(tau)*log(unlist(tau)))
d <- length(unlist(fm.out)) + length(tau)-1
fm.theta <- rep(0,G)
for(i in 1:G) {
fm.out[i,] <- fmM[[i]]$coef
int.out[,i] <- fmM[[i]]$sp.intercept
##dat.tau[,i] <- rep(tau[,i],each=n)
fm.theta[i] <- fmM[[i]]$theta
}
names(pi) <- paste("G",1:G,sep=".")
t.pi <- additive.logistic(pi,TRUE)
parms <- c(t.pi[1:(G-1)],fm.theta,unlist(fm.out),unlist(int.out))
logL.full <- logL
logL <- logLmix.gaussian(parms,first.fit,G,S,sp,sp.name)
options(warn=temp.warn)
if(full.model) return(list(logl=logL,aic = -2*logL + 2*d,tau=round(tau,4),pi=pi,bic=-2*logL + log(S)*d,ICL= -2*logL + log(S)*d +2*EN,coef=fm.out,sp.intercept=int.out,theta=fm.theta,fmM=fmM,model.tau=dat.tau,covar=0,aic.full= -2*logL.full + 2*d,bic.full= -2*logL.full + log(S)*d))
return(list(logl=logL,aic = -2*logL + 2*d,tau=round(tau,4),pi=pi,bic=-2*logL + log(S)*d,ICL= -2*logL + log(S)*d +2*EN,coef=fm.out,sp.intercept=int.out,theta=fm.theta,covar=0,aic.full= -2*logL.full + 2*d,bic.full= -2*logL.full + log(S)*d))
}
"fitmix.nbinom" <-
function (form, datsp, sp, G = 2, ite.max = 500, trace = TRUE,
full.model = FALSE)
{
temp.warn <- getOption("warn")
options(warn = -1)
sp.name <- unique(sp)
S <- length(unique(sp))
n <- length(which(sp == sp.name[1]))
cat("Fitting Group", G, "\n")
if (trace)
cat("Iteration | LogL \n")
dat.tau <- 0
pi <- rep(0, G)
ite <- 1
logL <- -99999999
old.logL <- -88888888
if (ite != 1)
t1 <- create.starting.values.nbinom(S, G, n, form, datsp)
t1 <- create.starting.values.nbinom.kmeans(S, G, n, form,
datsp)
pi <- t1$pi
fmM <- t1$fmM
tau <- t1$tau
dat.tau <- t1$dat.tau
first.fit <- t1$first.fit
while (abs(logL - old.logL) > 1e-04 & ite <= ite.max) {
old.logL <- logL
for (i in 1:G) {
pi[i] <- sum(tau[, i])/S
}
if (any(pi == 0)) {
cat("pi has gone to zero - restarting fitting \n")
t1 <- create.starting.values.nbinom(S, G, n, form,
datsp)
pi <- t1$pi
fmM <- t1$fmM
tau <- t1$tau
dat.tau <- t1$dat.tau
first.fit <- t1$first.fit
ite <- 1
}
fmM <- lapply(1:G, weighted.glm.nbinom, first.fit, tau,
n, fmM, sp)
for (j in 1:S) {
tmp <- rep(0, G)
for (g in 1:G) tmp[g] <- fmM[[g]]$sp.intercept[j]
tmp <- sum(tmp * tau[j, ])
for (g in 1:G) fmM[[g]]$sp.intercept[j] <- tmp
}
logL <- 0
tmp.like <- matrix(0, S, G)
est.tau <- lapply(1:S, estimate.pi.nbinom, sp, sp.name,
datsp, fmM, pi, G, first.fit)
for (j in 1:S) {
if (is.atomic(est.tau[[j]])) {
print(est.tau[[j]])
}
else {
tau[j, ] <- est.tau[[j]]$tau
logL <- logL + est.tau[[j]]$sum.like
}
}
if (trace)
cat(ite, " | ", logL, "\n")
ite <- ite + 1
}
fm.out <- data.frame(matrix(0, G, length(fmM[[1]]$coef)))
int.out <- rep(0, S)
names(fm.out) <- names(fmM[[1]]$coef)
tau <- data.frame(tau)
names(tau) <- paste("grp.", 1:G, sep = "")
EN <- -sum(unlist(tau) * log(unlist(tau)))
d <- length(unlist(fm.out)) + length(tau) - 1
fm.theta <- rep(0, G)
int.out <- fmM[[1]]$sp.intercept
for (i in 1:G) {
fm.out[i, ] <- fmM[[i]]$coef
}
names(pi) <- paste("G", 1:G, sep = ".")
t.pi <- additive.logistic(pi, TRUE)
parms <- c(t.pi[1:(G - 1)], unlist(fm.out), int.out, rep(1,
S))
logL.full <- logL
logL <- logLmix.nbinom(parms, first.fit, G, S, sp, sp.name)
options(warn = temp.warn)
if (full.model)
return(list(logl = logL, aic = -2 * logL + 2 * d, tau = round(tau,
4), pi = pi, bic = -2 * logL + log(S) * d, ICL = -2 *
logL + log(S) * d + 2 * EN, coef = fm.out, sp.intercept = int.out,
theta = fm.theta, fmM = fmM, model.tau = dat.tau,
covar = 0, aic.full = -2 * logL.full + 2 * d, bic.full = -2 *
logL.full + log(S) * d, pars = parms))
return(list(logl = logL, aic = -2 * logL + 2 * d, tau = round(tau,
4), pi = pi, bic = -2 * logL + log(S) * d, ICL = -2 *
logL + log(S) * d + 2 * EN, coef = fm.out, sp.intercept = int.out,
theta = fm.theta, covar = 0, aic.full = -2 * logL.full +
2 * d, bic.full = -2 * logL.full + log(S) * d, pars = parms))
}
"fitmix.nbinom.cpp" <-
function (form,datsp,sp,G=2,pars=NA,trace=TRUE,calc.hes=FALSE)
{
if(!is.numeric(sp)){
sp <- as.integer(factor(sp))
}
sp.name <- unique(sp)
S <- length(unique(sp))
n <- length(which(sp==sp.name[1]))
X <- model.matrix(form, data = datsp[sp==sp.name[1],])
## offset <- model.offset(form, data = datsp[sp==sp.name[1],])
offset <- model.frame(form, data = datsp[sp==sp.name[1],])
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,n)
##offset <- rep(0,n)
##X <- model.frame(form, data = datsp)
# X <- X[sp==sp.name[1],]
#y <- model.response(form, data = datsp)
y <- datsp$obs
if(is.na(pars[1])) {
##pars <- rep(0.01,G-1+(ncol(X)*G))
sp.int <- rep(0.5,S)
sp.dispersion <- rep(1,S)
fm <- matrix(runif(ncol(X)*G,-1,1),G,ncol(X))
pars <- c(runif(G-1),unlist(fm),sp.int,sp.dispersion)
}
## hes <- rep(0,length(pars)^2)
gradient <- rep(0,length(pars))
tau <- matrix(0,S,G) ##must leave this in as defines S & G
loglike <- try(.Call("SpeciesMix",pars,y,X,sp,tau,gradient,offset,as.integer(2),PACKAGE="SpeciesMix"))
calc.deriv <- function(p){
gradient <- rep(0,length(pars))
ll <- .Call("Calculate_Gradient",p,y,X,sp,tau,gradient,offset,as.integer(2),PACKAGE="SpeciesMix")
return(gradient)
}
r.deriv <- function(p){ logLmix.nbinom(p,list(y=y,x=model.matrix(form, data = datsp)),G,S,sp,sp.name,out.tau=FALSE)}
#r.grad <- nd2(pars,r.deriv)
##print(r.grad)
hes <- 0
covar <- 0
if(calc.hes){
hes <- nd2(pars,calc.deriv)
dim(hes) <- rep(length(pars),2)
dim(hes) <- rep(length(pars),2)
covar <- try(solve(hes))
#rownames(hes) <- colnames(hes) <- c(paste("G.",1:(G-1),sep=""),paste("G",1:G,rep(colnames(X),each=G),sep="."))
}
if(!is.numeric(loglike)) loglike <- 0
pi <- pars[1:(G-1)]
#coef <- pars[ (G):length(pars)]
coef <- pars[-1*(1:((G-1)))] ## remove pi
sp.int <- coef[(length(coef)-(2*S-1)):(length(coef)-S)]
sp.dispersion <- coef[(length(coef)-(S-1)):length(coef)]
fm <- coef[-1*((length(coef)-(2*S-1)):length(coef))]
offset <- model.frame(form, data = datsp)
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,nrow(datsp))
r.logl <- logLmix.nbinom(pars,list(y=y,x=model.matrix(form, data = datsp),offset=offset),G,S,sp,sp.name,out.tau=TRUE)
print(r.logl$logl)
pi <- additive.logistic(pi)
names(pi) <- paste("G.",1:G,sep="")
coef <- matrix(fm,G,ncol(X))
rownames(coef) <- paste("G.",1:G,sep="")
colnames(coef) <- colnames(X)
AIC <- 2*loglike + 2*length(pars)
BIC <- 2*loglike + log(S)*length(pars)
##list(logl=loglike,r.logl=r.logl$logl,pi=pi,coef=coef,tau=round(exp(r.logl$tau),4),aic=AIC,bic=BIC,hessian=hes,gradient=gradient)
list(logl=loglike,pi=pi,coef=coef,sp.intercept=sp.int,sp.dispersion=sp.dispersion,tau=round(exp(r.logl$tau),4),aic=AIC,bic=BIC,hessian=hes,gradient=gradient,covar=covar)#,r.grad=r.grad)
}
"fitmix.tweedie" <-
function (form, datsp, sp, G = 2, ite.max = 500, trace = TRUE,
full.model = FALSE)
{
temp.warn <- getOption("warn")
options(warn = -1)
sp.name <- unique(sp)
S <- length(unique(sp))
n <- length(which(sp == sp.name[1]))
cat("Fitting Group", G, "\n")
if (trace)
cat("Iteration | LogL \n")
dat.tau <- 0
pi <- rep(0, G)
ite <- 1
logL <- -99999999
old.logL <- -88888888
if (ite != 1)
t1 <- create.starting.values.tweedie(S, G, n, form, datsp)
t1 <- create.starting.values.tweedie.kmeans(S, G, n, form,
datsp)
pi <- t1$pi
fmM <- t1$fmM
tau <- t1$tau
dat.tau <- t1$dat.tau
first.fit <- t1$first.fit
while (abs(logL - old.logL) > 1e-04 & ite <= ite.max) {
old.logL <- logL
for (i in 1:G) {
pi[i] <- sum(tau[, i])/S
}
if (any(pi == 0)) {
cat("pi has gone to zero - restarting fitting \n")
t1 <- create.starting.values.tweedie(S, G, n, form,
datsp)
pi <- t1$pi
fmM <- t1$fmM
tau <- t1$tau
dat.tau <- t1$dat.tau
first.fit <- t1$first.fit
ite <- 1
}
fmM <- lapply(1:G, weighted.glm.tweedie, first.fit, tau,
n, fmM, sp)
for (j in 1:S) {
tmp <- rep(0, G)
for (g in 1:G) tmp[g] <- fmM[[g]]$sp.intercept[j]
tmp <- sum(tmp * tau[j, ])
for (g in 1:G) fmM[[g]]$sp.intercept[j] <- tmp
}
logL <- 0
tmp.like <- matrix(0, S, G)
est.tau <- lapply(1:S, estimate.pi.tweedie, sp, sp.name,
datsp, fmM, pi, G, first.fit)
for (j in 1:S) {
if (is.atomic(est.tau[[j]])) {
print(est.tau[[j]])
}
else {
tau[j, ] <- est.tau[[j]]$tau
logL <- logL + est.tau[[j]]$sum.like
}
}
if (trace)
cat(ite, " | ", logL, "\n")
ite <- ite + 1
}
fm.out <- data.frame(matrix(0, G, length(fmM[[1]]$coef)))
int.out <- fmM[[1]]$sp.intercept
names(fm.out) <- names(fmM[[1]]$coef)
tau <- data.frame(tau)
names(tau) <- paste("grp.", 1:G, sep = "")
EN <- -sum(unlist(tau) * log(unlist(tau)))
d <- length(unlist(fm.out)) + length(tau) - 1
fm.phi <- fm.p <- rep(0, G)
for (i in 1:G) {
fm.out[i, ] <- fmM[[i]]$coef
fm.phi[i] <- fmM[[i]]$phi
fm.p[i] <- fmM[[i]]$p
}
names(pi) <- paste("G", 1:G, sep = ".")
t.pi <- additive.logistic(pi, TRUE)
parms <- c(t.pi[1:(G - 1)], unlist(fm.out), int.out, rep(2,
S))
names(parms) <- c(rep("pi", G - 1), rep("coef", length(unlist(fm.out))),
rep("int", length(int.out)), rep("phi", S))
logL.full <- logL
logL <- logLmix.tweedie(parms, first.fit, G, S, sp, sp.name)
print(logL)
options(warn = temp.warn)
if (full.model)
return(list(logl = logL, aic = -2 * logL + 2 * d, tau = round(tau,
4), pi = pi, bic = -2 * logL + log(S) * d, ICL = -2 *
logL + log(S) * d + 2 * EN, coef = fm.out, sp.intercept = int.out,
phi = fm.phi, p = fm.p, fmM = fmM, model.tau = dat.tau,
covar = 0, aic.full = -2 * logL.full + 2 * d, bic.full = -2 *
logL.full + log(S) * d, pars = parms))
return(list(logl = logL, aic = -2 * logL + 2 * d, tau = round(tau,
4), pi = pi, bic = -2 * logL + log(S) * d, ICL = -2 *
logL + log(S) * d + 2 * EN, coef = fm.out, sp.intercept = int.out,
phi = fm.phi, p = fm.p, covar = 0, aic.full = -2 * logL.full +
2 * d, bic.full = -2 * logL.full + log(S) * d, pars = parms))
}
"fitmix.tweedie.cpp" <-
function (form, datsp, sp, G = 2, pars = NA, trace = TRUE, calc.hes = FALSE)
{
if (!is.numeric(sp)) {
sp <- as.integer(factor(sp))
}
sp.name <- unique(sp)
S <- length(unique(sp))
n <- length(which(sp == sp.name[1]))
X <- model.matrix(form, data = datsp[sp == sp.name[1], ])
offset <- model.frame(form, data = datsp[sp==sp.name[1],])
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,n)
y <- datsp$obs
if (is.na(pars[1])) {
sp.int <- rep(0.5, S)
sp.phi <- rep(1, S)
sp.p <- rep(1.6, S)
fm <- matrix(rep(0.5, ncol(X) * G), G, ncol(X))
pars <- c(rep(0.5, G - 1), unlist(fm), sp.int, sp.phi)
}
pars.og <- pars
pars.og[1] <- pars.og[1] * 0.99
gradient <- rep(0, length(pars))
tau <- matrix(0, S, G)
loglike <- try(.Call("SpeciesMix", pars, y, X, sp, tau, gradient,
offset, as.integer(3),PACKAGE="SpeciesMix"))
calc.deriv <- function(p) {
gradient <- rep(0, length(pars))
ll <- .Call("Calculate_Gradient", p, y, X, sp, tau, gradient,
offset, as.integer(3),PACKAGE="SpeciesMix")
return(gradient)
}
hes <- 0
if (calc.hes) {
hes <- jacobian(calc.deriv, pars, method = "simple")
dim(hes) <- rep(length(pars), 2)
dim(hes) <- rep(length(pars), 2)
}
if (!is.numeric(loglike))
loglike <- 0
pi <- pars[1:(G - 1)]
coef <- pars[-1 * (1:((G - 1)))]
sp.int <- coef[(length(coef) - (2 * S - 1)):(length(coef) -
S)]
sp.dispersion <- coef[(length(coef) - (S - 1)):length(coef)]
sp.phi <- sp.dispersion[1:S]
sp.p <- rep(1.6, S)
fm <- coef[-1 * ((length(coef) - (2 * S - 1)):length(coef))]
offset <- model.frame(form, data = datsp)
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,nrow(datsp))
names(gradient) <- names(pars) <- c(rep("pi", G - 1), rep("coef",
length(unlist(fm))), rep("int", length(sp.int)), rep("phi",
S))
r.deriv <- function(p) {
logLmix.tweedie(p, list(y = y, x = model.matrix(form,
data = datsp)), G, S, sp, sp.name, out.tau = FALSE)
}
r.logl <- logLmix.tweedie(pars, list(y = y, x = model.matrix(form,
data = datsp),offset=offset), G, S, sp, sp.name, out.tau = TRUE)
print(r.logl$logl)
pi <- additive.logistic(pi)
names(pi) <- paste("G.", 1:G, sep = "")
coef <- matrix(fm, G, ncol(X))
rownames(coef) <- paste("G.", 1:G, sep = "")
colnames(coef) <- colnames(X)
AIC <- 2 * loglike + 2 * length(pars)
BIC <- 2 * loglike + log(S) * length(pars)
list(logl = loglike, pi = pi, coef = coef, sp.intercept = sp.int,
phi = sp.phi, p = sp.p, tau = round(exp(r.logl$tau),
4), aic = AIC, bic = BIC, hessian = hes, gradient = gradient)
}
"glm.fit.nbinom" <-
function (x,y,offset=NULL,weights=NULL,mustart=NULL,est.var=FALSE)
{
X <- x
if(is.null(offset)) offset <- 0
if(is.null(weights)) weights <- rep(1,length(y))
gradient <- rep(0,ncol(X)+1)
if(is.null(mustart)){ pars <- gradient+1} else{pars <- mustart}
fitted.values <- rep(0,length(y))
logl <- .Call("Neg_Bin",pars,X,y,weights,offset,gradient,fitted.values,PACKAGE="SpeciesMix")
vcov <- 0
se <- rep(0,length(pars))
if(est.var) {
calc.deriv <- function(p){
gradient <- rep(0,length(pars))
ll <- .Call("Neg_Bin_Gradient",p,X,y,weights,offset,gradient,PACKAGE="SpeciesMix")
return(gradient)
}
hes <- nd2(pars,calc.deriv)
dim(hes) <- rep(length(pars),2)
vcov <- try(solve(hes))
se <- try(sqrt(diag(vcov)))
colnames(vcov) <- rownames(vcov) <- c("theta",colnames(X))
}
names(pars) <- names(se) <- names(gradient) <- c("theta",colnames(X))
return(list(logl=logl,coef=pars[-1],theta=pars[1],se=se[-1],se.theta=se[1],fitted=fitted.values,gradient=gradient,vcov=vcov))
}
"glm.nbinom" <-
function (form,data,weights=NULL,mustart=NULL,est.var=FALSE)
{
X <- model.matrix(form,data)
t1 <- model.frame(form,data)
y <- model.response(t1)
offset <- model.offset(t1)
if(is.null(offset)) offset <- 0
if(is.null(weights)) weights <- rep(1,length(y))
gradient <- rep(0,ncol(X)+1)
if(is.null(mustart)){ pars <- gradient+1}else{pars <- mustart}
fitted.values <- rep(0,length(y))
logl <- .Call("Neg_Bin",pars,X,y,weights,offset,gradient,fitted.values,PACKAGE="SpeciesMix")
vcov <- 0
se <- rep(0,length(pars))
if(est.var) {
calc.deriv <- function(p){
gradient <- rep(0,length(pars))
ll <- .Call("Neg_Bin_Gradient",p,X,y,weights,offset,gradient,PACKAGE="SpeciesMix")
return(gradient)
}
hes <- nd2(pars,calc.deriv)
dim(hes) <- rep(length(pars),2)
vcov <- try(solve(hes))
se <- try(sqrt(diag(vcov)))
colnames(vcov) <- rownames(vcov) <- c("theta",colnames(X))
}
names(pars) <- names(se) <- names(gradient) <- c("theta",colnames(X))
return(list(logl=logl,coef=pars[-1],theta=pars[1],se=se[-1],se.theta=se[1],fitted=fitted.values,gradient=gradient,vcov=vcov))
}
"ldTweedie.lp" <-
function ( parms, y, X.p, offsetty, phi, p, wts=rep( 1, length( y)))
{
mu <- exp( X.p %*% parms[1:ncol( X.p)] + offsetty)
if( is.null( phi) & is.null( p)){
phi <- parms[ncol( X.p) + 1]
p <- parms[ncol( X.p) + 2]
}
if( is.null( phi) & !is.null( p))
phi <- parms[ncol( X.p)+1]
if( !is.null( phi) & is.null( p))
p <- parms[ncol( X.p)+1]
lambda <- ( mu^( 2-p)) / ( phi*(2-p))
alpha <- ( 2-p) / ( p-1)
tau <- phi*(p-1)*mu^(p-1)
mu.Z <- alpha * tau
return( -sum( wts * dPoisGam( y, lambda=lambda, mu.Z=mu.Z, alpha=alpha, LOG=TRUE)))
}
"ldTweedie.lp.deriv" <-
function ( parms, y, X.p, offsetty, phi, p, wts=rep( 1, length( y)))
{
mu <- exp( X.p %*% parms[1:ncol( X.p)] + offsetty)
p.flag <- phi.flag <- FALSE
if( is.null( phi) & is.null( p)){
p.flag <- phi.flag <- TRUE
phi <- parms[ncol( X.p) + 1]
p <- parms[ncol( X.p) + 2]
}
if( is.null( phi) & !is.null( p)){
phi <- parms[ncol( X.p)+1]
phi.flag <- TRUE
}
if( !is.null( phi) & is.null( p)){
p <- parms[ncol( X.p)+1]
p.flag <- TRUE
}
lambda <- ( mu^( 2-p)) / ( phi*(2-p))
alpha <- ( 2-p) / ( p-1)
tau <- phi*(p-1)*mu^(p-1)
mu.Z <- alpha * tau
dTweedparms <- -wts * dPoisGamDerivs( y, lambda=lambda, mu.Z=mu.Z, alpha=alpha)
DTweedparmsDmu <- matrix( c( ( mu^(1-p)) / phi, alpha*phi*( ( p-1)^2)*( mu^(p-2)), rep( 0, length( mu))), nrow=3, byrow=T)
tmp <- rowSums( dTweedparms * t( DTweedparmsDmu))
tmp <- tmp * mu
tmp <- apply( X.p, 2, function( x) x*tmp)
derivs <- colSums( tmp)
if( phi.flag){
DTweedparmsDphi <- matrix( c( -( ( mu^(2-p)) / ( ( phi^2)*(2-p))), alpha*( p-1)*( mu^( p-1)), rep( 0, length( mu))), nrow=3, byrow=T)
tmpPhi <- rowSums( dTweedparms * t( DTweedparmsDphi))#vectorised way of doing odd calculation
derivs <- c( derivs, sum( tmpPhi))
names( derivs)[length( derivs)] <- "phi"
}
if( p.flag){
dalphadp <- -( 1+alpha) / ( p-1)
DTweedparmsDp <- matrix( c( lambda*( 1/(2-p) - log( mu)), mu.Z*( dalphadp/alpha + 1/( p-1) + log( mu)), rep( dalphadp, length( y))), nrow=3, byrow=T)
tmpP <- rowSums( dTweedparms * t( DTweedparmsDp))
derivs <- c( derivs, sum( tmpP))
names( derivs)[length( derivs)] <- "p"
}
return( derivs)
}
"Lmix" <-
function (pars,y,x,G)
{
fm <- pars[-1*(1:(G-1))]
pi <- pars[(1:(G-1))]
dim(fm) <- c(G,(length(pars)-(G-1))/G)
pi <- additive.logistic(pi)
S <- dim(y)[2]
link.fun <- make.link("logit")
like <- 1
lf <- matrix(0,G,S)
lg <- rep(0,S)
dBi <- array(0,dim=c(S,G,dim(fm)[2]))
for(s in 1:S){
tmp.like <- rep(0,G)
for(g in 1:G){
lpre <- x%*%fm[g,]
for(j in 1:dim(fm)[2]){
dBi[s,g,j] <- sum((y[,s]-link.fun$linkinv(lpre))*x[,j])
}
tmp.like[g] <- pi[g]*prod(dbinom(y[,s],1,link.fun$linkinv(lpre),log=F))
lf[g,s] <- prod(dbinom(y[,s],1,link.fun$linkinv(lpre),log=F))
}
lg[s] <- sum(tmp.like)
like <- like*sum(tmp.like)
}
#print(dBi)
# print(log(lg))
# print(log(lf))
dl.dpi <- rep(0,G)
for(g in 1:G) dl.dpi[g] <- ( sum( exp(-log(lg) + log(lf[g,]))))
der <- matrix(0,dim(fm)[1],dim(fm)[2])
for(j in 1:dim(fm)[2]){
for(g in 1:G){
for(s in 1:S){
der[g,j] <- der[g,j]+ exp( -log(lg[s]) + log(pi[g]) + log(lf[g,s])) * dBi[s,g,j]
#if(g==1 & j == 1) print(c(der[g,j],-log(lg[s]),pi[g],log(lf[g,s]),dBi[s,g,j]))
}
}
}
dpi.deta <- matrix(0,G-1,G)
ad.trans <- 1+sum(exp(pars[1:(G-1)]))
for(i in 1:(G-1))
for(g in 1:G-1){
if(i==g) {dpi.deta[i,g] <- exp(pars[i])*exp(pars[g])/ad.trans^2}
else{ dpi.deta[i,g] <- exp(pars[i])/ad.trans - exp(2*pars[g])/ad.trans^2}
}
dpi.deta[,G] <- -rowSums(dpi.deta)
print(dpi.deta)
print(dl.dpi)
d1 <- dpi.deta%*%dl.dpi
print(d1)
list(like,0-der)
}
"logLike.pars" <-
function (pi,coef,sp.form,sp.data,covar.data)
{
G <- length(pi)
pars <- c(additive.logistic(pi,T)[1:(G-1)],unlist(coef))
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=unlist(sp.data),covar.data)
names(data)[1] <- as.character(sp.form)[2]
first.fit <- list(y=data[,1],x=model.matrix(sp.form,data=data))
logl <- -logLmix(pars,first.fit,G,S,sp,sp.name)
logl
}
"logLmix" <-
function (pars,first.fit,G,S,sp,spname,out.tau=FALSE)
{
tau <- matrix(0,S,G)
##tau,out.tau=FALSE
if(G>1){
fm <- pars[-1*(1:(G-1))]
pi <- pars[(1:(G-1))]
## fm <- tau[-1*((length(tau)-(G-2)):length(tau))]
#pi <- tau[((length(tau)-(G-2)):length(tau))]
dim(fm) <- c(G,(length(pars)-(G-1))/G)
##pi[G] <- 1-sum(pi)
pi <- additive.logistic(pi)
} else{
fm <- tau[1:(length(pars)-1)]
dim(fm) <- c(1,length(fm))
pi <- 1
}
link.fun <- make.link("logit")
log.like <- 0
for(j in 1:S){
sel.sp <- which(sp==spname[j])
tmp.like <- rep(0,G)
for(i in 1:G){
if(length(fm[i,])==1){lpre <- first.fit$x[sel.sp,]*fm[i,]
}else{ lpre <- first.fit$x[sel.sp,]%*%fm[i,]}
tmp.like[i] <- sum(dbinom(first.fit$y[sel.sp],1,link.fun$linkinv(lpre),log=T))
}
eps <- max(tmp.like)
log.like <- log.like + (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
tau[j,] <- log(pi) + tmp.like - (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
}
if(out.tau)return(list(logl=log.like,tau=tau))
log.like
}
"logLmix.gaussian" <-
function (pars,first.fit,G,S,sp,spname,out.tau=FALSE)
{
tau <- matrix(0,S,G)
##tau,out.tau=FALSE
if(G>1){
fm <- pars[-1*(1:((G-1)+G))] ## remove pi
sp.int <- fm[(length(fm)-(S*G-1)):length(fm)]
fm <- fm[-1*((length(fm)-(S*G-1)):length(fm))]
pi <- pars[(1:(G-1))]
theta <- pars[G:(G+G-1)]
## fm <- tau[-1*((length(tau)-(G-2)):length(tau))]
#pi <- tau[((length(tau)-(G-2)):length(tau))]
dim(sp.int) <- c(S,G)
dim(fm) <- c(G,length(fm)/G)
##pi[G] <- 1-sum(pi)
pi <- additive.logistic(pi)
} else{
return(0)
fm <- tau[1:(length(pars)-1)]
dim(fm) <- c(1,length(fm))
pi <- 1
}
log.like <- 0
for(j in 1:S){
sel.sp <- which(sp==spname[j])
tmp.like <- rep(0,G)
for(i in 1:G){
lpre <- cbind(1,first.fit$x[sel.sp,])%*%c(sp.int[j,i],fm[i,])
##tmp.like[i] <- sum(dpois(first.fit$y[sel.sp],exp(lpre),log=T))
tmp.like[i] <- sum(dnorm(first.fit$y[sel.sp],mean=lpre,sd=theta[i],log=T))
eps <- max(tmp.like)
log.like <- log.like + (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
tau[j,] <- log(pi) + tmp.like - (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
}
}
if(out.tau)return(list(logl=log.like,tau=tau))
log.like
}
"logLmix.nbinom" <-
function (pars,first.fit,G,S,sp,spname,out.tau=FALSE)
{
tau <- matrix(0,S,G)
##tau,out.tau=FALSE
if(G>1){
fm <- pars[-1*(1:((G-1)))] ## remove pi
##sp.int <- fm[(length(fm)-(S*G-1)):length(fm)]
sp.int <- fm[(length(fm)-(2*S-1)):(length(fm)-S)]
sp.dispersion <- fm[(length(fm)-(S-1)):length(fm)]
##fm <- fm[-1*((length(fm)-(S*G-1)):length(fm))]
fm <- fm[-1*((length(fm)-(2*S-1)):length(fm))]
pi <- pars[(1:(G-1))]
theta <- pars[G:(G+G-1)]
## fm <- tau[-1*((length(tau)-(G-2)):length(tau))]
#pi <- tau[((length(tau)-(G-2)):length(tau))]
##dim(sp.int) <- c(S,G)
dim(fm) <- c(G,length(fm)/G)
##pi[G] <- 1-sum(pi)
pi <- additive.logistic(pi)
} else{
return(0)
fm <- tau[1:(length(pars)-1)]
dim(fm) <- c(1,length(fm))
pi <- 1
}
log.like <- 0
for(j in 1:S){
sel.sp <- which(sp==spname[j])
tmp.like <- rep(0,G)
for(i in 1:G){
lpre <- cbind(1,first.fit$x[sel.sp,])%*%c(sp.int[j],fm[i,])+first.fit$offset[sel.sp]##,i],fm[i,])
##tmp.like[i] <- sum(dpois(first.fit$y[sel.sp],exp(lpre),log=T))
tmp.like[i] <- sum(dnbinom(first.fit$y[sel.sp],mu=exp(lpre),size=sp.dispersion[j],log=T))
}
# print(tmp.like)
eps <- max(tmp.like)
log.like <- log.like + (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
tau[j,] <- log(pi) + tmp.like - (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
}
if(out.tau)return(list(logl=log.like,tau=tau))
log.like
}
"logLmix.tweedie" <-
function (pars,first.fit,G,S,sp,spname,out.tau=FALSE)
{
tau <- matrix(0,S,G)
##tau,out.tau=FALSE
if(G>1){
pi <- pars[which(names(pars)=="pi")]
phi <- pars[which(names(pars)=="phi")]
p <- rep(1.6,S)
fm <- pars[which(names(pars)=="coef")]
sp.int <- pars[which(names(pars)=="int")]
##fm <- pars[-1*(1:((G-1)+G))] ## remove pi
##sp.int <- fm[(length(fm)-(S*G-1)):length(fm)]
##fm <- fm[-1*((length(fm)-(S*G-1)):length(fm))]
##pi <- pars[(1:(G-1))]
##theta <- pars[G:(G+G-1)]
## fm <- tau[-1*((length(tau)-(G-2)):length(tau))]
#pi <- tau[((length(tau)-(G-2)):length(tau))]
##dim(sp.int) <- c(S,G)
dim(fm) <- c(G,length(fm)/G)
##pi[G] <- 1-sum(pi)
pi <- additive.logistic(pi)
} else{
return(0)
fm <- tau[1:(length(pars)-1)]
dim(fm) <- c(1,length(fm))
pi <- 1
}
log.like <- 0
for(j in 1:S){
sel.sp <- which(sp==spname[j])
tmp.like <- rep(0,G)
for(i in 1:G){
lpre <- cbind(1,first.fit$x[sel.sp,])%*%c(sp.int[j],fm[i,])+first.fit$offset[sel.sp]
##tmp.like[i] <- sum(dpois(first.fit$y[sel.sp],exp(lpre),log=T))
tmp.like[i] <- sum(dTweedie(y=first.fit$y[sel.sp],mu=exp(lpre),phi=phi[j],p=p[j],LOG=T))
}
##print(tmp.like)
eps <- max(tmp.like)
log.like <- log.like + (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
tau[j,] <- log(pi) + tmp.like - (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
}
if(out.tau)return(list(logl=log.like,tau=tau))
log.like
}
"mix.residuals" <-
function (fmM,form,datsp,sp)
{
cat("calculating residuals \n")
link.fun <- make.link("logit")
x <- model.matrix(form,data=datsp)
spname <- unique(sp)
S <- length(spname)
G <- ncol(fmM$tau)
PIT <- matrix(NA,S,G)
for(g in 1:G){
for(s in 1:length(spname)){
sel.sp <- which(sp==spname[s])
t.obs <- sum(datsp$obs[sel.sp])
pre <- link.fun$linkinv(x[sel.sp,]%*%fmM$coef[g,])
obs <- datsp$obs[sel.sp]
dis <- distr.binom(pre)
# PIT[s,g] <- qnorm(sum(dis[1:(length(dis)-1)],dis[length(dis)]/2))
nSucc <- sum( obs)
transfo <- sum( dis[1:nSucc],dis[nSucc+1]/2)
transfo <- min( transfo, 1)
##transfor <- max( transfo, 0)
PIT[s,g] <- qnorm( transfo)
}
}
PIT
}
"nd2" <-
function( x0, f, m=NULL, D.accur=4, ...) {
# A function to compute highly accurate first-order derivatives
# From Fornberg and Sloan (Acta Numerica, 1994, p. 203-267; Table 1, page 213)
# Adapted by Scott Foster from code nicked off the net 2007
D.n <- length( x0)
if ( is.null( m)) {
D.f0 <- f(x0, ...)
m <- length( D.f0)
}
if ( D.accur == 2) {
D.w <- tcrossprod( rep( 1, m),c( -1/2, 1/2))
D.co <- c( -1, 1)
}
else {
D.w <- tcrossprod( rep( 1, m),c( 1/12, -2/3, 2/3, -1/12))
D.co <- c( -2, -1, 1, 2)
}
D.n.c <- length( D.co)
macheps <- .Machine$double.eps
D.h <- macheps^( 1/3)*abs( x0)
D.deriv <- matrix( NA, nrow=m, ncol=D.n)
for ( ii in 1:D.n) {
D.temp.f <- matrix( 0, m, D.n.c)
for ( jj in 1:D.n.c) {
D.xd <- x0+D.h[ii]*D.co[jj]*( 1:D.n == ii)
D.temp.f[,jj] <- f( D.xd, ...)
}
D.deriv[,ii] <- rowSums( D.w*D.temp.f)/D.h[ii]
}
return( as.double( D.deriv))
}
"nH2" <-
function( pt, fun, accur=c(4,4), type="H.Diag", ...) {
# A function to compute highly accurate second order Hessian diags and other off-diags
# partially from Fornberg and Sloan (Acta Numerica, 1994, p. 203-267; Table 1, page 213)
# Adapted by Scott Foster from code nicked off the net 2007
H.n <- length( pt)
derivs <- function( d.x0, ...) { nd2( x0=d.x0, f=fun, m=1, D.accur=accur[2], ...) }
Hes <- nd2( x0=pt, f=derivs, D.accur=accur[2], ...)
Hes <- matrix( Hes, nrow=length( pt))
Hes <- ( Hes+t( Hes))/2
if ( type == "H.Diag") {
macheps <- .Machine$double.eps
H.h <- macheps^(1/4)*abs(pt)
H.f0 <- fun( pt, ...)
H.m <- length( H.f0)
if ( accur[1] == 2) {
H.w <- tcrossprod( rep( 1, H.m), c( 1, -2, 1))
H.co <- c( -1, 0, 1)
}
else {
H.w <- tcrossprod( rep( 1, H.m), c( -1/12, 4/3, -5/2, 4/3, -1/12))
H.co <- c( -2, -1, 0, 1, 2)
}
H.n.c <- length( H.co)
Hes.diag <- double( length=H.n)
for ( ii in 1:H.n) {
H.temp.f <- matrix( 0, H.m, H.n.c)
for ( jj in 1:H.n.c) {
if ( H.co[jj] != 0) {
H.xd <- pt+H.h[ii]*H.co[jj]*( 1:H.n == ii)
H.temp.f[,jj] <- fun( H.xd, ...)
}
else
H.temp.f[,jj] <- H.f0
}
Hes.diag[ii] <- rowSums( H.w*H.temp.f)/( H.h[ii]^2)
}
diag( Hes) <- Hes.diag
}
return( Hes)
}
".onLoad" <-
function (libname, pkgname)
{
##require(MASS)
# Generic DLL loader
dll.path <- file.path( libname, pkgname, 'libs')
if( nzchar( subarch <- .Platform$r_arch))
dll.path <- file.path( dll.path, subarch)
this.ext <- paste( sub( '.', '[.]', .Platform$dynlib.ext,
fixed=TRUE), '$', sep='')
dlls <- dir( dll.path, pattern=this.ext, full.names=FALSE)
names( dlls) <- dlls
if( length( dlls))
lapply( dlls, function( x) library.dynam( sub( this.ext, '', x),
package=pkgname, lib.loc=libname))
}
"predict.archetype" <-
function (object,new.obs,...)
{
mixture.model <- object
if(class(mixture.model)[2]=="bernoulli"){
G <- length(mixture.model$pi)
covar <- mixture.model$covar[-(1:(G-1)),-(1:(G-1))]
coef <- mixture.model$coef
model.fm <- as.formula(mixture.model$formula)
model.fm[[2]] <- NULL
X <- model.matrix(model.fm,new.obs)
link.fun <- make.link("logit")
outvar <- matrix(NA,dim(X)[1],G)
outpred <- matrix(NA,dim(X)[1],G)
colnames(outvar) <- colnames(outpred) <- paste("G",1:G,sep=".")
for(g in 1:G){
lp <- as.numeric(X%*%coef[g,])
outpred[,g] <- link.fun$linkinv(lp)
dhdB <- (exp(lp)/(1+exp(lp)))*X - exp(lp)^2/((1+exp(lp))^2)*X
c2 <- covar[seq(g,dim(covar)[1],G),seq(g,dim(covar)[1],G)]
for(k in 1:dim(X)[1]){
outvar[k,g] <- (dhdB[k,]%*%c2)%*%(dhdB[k,])
}
}
}
if(class(mixture.model)[2]=="negbin"|class(mixture.model)[2]=="tweedie"){
G <- length(mixture.model$pi)
##covar <- mixture.model$covar[-1*c(1:(G-1),(dim(mixture.model$covar)[1]-2*length(mixture.model$sp.intercept)+1):dim(mixture.model$covar)[1]),-1*c(1:(G-1),(dim(mixture.model$covar)[1]-2*length(mixture.model$sp.intercept)+1):dim(mixture.model$covar)[1])] # remove pi,sp.int,sp.disp from matrix
covar <- mixture.model$covar[-1*c(1:(G-1),(dim(mixture.model$covar)[1]-length(mixture.model$sp.intercept)+1):dim(mixture.model$covar)[1]),-1*c(1:(G-1),(dim(mixture.model$covar)[1]-length(mixture.model$sp.intercept)+1):dim(mixture.model$covar)[1])] # remove pi,sp.disp from matrix
sp.int <- mixture.model$sp.intercept
coef <- mixture.model$coef
model.fm <- as.formula(mixture.model$formula)
model.fm[[2]] <- NULL
X <- cbind(model.matrix(model.fm,new.obs),1)
offset <- model.frame(model.fm, data = new.obs)
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,nrow(X))
outvar <- matrix(NA,dim(X)[1],G)
outpred <- matrix(NA,dim(X)[1],G)
colnames(outvar) <- colnames(outpred) <- paste("G",1:G,sep=".")
for(g in 1:G){
s.outvar <- matrix(NA,dim(X)[1],length(sp.int))
s.outpred <- matrix(NA,dim(X)[1],length(sp.int))
for(s in 1:length(sp.int)){
lp <- as.numeric(X%*%c(coef[g,],sp.int[s])+offset)
s.outpred[,s] <- exp(lp)
dhdB <- exp(lp)*X
c2 <- covar[c(seq(g,G*(dim(X)[2]-1),G),G*(dim(X)[2]-1)+s),c(seq(g,G*(dim(X)[2]-1),G),G*(dim(X)[2]-1)+s)]
for(k in 1:dim(X)[1]){
s.outvar[k,s] <- (dhdB[k,]%*%c2)%*%(dhdB[k,])
}
}
outpred[,g] <- apply(s.outpred*rep(mixture.model$tau[,g],each=dim(X)[1]),1,mean)/sum(mixture.model$tau[,g])
outvar[,g]<- apply(s.outvar*rep(mixture.model$tau[,g],each=dim(X)[1]),1,mean)/sum(mixture.model$tau[,g])
}
}
list(fit=outpred,se.fit=sqrt(outvar))
}
"print.archetype" <-
function (x,...)
{
cat("\nMixing probabilities\n")
print(x$pi)
cat("\nCoefficents\n")
print(x$coef)
if(!is.na(x$se[1])){
cat("\nStandard Errors of coefficents\n")
print(x$se)
}
cat("\nPosterior Probabilities\n")
print(x$tau)
}
"pTweedie" <-
function ( quant, mu, phi, p)
{
lambda <- ( mu^( 2-p)) / ( phi*(2-p))
alpha <- ( 2-p) / ( p-1)
tau <- phi*(p-1)*mu^(p-1)
mu.Z <- alpha * tau
ps <- 0
# ps <- pPoisGam( quant, lambda, mu.Z, alpha)
# require( tweedie)
# ps <- ptweedie( as.numeric( q), mu=as.numeric( mu), phi=as.numeric( phi), power=as.numeric( p))
return( ps)
}
"rPoisGam" <-
function ( n, lambda, mu.Z, alpha)
{
mu.N <- lambda
#simulate n random variables from the same compound poisson distribution
my.fun <- function (parms)
return( rgamma( n=1, scale=parms[3], shape=parms[1]*parms[2]))
# return( sum( rgamma( n=parms[1], scale=parms[3], shape=parms[2])))
tmp <- c( n, length( mu.N), length(mu.Z), length( alpha))
names( tmp) <- c( "n","mu.N","mu.Z","alpha")
if( !all( is.element( tmp[-1], c( 1, tmp[1])))) {
print( "rPoisGam: error -- length of arguments are not compatible")
return( tmp)
}
if( tmp["mu.N"]==1)
mu.N <- rep( mu.N, tmp["n"])
if( tmp["mu.Z"]==1)
mu.Z <- rep( mu.Z, tmp["n"])
if( tmp["alpha"]==1)
alpha <- rep( alpha, tmp["n"])
np <- matrix( rpois( n, mu.N), ncol=1)
beta <- mu.Z / alpha
y <- apply( cbind( np, alpha, beta), 1, my.fun)
return( y)
}
"rTweedie" <-
function ( n, mu, phi, p)
{
lambda <- ( mu^( 2-p)) / ( phi*(2-p))
alpha <- ( 2-p) / ( p-1)
tau <- phi*(p-1)*mu^(p-1)
mu.Z <- alpha * tau
rans <- rPoisGam( n, lambda, mu.Z, alpha)
return( rans)
}
"SpeciesMix" <-
function (sp.form,sp.data,covar.data,G=2, pars=NA, em.prefit=TRUE,em.steps=3, em.refit = 1 , dist="bernoulli",est.var = FALSE,residuals=FALSE,trace=TRUE)
{
if(dist=="bernoulli") return(SpeciesMix.bernoulli(sp.form,sp.data,covar.data,G, pars, em.prefit,em.steps, em.refit ,est.var,residuals,trace,r1=FALSE))
## if(dist=="bernoulli-random") return(SpeciesMix.bernoulli(sp.form,sp.data,covar.data,G, pars, em.prefit,em.steps, em.refit ,est.var,residuals,trace,r1=TRUE))
if(dist=="negbin") return(SpeciesMix.nbinom(sp.form,sp.data,covar.data,G, pars, em.prefit,em.steps, em.refit ,est.var,residuals,trace))
if(dist=="tweedie") return(SpeciesMix.tweedie(sp.form,sp.data,covar.data,G, pars, em.prefit,em.steps, em.refit ,est.var,residuals,trace))
##if(dist=="gaussian") return(SpeciesMix.gaussian(sp.form,sp.data,covar.data,G, pars, em.prefit,em.steps, em.refit ,est.var,residuals,trace))
print("incorrect distribution type, options are bernoulli, negbin or tweedie")
return(0)
}
"SpeciesMix.bernoulli" <-
function (sp.form,sp.data,covar.data,G=2, pars=NA, em.prefit=TRUE,em.steps=4, em.refit = 1 , est.var = FALSE,residuals=FALSE,trace=TRUE,r1=FALSE)
{
t.covar.data <- covar.data
t.sp.data <- sp.data
sp.form <- update.formula(sp.form,obs~1+.)
if(em.prefit | G==1){
prefit <- SpeciesMix.em(sp.form,sp.data,covar.data,G,ite.max=em.steps,em.refit=em.refit,r1)
pars <- c(additive.logistic(prefit$pi,T)[1:(G-1)],unlist(prefit$coef))
}
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
if(G==1) return(prefit)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=as.numeric(unlist(sp.data)),covar.data)
names(data)[1] <- as.character(sp.form)[2]
fmM.out <- fitmix.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var,r1)
while(fmM.out$logl==0) {
prefit <- SpeciesMix.em(sp.form,t.sp.data,t.covar.data,G,ite.max=em.steps,em.refit=1,r1)
pars <- c(additive.logistic(prefit$pi,T)[1:(G-1)],unlist(prefit$coef))
fmM.out <- fitmix.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var,r1)
}
rownames(fmM.out$tau) <- sp.name ## add names to taus
fmM.out$se <- NA
if(est.var){
fmM.out$covar <- try(solve(fmM.out$hessian))
if(class(fmM.out$covar)!="try-error"){
##colnames(fmM.out$covar) <- rownames(fmM.out$covar) <- names(fmM.out$gradient)
tmp <- sqrt(diag(fmM.out$covar))
tmp <- tmp[(G):length(tmp)]
fmM.out$se <- matrix(tmp,G,ncol(fmM.out$coef))
colnames(fmM.out$se) <- colnames(fmM.out$coef)
rownames(fmM.out$se) <- rownames(fmM.out$coef)
}
}
if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
fmM.out$formula <- sp.form
class(fmM.out) <- c("archetype","bernoulli")
fmM.out
}
"SpeciesMix.em" <-
function (sp.form,sp.data,covar.data,G=2,em.refit=1,ite.max=500, est.var = FALSE,residuals=FALSE,trace=TRUE,r1=FALSE)
{
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=unlist(sp.data),covar.data)
fmM.out <- fitMix(sp.form,data,sp,G,ite.max,trace,r1)
if(em.refit>1)
for(i in 2:em.refit){
fmM <- fitMix(sp.form,data,sp,G,ite.max,trace,r1)
if(fmM$logl>fmM.out$logl) fmM.out <- fmM
}
if(est.var){
var <- 0
t.pi <- additive.logistic(fmM.out$pi,TRUE)
parms <- c(t.pi[1:(G-1)],unlist(fmM.out$coef))
first.fit <- list(y=data[,1],x=model.matrix(sp.form,data=data))
fun.est.var <- function(x){-logLmix(x,first.fit,G,S,sp,sp.name)}
var <- solve( nH2( pt=parms, fun=fun.est.var))
colnames( var) <- rownames( var) <- names( parms)
fmM.out$covar <- var
}
if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
fmM.out$formula <- sp.form
fmM.out
}
"SpeciesMix.em.gaussian" <-
function (sp.form,sp.data,covar.data,G=2,em.refit=1,ite.max=500, est.var = FALSE,residuals=FALSE,trace=TRUE)
{
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=unlist(sp.data),covar.data)
fmM.out <- fitmix.gaussian(sp.form,data,sp,G,ite.max,trace)
if(em.refit>1)
for(i in 2:em.refit){
fmM <- fitmix.gaussian(sp.form,data,sp,G,ite.max,trace)
if(fmM$logl>fmM.out$logl) fmM.out <- fmM
}
##est.var <- FALSE
if(est.var){
var <- 0
t.pi <- additive.logistic(fmM.out$pi,TRUE)
##parms <- c(t.pi[1:(G-1)],unlist(fmM.out$coef))
parms <- c(t.pi[1:(G-1)],fmM.out$theta,unlist(fmM.out$coef),unlist(fmM.out$sp.intercept))
first.fit <- list(y=data[,1],x=model.matrix(sp.form,data=data)[,-1])
fun.est.var <- function(x){-logLmix.gaussian(x,first.fit,G,S,sp,sp.name)}
var <- solve( nH2( pt=parms, fun=fun.est.var))
colnames( var) <- rownames( var) <- names( parms)
fmM.out$covar <- var
}
residuals <- FALSE
if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
fmM.out$formula <- sp.form
fmM.out
}
"SpeciesMix.em.nbinom" <-
function (sp.form,sp.data,covar.data,G=2,em.refit=1,ite.max=500, est.var = FALSE,residuals=FALSE,trace=TRUE)
{
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=unlist(sp.data),covar.data)
fmM.out <- fitmix.nbinom(sp.form,data,sp,G,ite.max,trace)
if(em.refit>1)
for(i in 2:em.refit){
fmM <- fitmix.nbinom(sp.form,data,sp,G,ite.max,trace)
if(fmM$logl>fmM.out$logl) fmM.out <- fmM
}
##est.var <- FALSE
if(est.var){
var <- 0
t.pi <- additive.logistic(fmM.out$pi,TRUE)
##parms <- c(t.pi[1:(G-1)],unlist(fmM.out$coef))
parms <- c(t.pi[1:(G-1)],unlist(fmM.out$coef),fmM.out$sp.intercept,rep(1,S))##unlist(fmM.out$sp.intercept))
first.fit <- list(y=data[,1],x=model.matrix(sp.form,data=data)[,-1])
# first.fit$x <- cbind(sp.mat,first.fit$x)
fun.est.var <- function(x){-logLmix.nbinom(x,first.fit,G,S,sp,sp.name)}
#
deriv <- nd2(parms,fun.est.var)
var <- solve( nH2( pt=parms, fun=fun.est.var))
colnames( var) <- rownames( var) <- names( parms)
fmM.out$covar <- var
}
residuals <- FALSE
if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
fmM.out$formula <- sp.form
fmM.out
}
"SpeciesMix.em.tweedie" <-
function (sp.form,sp.data,covar.data,G=2,em.refit=1,ite.max=500, est.var = FALSE,residuals=FALSE,trace=TRUE)
{
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=unlist(sp.data),covar.data)
fmM.out <- fitmix.tweedie(sp.form,data,sp,G,ite.max,trace)
if(em.refit>1)
for(i in 2:em.refit){
fmM <- fitmix.tweedie(sp.form,data,sp,G,ite.max,trace)
if(fmM$logl>fmM.out$logl) fmM.out <- fmM
}
est.var <- FALSE
if(est.var){
var <- 0
t.pi <- additive.logistic(fmM.out$pi,TRUE)
##parms <- c(t.pi[1:(G-1)],unlist(fmM.out$coef))
parms <- c(t.pi[1:(G-1)],unlist(fmM.out$coef),unlist(fmM.out$sp.intercept),fmM.out$theta)
first.fit <- list(y=data[,1],x=model.matrix(sp.form,data=data)[,-1])
fun.est.var <- function(x){-logLmix.tweedie(x,first.fit,G,S,sp,sp.name)}
var <- solve( nH2( pt=parms, fun=fun.est.var))
colnames( var) <- rownames( var) <- names( parms)
fmM.out$covar <- var
}
residuals <- FALSE
if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
fmM.out$formula <- sp.form
fmM.out
}
"SpeciesMix.gaussian" <-
function (sp.form,sp.data,covar.data,G=2, pars=NA, em.prefit=TRUE,em.steps=4, em.refit = 1 , est.var = FALSE,residuals=FALSE,trace=TRUE)
{
t.covar.data <- covar.data
t.sp.data <- sp.data
sp.form <- update.formula(sp.form,obs~1+.)
em.steps=100
##if(em.prefit | G==1){
prefit <- SpeciesMix.em.gaussian(sp.form,sp.data,covar.data,G,ite.max=em.steps,est.var=est.var,em.refit=em.refit)
return(prefit)
##pars <- c(additive.logistic(prefit$pi,T)[1:(G-1)],unlist(prefit$coef))
##}
##S <- dim(sp.data)[2]
##if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
##n <- dim(sp.data)[1]
##sp <- rep(sp.name,each=n)
##if(G==1) return(prefit)
##var.names <- colnames(covar.data)
##covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
##names(covar.data) <- var.names
##sp.data <- as.data.frame(sp.data)
##data <- data.frame(obs=as.numeric(unlist(sp.data)),covar.data)
##names(data)[1] <- as.character(sp.form)[2]
##fmM.out <- fitmix.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var)
##while(fmM.out$logl==0) {
## prefit <- SpeciesMix.em(sp.form,t.sp.data,t.covar.data,G,ite.max=em.steps,em.refit=1)
## pars <- c(additive.logistic(prefit$pi,T)[1:(G-1)],unlist(prefit$coef))
##fmM.out <- fitmix.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var)
##}
##rownames(fmM.out$tau) <- sp.name ## add names to taus
##if(est.var){
## fmM.out$covar <- try(solve(fmM.out$hessian))
##}
## if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
##fmM.out$formula <- sp.form
##class(fmM.out) <- "archetype"
##fmM.out
}
"SpeciesMix.nbinom" <-
function (sp.form,sp.data,covar.data,G=2, pars=NA, em.prefit=TRUE,em.steps=3, em.refit = 1 , est.var = FALSE,residuals=FALSE,trace=TRUE)
{
t.covar.data <- covar.data
t.sp.data <- sp.data
sp.form <- update.formula(sp.form,obs~1+.)
#em.steps=100
if(em.prefit | G==1){
prefit <- SpeciesMix.em.nbinom(sp.form,sp.data,covar.data,G,ite.max=em.steps,est.var=FALSE,em.refit=em.refit)
##return(prefit)
pars <- prefit$pars#c(additive.logistic(prefit$pi,T)[1:(G-1)],unlist(prefit$coef))
}
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
if(G==1) return(prefit)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=as.numeric(unlist(sp.data)),covar.data)
names(data)[1] <- as.character(sp.form)[2]
sp.form <- update.formula(sp.form,obs~-1+.)
fmM.out <- fitmix.nbinom.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var)
#while(fmM.out$logl==0) {
# prefit <- SpeciesMix.em.nbinom(sp.form,sp.data,covar.data,G,ite.max=em.steps,est.var=est.var,em.refit=em.refit)
# pars <- prefit$pars
# fmM.out <- fitmix.nbinom.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var)
#}
rownames(fmM.out$tau) <- sp.name ## add names to taus
fmM.out$se <- NA
if(est.var){
fmM.out$covar <- try(solve(fmM.out$hessian))
if(class(fmM.out$covar)!="try-error"){
#colnames(fmM.out$covar) <- rownames(fmM.out$covar) <- names(fmM.out$gradient)
tmp <- sqrt(diag(fmM.out$covar))
tmp <- tmp[-1*(1:((G-1)))]
tmp <- tmp[-1*((length(tmp)-(2*S-1)):length(tmp))]
fmM.out$se <- matrix(tmp,G,ncol(fmM.out$coef))
colnames(fmM.out$se) <- colnames(fmM.out$coef)
rownames(fmM.out$se) <- rownames(fmM.out$coef)
}
##t1 <- sqrt(diag(fmM.out$covar))
##t1 <- t1[-(1:G-1)]
##fmM.out$se.coef <- t1[1:length(fmM.out$coef)]
##t1 <- t1[-(1:length(fmM.out$coef))]
##dim(fmM.out$se.coef) <- dim(fmM.out$coef)
##fmM.out$se.int <- t1[1:S]
##fmM.out$se.disp <- t1[-(1:S)]
}
if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
fmM.out$formula <- sp.form
class(fmM.out) <- c("archetype","negbin")
fmM.out
}
"SpeciesMix.tweedie" <-
function (sp.form,sp.data,covar.data,G=2, pars=NA, em.prefit=TRUE,em.steps=3, em.refit = 1 , est.var = FALSE,residuals=FALSE,trace=TRUE)
{
t.covar.data <- covar.data
t.sp.data <- sp.data
sp.form <- update.formula(sp.form,obs~1+.)
pars <- NA
if(em.prefit | G==1){
prefit <- SpeciesMix.em.tweedie(sp.form,sp.data,covar.data,G,ite.max=em.steps,est.var=est.var,em.refit=em.refit)
##return(prefit)
pars <- prefit$pars##pars <- c(additive.logistic(prefit$pi,T)[1:(G-1)],unlist(prefit$coef))
}
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
if(G==1) return(prefit)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=as.numeric(unlist(sp.data)),covar.data)
names(data)[1] <- as.character(sp.form)[2]
sp.form <- update.formula(sp.form,obs~-1+.)
fmM.out <- fitmix.tweedie.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var)
##while(fmM.out$logl==0) {
## prefit <- SpeciesMix.em(sp.form,t.sp.data,t.covar.data,G,ite.max=em.steps,em.refit=1)
## pars <- c(additive.logistic(prefit$pi,T)[1:(G-1)],unlist(prefit$coef))
##fmM.out <- fitmix.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var)
##}
rownames(fmM.out$tau) <- sp.name ## add names to taus
fmM.out$se <- NA
if(est.var){
##fmM.out$covar <- try(solve(fmM.out$hessian))
fmM.out$covar <- try(solve(fmM.out$hessian))
if(class(fmM.out$covar)!="try-error"){
colnames(fmM.out$covar) <- rownames(fmM.out$covar) <- names(fmM.out$gradient)
tmp <- sqrt(diag(fmM.out$covar))
tmp <- tmp[-1*(1:((G-1)))]
tmp <- tmp[-1*((length(tmp)-(2*S-1)):length(tmp))]
fmM.out$se <- matrix(tmp,G,ncol(fmM.out$coef))
colnames(fmM.out$se) <- colnames(fmM.out$coef)
rownames(fmM.out$se) <- rownames(fmM.out$coef)
}
}
if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
fmM.out$formula <- sp.form
class(fmM.out) <- c("archetype","tweedie")
fmM.out
}
"tglm" <-
function ( mean.form, data, wts=NULL, phi=NULL, p=NULL, inits=NULL, vcov=TRUE, residuals=TRUE, trace=1, iter.max=150)
{
if( is.null( wts))
wts <- rep( 1, nrow( data))
e<-new.env()
e$wts <- wts
environment(mean.form) <- e
temp.p <- model.frame( mean.form, data=as.data.frame( data), weights=wts)
y <- model.response( temp.p)
names( y) <- NULL
X.p <- model.matrix( mean.form, data)
offset.p <- model.offset( temp.p)
if ( is.null( offset.p))
offset.p <- rep( 0, nrow( temp.p))
wts1 <- as.vector( model.weights( temp.p))
fm1 <- NULL
if( is.null( inits)){
inits <- rep( 0, ncol( X.p))
if( trace!=0)
print( "Obtaining initial mean values from log-linear Poisson model -- this might be stupid")
abit <- 1
ystar <- round( y, digits=0)
fm1 <- glm( ystar~-1+X.p, family=poisson( link="log"), weights=wts1)
inits <- fm1$coef
}
if( is.null( phi) & length( inits)==ncol( X.p)){
if( is.null( fm1))
inits <- c( inits, 1)
else{
if( trace!=0)
print( "Obtaining initial dispersion from the smaller of the Pearson or Deviance estimator (or 25)")
if( is.null( p) & length( inits)==ncol( X.p))
ptemp <- 1.9
else
ptemp <- p
disDev <- fm1$deviance/( length( y) - length( fm1$coef))
disPear <- sum( ( wts1 * ( y-fm1$fitted)^2) / ( fm1$fitted^ptemp)) / ( length( y) - length( fm1$coef))
dis <- min( disDev, disPear, 25)
inits <- c( inits, dis)
}
}
if( is.null( p) & length( inits)==ncol( X.p) + is.null( phi))
inits <- c( inits, 1.6)
if( length( inits) != ncol( X.p) + is.null( phi) + is.null( p)) {
print( "Initial values supplied are of the wrong length -- please check")
tmp <- c( length( inits), ncol( X.p) + is.null( phi) + is.null( p))
names( tmp) <- c( "inits", "nParams")
return( tmp)
}
fmTGLM <- tglm.fit( x=X.p, y=y, wts=wts1, offset=offset.p, inits=inits, phi=phi, p=p, vcov=vcov, residuals=residuals, trace=trace, iter.max=iter.max)
return( fmTGLM)
}
"tglm.fit" <-
function ( x, y, wts=NULL, offset=rep( 0, length( y)), inits, phi=NULL, p=NULL, vcov=TRUE, residuals=TRUE, trace=1, iter.max=150)
{
if( trace!=0){
print( "Estimating parameters")
if( is.null( phi) & is.null( p))
cat("iter:", "-logl", colnames(x), "phi", "p", "\n", sep = "\t")
if( is.null( phi) & !is.null( p))
cat("iter:", "-logl", colnames(x), "phi", "\n", sep = "\t")
if( !is.null( phi) & is.null( p))
cat("iter:", "-logl", colnames(x), "p", "\n", sep = "\t")
if( !is.null( phi) & !is.null( p))
cat("iter:", "-logl", colnames(x), "\n", sep = "\t")
}
eps <- 1e-5
my.lower <- rep( -Inf, ncol( x))
my.upper <- rep( Inf, ncol( x))
if( is.null( phi))
{my.lower <- c( my.lower, eps); my.upper <- c( my.upper, Inf)}
if( is.null( p))
{my.lower <- c( my.lower, 1+eps); my.upper <- c( my.upper, 2-eps)}
fm <- nlminb( start=inits, objective=ldTweedie.lp, gradient=ldTweedie.lp.deriv, hessian=NULL, lower=my.lower, upper=my.upper, y=y, X.p=x, offsetty=offset, phi=phi, p=p, control=list(trace=trace, iter.max=iter.max), wts=wts)
parms <- fm$par
tmp <- colnames( x)
if( !is.null( phi) & !is.null( p))
tmp <- tmp
if( !is.null( phi) & is.null( p))
tmp <- c( tmp, "p")
if( is.null( phi) & !is.null( p))
tmp <- c( tmp, "phi")
if( is.null( phi) & is.null( p))
tmp <- c( tmp, "phi", "p")
names( parms) <- tmp
if( vcov){
if( trace!=0)
print( "Calculating variance matrix of estimates")
vcovar <- nd2(x0=parms, f=ldTweedie.lp.deriv, y=y, X.p=x, offsetty=offset, phi=phi, p=p, wts=wts)
vcovar <- 0.5 * ( vcovar + t( vcovar))
vcovar <- solve( vcovar)
rownames( vcovar) <- colnames( vcovar) <- names( parms)
}
else{
if( trace!=0)
print( "Not calculating variance matrix of estimates")
vcovar <- NULL
}
scores <- -ldTweedie.lp.deriv( parms=parms, y=y, X.p=x, offsetty=offset, phi=phi, p=p, wts=wts)
if( trace !=0)
print( "Calculating means")
mu <- exp( x %*% parms[1:ncol( x)] + offset)
if( residuals){
if( trace!=0)
print( "Calculating quantile residuals")
if( is.null( phi)){
phi1 <- parms[ncol( x)+1]
if( is.null( p))
p1 <- parms[ncol(x)+2]
else
p1 <- p
}
else{
phi1 <- phi
if( is.null( p))
p1 <- parms[ncol( x)+1]
else
p1 <- p
}
resids <- matrix( rep( pTweedie( y, mu, phi1, p1), 2), ncol=2)
resids[y==0,1] <- 0.5 * dTweedie( y[y==0], mu[y==0], phi1, p1, LOG=FALSE)
nzero <- sum( y==0)
resids[y==0,2] <- runif( nzero, min=rep( 0, nzero), max=2*resids[y==0,1])
resids <- qnorm( resids)
colnames( resids) <- c("expect","random")
}
else{
if( trace!=0)
print( "Not calculating quantile residuals")
resids <- NULL
}
if( trace!=0)
print( "Done")
ICadj <- 0
if( !is.null( phi))
ICadj <- ICadj+1
if( !is.null( p))
ICadj <- ICadj+1
AIC <- -2*(-fm$objective) + 2*(length( parms)-ICadj)
BIC <- -2*(-fm$objective) + log( nrow( x))*(length( parms)-ICadj)
res <- list( coef=parms, logl=-fm$objective, scores=scores, vcov=vcovar, conv=fm$convergence, message=fm$message, niter=fm$iterations, evals=fm$evaluations, call=match.call(), fitted=mu, residuals=resids, AIC=AIC, BIC=BIC)
class( res) <- "tglm"
return( res)
}
"weighted.glm" <-
function (g,first.fit,tau,n,fmM,sp)
{
dat.tau <- rep(tau[,g],each=n)
##lpre <- first.fit$x%*%fmM[[g]]$coef
##f.mix <- glm.fit(x=first.fit$x,y=first.fit$y,weights=dat.tau[,g],family=binomial(),etastart=fmM[[g]]$linear.predictors)
f.mix <- glm.fit(x=first.fit$x,y=first.fit$y,weights=dat.tau,family=binomial(),start=fmM[[g]]$coef)
return(list(coef=f.mix$coef))
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
##list(linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
}
"weighted.glm.gaussian" <-
function (g,first.fit,tau,n,fmM,sp)
{
dat.tau <- rep(tau[,g],each=n)
sp.name <- unique(sp)
X <- first.fit$x
sp.mat <- matrix(0,dim(X)[1],length(sp.name))
for(i in 1:length(sp.name)){
sp.mat[sp==sp.name[i],i] <- 1
}
X <- cbind(sp.mat,X)
f.mix <- glm.fit(x=X,y=first.fit$y,weights=dat.tau,family=gaussian())
#f.mix$theta <- 1
## f.mix <- glm.fit.nbinom(x=X,y=first.fit$y,weights=dat.tau)
sp.intercept <- f.mix$coef[1:length(sp.name)]
sp.intercept[is.na(sp.intercept)] <- 0
return(list(coef=f.mix$coef[-1:-length(sp.name)],theta=sqrt(f.mix$deviance/f.mix$df.residual),sp.intercept=sp.intercept,fitted=f.mix$fitted))#,lpre=f.mix$linear.predictors))
}
"weighted.glm.nbinom" <-
function (g,first.fit,tau,n,fmM,sp)
{
dat.tau <- rep(tau[,g],each=n)
sp.name <- unique(sp)
X <- first.fit$x
sp.mat <- matrix(0,dim(X)[1],length(sp.name))
for(i in 1:length(sp.name)){
sp.mat[sp==sp.name[i],i] <- 1
}
X <- cbind(sp.mat,X)
#f.mix <- glm.fit(x=X,y=first.fit$y,weights=dat.tau,family=poisson())
#f.mix$theta <- 1
f.mix <- glm.fit.nbinom(x=X,y=first.fit$y,offset=first.fit$offset,weights=dat.tau)
sp.intercept <- f.mix$coef[1:length(sp.name)]
sp.intercept[is.na(sp.intercept)] <- 0
return(list(coef=f.mix$coef[-1:-length(sp.name)],theta=f.mix$theta,sp.intercept=sp.intercept,fitted=f.mix$fitted))#,lpre=f.mix$linear.predictors))
}
"weighted.glm.tweedie" <-
function (g, first.fit, tau, n, fmM, sp)
{
dat.tau <- rep(tau[, g], each = n)
sp.int.offset <- rep(fmM[[g]]$sp.intercept, each = n)+first.fit$offset
sp.name <- unique(sp)
X <- first.fit$x
f.mix <- tglm.fit(x = X, y = first.fit$y, wts = dat.tau,
offset = sp.int.offset, vcov = FALSE, residuals = FALSE,
trace = 0, inits = fmM[[g]]$coef, phi = fmM[[g]]$phi,
p = 1.6)
return(list(coef = f.mix$coef, phi = fmM[[g]]$phi, p = 1.6,
sp.intercept = fmM[[g]]$sp.intercept, fitted = f.mix$fitted))
}
Try the SpeciesMix package in your browser
Any scripts or data that you put into this service are public.
SpeciesMix documentation built on May 2, 2019, 4:22 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# This is package SpeciesMix
"additive.logistic" <-
function (x,inv=FALSE)
{
if(inv){
x <- log(x/x[length(x)])
return(x)
}
x.t <- exp(x)
x.t <- x.t/(1+sum(x.t))
x.t[length(x.t)+1] <- 1-sum(x.t)
return(x.t)
}
"apply.glm" <-
function (i,form,datsp,tau,n)
{
dat.tau <- rep(tau[,i],each=n)
x <- model.matrix(as.formula(form),data=datsp)
y <- datsp$obs
##dat.tau <- get("dat.tau")
##datsp <- get("datsp")
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
f.mix <- glm.fit(x,y,dat.tau,family=binomial())
return(list(coef=f.mix$coef))
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="binomial",x=T,y=T)
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
}
"apply.glm.gaussian" <-
function (i,form,datsp,tau,n)
{
dat.tau <- rep(tau[,i],each=n)
x <- model.matrix(as.formula(form),data=datsp)
y <- datsp$obs
environment(form) <- environment()
##dat.tau <- get("dat.tau")
##datsp <- get("datsp")
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
##f.mix <- glm.nb(as.formula(form),data=datsp,weights=dat.tau)
f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="gaussian")
##f.mix <- glm.nbinom(as.formula(form),data=datsp,weights=dat.tau)
sp.int <- rep(f.mix$coef[1],dim(tau)[1])
##return(list(coef=f.mix$coef[-1],theta=f.mix$theta))
return(list(coef=f.mix$coef[-1],theta=sqrt(f.mix$deviance/f.mix$df.residual),sp.intercept=sp.int))
##return(list(coef=f.mix$coef,theta=f.mix$theta))
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="binomial",x=T,y=T)
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
}
"apply.glm.nbinom" <-
function (i,form,datsp,tau,n)
{
dat.tau <- rep(tau[,i],each=n)
x <- model.matrix(as.formula(form),data=datsp)
y <- datsp$obs
environment(form) <- environment()
##dat.tau <- get("dat.tau")
##datsp <- get("datsp")
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
##f.mix <- glm.nb(as.formula(form),data=datsp,weights=dat.tau)
#f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="poisson")
f.mix <- glm.nbinom(as.formula(form),data=datsp,weights=dat.tau)
sp.int <- rep(f.mix$coef[1],dim(tau)[1])
##return(list(coef=f.mix$coef[-1],theta=f.mix$theta))
return(list(coef=f.mix$coef[-1],theta=f.mix$theta,sp.intercept=sp.int))
##return(list(coef=f.mix$coef,theta=f.mix$theta))
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="binomial",x=T,y=T)
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
}
"apply.glm.tweedie" <-
function (i,form,datsp,tau,n)
{
dat.tau <- rep(tau[,i],each=n)
x <- model.matrix(as.formula(form),data=datsp)
y <- datsp$obs
environment(form) <- environment()
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
f.mix <- tglm(as.formula(form),data=datsp,wts=dat.tau,vcov=FALSE,residuals=FALSE,trace=0)
sp.int <- rep(f.mix$coef[1],dim(tau)[1])
return(list(coef=f.mix$coef[c(-1,-(length(f.mix$coef)-1),-(length(f.mix$coef)))],phi=f.mix$coef["phi"],p=f.mix$coef["p"],sp.intercept=sp.int))
}
"artificial.data" <-
function (formula, data, theta, S, dist = "bernoulli")
{
X <- model.matrix(formula, data)
out <- matrix(0, dim(X)[1], S)
k <- dim(theta)[1]
sp.int <- rep(0, S)
group <- rep(0, S)
for (s in 1:S) {
g <- ceiling(runif(1) * k)
if (dist == "bernoulli") {
lgtp <- X %*% theta[g, ]
p <- exp(lgtp)/(1 + exp(lgtp))
out[, s] <- rbinom(dim(X)[1], 1, p)
}
if (dist == "negbin") {
tmp <- rep(1e+05, dim(X)[1])
while (max(tmp, na.rm = T) > 5000 | sum(tmp) < 100) {
theta[g, 1] <- runif(1, -15, 5)
sp.int[s] <- theta[g, 1]
lgtp <- X %*% theta[g, ]
p <- exp(lgtp)
tmp <- rnbinom(dim(X)[1], mu = p, size = 1)
}
out[, s] <- tmp
}
if (dist == "tweedie") {
tmp <- rep(6e+05, dim(X)[1])
while (max(tmp, na.rm = T) > 5e+05 | sum(tmp) < 100) {
theta[g, 1] <- runif(1, -15, 5)
theta[g, 1] <- runif(1, 1, 5)
sp.int[s] <- theta[g, 1]
lgtp <- X %*% theta[g, ]
p <- exp(lgtp)
tmp <- rTweedie(dim(X)[1], mu = p, phi = 2, p = 1.6)
}
out[, s] <- tmp
}
if (dist == "gaussian") {
tmp <- rep(1e+05, dim(X)[1])
while (max(tmp, na.rm = T) > 50000 | sum(tmp) < 100) {
theta[g, 1] <- runif(1, 100, 500)
lgtp <- X %*% theta[g, ]
p <- (lgtp)
tmp <- rnorm(dim(X)[1], mean = p, sd = 1)
}
out[, s] <- tmp
}
group[s] <- g
}
pi <- tapply(group, group, length)/S
if (dist == "negbin")
return(list(pa = out, group = group, pi = pi, sp.int = sp.int))
if (dist == "tweedie")
return(list(pa = out, group = group, pi = pi, sp.int = sp.int))
list(pa = out, group = group, pi = pi)
}
"clusterSelect" <-
function (sp.form,sp.data,covar.data,G=1:10,em.prefit=TRUE, em.steps=4 ,em.refit=3,est.var=FALSE,trace=TRUE)
{
my.fun <- function(g,form,sp.data,covar.data){
cat("Fitting group",g,"\n")
try(SpeciesMix(form,sp.data,covar.data,g,em.prefit=em.prefit,em.steps=em.steps,em.refit=em.refit,est.var=est.var,trace=trace))
}
# if(mc){ out <- mclapply(G,my.fun,form,dat,mc.preschedule = FALSE, mc.set.seed = TRUE, mc.silent = FALSE, mc.cores = set.cores)} else
{ out <- lapply(G,my.fun,sp.form,sp.data,covar.data)}
aic <- rep(0,length(G))
bic <- rep(0,length(G))
fm <- list()
for(i in 1:length(G))
if(!is.atomic(out[[i]])){
aic[i] <- out[[i]]$aic
bic[i] <- out[[i]]$bic
fm[[i]] <- list(logl=out[[i]]$logl,coef=out[[i]]$coef,tau=out[[i]]$tau,pi=out[[i]]$pi,covar=out[[i]]$covar)
}
return(list(aic=aic,bic=bic,fm=fm))
}
"create.starting.values" <-
function (S,G,n,form,datsp)
{
##apply.glm <- function(i,form,datsp,tau,n,dat.tau){
## dat.tau <- rep(tau[,i],each=n)
## dat.tau <- get("dat.tau")
## datsp <- get("datsp")
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
## f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="binomial",x=T,y=T)
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
## list(coef=f.mix$coef)
##}
environment(form) <- environment()
tau <- matrix(runif(S*G),S,G)
tau <- (tau/rowSums(tau))
fmM <- list()
for(i in 1:G){
##dat.tau[,i] <- rep(tau[,i],each=n)
pi[i] <- sum(tau[,i])/S
}
##dat.tau <- rep(0,dim(datsp)[1])
fmM <- lapply(1:G,apply.glm,form,datsp,tau,n)
first.fit <- list(x=model.matrix(as.formula(form),data=datsp),y=datsp$obs,formula=form)
##return(list(pi=pi,fmM=fmM,tau=tau,dat.tau=dat.tau,first.fit=first.fit))
return(list(pi=pi,fmM=fmM,tau=tau,first.fit=first.fit))
}
"create.starting.values.gaussian" <-
function (S,G,n,form,datsp,mc=FALSE,set.cores=2)
{
##apply.glm <- function(i,form,datsp,tau,n,dat.tau){
## dat.tau <- rep(tau[,i],each=n)
## dat.tau <- get("dat.tau")
## datsp <- get("datsp")
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
## f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="binomial",x=T,y=T)
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
## list(coef=f.mix$coef)
##}
environment(form) <- environment()
tau <- matrix(runif(S*G),S,G)
tau <- (tau/rowSums(tau))
fmM <- list()
for(i in 1:G){
##dat.tau[,i] <- rep(tau[,i],each=n)
pi[i] <- sum(tau[,i])/S
}
##dat.tau <- rep(0,dim(datsp)[1])
fmM <- lapply(1:G,apply.glm.gaussian,form,datsp,tau,n)
first.fit <- list(x=model.matrix(as.formula(form),data=datsp)[,-1],y=datsp$obs,formula=form)
##return(list(pi=pi,fmM=fmM,tau=tau,dat.tau=dat.tau,first.fit=first.fit))
return(list(pi=pi,fmM=fmM,tau=tau,first.fit=first.fit))
}
"create.starting.values.nbinom" <-
function (S,G,n,form,datsp,mc=FALSE,set.cores=2)
{
##apply.glm <- function(i,form,datsp,tau,n,dat.tau){
## dat.tau <- rep(tau[,i],each=n)
## dat.tau <- get("dat.tau")
## datsp <- get("datsp")
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
## f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="binomial",x=T,y=T)
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
## list(coef=f.mix$coef)
##}
environment(form) <- environment()
tau <- matrix(runif(S*G),S,G)
tau <- (tau/rowSums(tau))
fmM <- list()
for(i in 1:G){
##dat.tau[,i] <- rep(tau[,i],each=n)
pi[i] <- sum(tau[,i])/S
}
##dat.tau <- rep(0,dim(datsp)[1])
fmM <- lapply(1:G,apply.glm.nbinom,form,datsp,tau,n)
offset <- model.frame(as.formula(form),data=datsp)
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,length(datsp$obs))
first.fit <- list(x=model.matrix(as.formula(form),data=datsp)[,-1],y=datsp$obs,offset=offset,formula=form)
##return(list(pi=pi,fmM=fmM,tau=tau,dat.tau=dat.tau,first.fit=first.fit))
return(list(pi=pi,fmM=fmM,tau=tau,first.fit=first.fit))
}
"create.starting.values.nbinom.kmeans" <-
function (S, G, n, form, datsp, tol = 0.1)
{
MM <- model.matrix(form, datsp)
offset <- model.frame(form,datsp)
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,nrow(MM))
first.fit <- list(x = model.matrix(as.formula(form), data = datsp)[,
-1], y = datsp$obs, offset=offset, formula = form)
if(class(first.fit$x)=="numeric") {first.fit$x <- matrix(first.fit$x,length(first.fit$x),1) }#deal with model matrix returning a vector for covar == 1
if (tol < 0 || tol >= 1)
stop("Minimum Prevalence % must be between 0 and 1")
sp.name <- 1:S
sp <- rep(sp.name, each = n)
starting.fitem <- list(intercepts = rep(0, S), alpha = rep(0,
S))
all.betas <- matrix(0, nrow = S, ncol = ncol(MM))
colnames(all.betas) <- colnames(MM)
for (j in 1:S) {
fit <- glm.nbinom(form, datsp[((j - 1) * n):(j * n -
1), ], est.var = FALSE)
starting.fitem$sp.intercepts[j] <- fit$coef[1]
all.betas[j, ] <- fit$coef
starting.fitem$theta[j] = fit$theta
}
all.betas[, 1] <- 0
cat("Clustering...\n")
fmmvnorm <- kmeans(x = all.betas, centers = G, iter.max = 100,
nstart = 50)
starting.fitem$coef <- fmmvnorm$centers
fmM <- list()
for (i in 1:G) {
B <- matrix(rep(fmmvnorm$centers[i, ], nrow(datsp)),
nrow(datsp), ncol(fmmvnorm$centers), byrow = T)
B[, 1] <- rep(starting.fitem$sp.intercepts, each = n)
fitted <- exp(rowSums(MM * B)+offset)
fmM[[i]] <- list(coef = fmmvnorm$centers[i, 2:ncol(fmmvnorm$centers)],
theta = mean(starting.fitem$theta[fmmvnorm$cluster ==
i]), sp.intercept = starting.fitem$sp.intercepts,
fitted = fitted)
}
tau <- matrix(0, S, G)
pi <- rep(1/G, G)
pi <- runif(G, 0.2, 0.8)
pi <- pi/sum(pi)
est.tau <- lapply(1:S, estimate.pi.nbinom, sp, sp.name, datsp,
fmM, pi, G, first.fit)
max.newTau <- 0.8
alpha <- (1 - max.newTau * G)/(max.newTau * (2 - G) - 1)
for (j in 1:S) {
newTau <- (2 * alpha * est.tau[[j]]$tau - alpha + 1)/(2 *
alpha - alpha * G + G)
tau[j, ] <- newTau
}
for (i in 1:G) {
pi[i] <- sum(tau[, i])/S
}
return(list(pi = pi, fmM = fmM, tau = tau, first.fit = first.fit))
}
"create.starting.values.tweedie" <-
function (S,G,n,form,datsp)
{
##apply.glm <- function(i,form,datsp,tau,n,dat.tau){
## dat.tau <- rep(tau[,i],each=n)
## dat.tau <- get("dat.tau")
## datsp <- get("datsp")
##f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau[,i],family="binomial",x=T,y=T)
## f.mix <- glm(as.formula(form),data=datsp,weights=dat.tau,family="binomial",x=T,y=T)
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
## list(coef=f.mix$coef)
##}
environment(form) <- environment()
tau <- matrix(runif(S*G),S,G)
tau <- (tau/rowSums(tau))
fmM <- list()
for(i in 1:G){
##dat.tau[,i] <- rep(tau[,i],each=n)
pi[i] <- sum(tau[,i])/S
}
##dat.tau <- rep(0,dim(datsp)[1])
offset <- model.frame(form,datsp)
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,length(datsp$obs))
fmM <- lapply(1:G,apply.glm.tweedie,form,datsp,tau,n)
first.fit <- list(x=model.matrix(as.formula(form),data=datsp)[,-1],y=datsp$obs,formula=form)
##return(list(pi=pi,fmM=fmM,tau=tau,dat.tau=dat.tau,first.fit=first.fit))
return(list(pi=pi,fmM=fmM,tau=tau,first.fit=first.fit))
}
"create.starting.values.tweedie.kmeans" <-
function (S, G, n, form, datsp, tol = 0.1)
{
MM <- model.matrix(form, datsp)
offset <- model.frame(form,datsp)
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,nrow(MM))
first.fit <- list(x = model.matrix(as.formula(form), data = datsp)[,
-1], y = datsp$obs, offset=offset,formula = form)
if(class(first.fit$x)=="numeric") {first.fit$x <- matrix(first.fit$x,length(first.fit$x),1) }#deal with model matrix returning a vector for covar == 1
if (tol < 0 || tol >= 1)
stop("Minimum Prevalence % must be between 0 and 1")
sp.name <- 1:S
sp <- rep(sp.name, each = n)
starting.fitem <- list(intercepts = rep(0, S), alpha = rep(0,
S))
all.betas <- matrix(0, nrow = S, ncol = ncol(MM))
colnames(all.betas) <- colnames(MM)
for (j in 1:S) {
fit <- tglm(form, datsp[((j - 1) * n):(j * n - 1), ],
vcov = FALSE, residuals = FALSE, trace = 0, p = 1.6)
starting.fitem$sp.intercepts[j] <- fit$coef[1]
all.betas[j, ] <- fit$coef[-length(fit$coef)]
starting.fitem$phi[j] = fit$coef["phi"]
}
all.betas[, 1] <- 0
cat("Clustering...\n")
fmmvnorm <- kmeans(x = all.betas, centers = G, iter.max = 100,
nstart = 50)
starting.fitem$coef <- fmmvnorm$centers
fmM <- list()
for (i in 1:G) {
B <- matrix(rep(fmmvnorm$centers[i, ], nrow(datsp)),
nrow(datsp), ncol(fmmvnorm$centers), byrow = T)
B[, 1] <- rep(starting.fitem$sp.intercepts, each = n)
fitted <- exp(rowSums(MM * B)+offset)
fmM[[i]] <- list(coef = fmmvnorm$centers[i, 2:ncol(fmmvnorm$centers)],
phi = mean(starting.fitem$phi[fmmvnorm$cluster ==
i]), p = 1.6, sp.intercept = starting.fitem$sp.intercepts,
fitted = fitted)
}
tau <- matrix(0, S, G)
pi <- rep(1/G, G)
pi <- runif(G, 0.2, 0.8)
pi <- pi/sum(pi)
est.tau <- lapply(1:S, estimate.pi.tweedie, sp, sp.name,
datsp, fmM, pi, G, first.fit)
max.newTau <- 0.8
alpha <- (1 - max.newTau * G)/(max.newTau * (2 - G) - 1)
for (j in 1:S) {
newTau <- (2 * alpha * est.tau[[j]]$tau - alpha + 1)/(2 *
alpha - alpha * G + G)
tau[j, ] <- newTau
}
for (i in 1:G) {
pi[i] <- sum(tau[, i])/S
}
return(list(pi = pi, fmM = fmM, tau = tau, first.fit = first.fit))
}
"distr.binom" <-
function( p){
nobs <- length( p)
new.dist <- old.dist <- rep( 0, nobs+1)
old.dist[1] <- 1-p[1]
old.dist[2] <- p[1]
for( ii in 2:nobs){
new.dist[1] <- old.dist[1]*(1-p[ii])
for( jj in 2:ii)
new.dist[jj] <- old.dist[jj-1]*p[ii] + old.dist[jj]*(1-p[ii])
new.dist[ii+1] <- old.dist[ii]*p[ii]
old.dist <- new.dist
}
return( new.dist)
}
"dPoisGam" <-
function ( y, lambda, mu.Z, alpha, LOG=TRUE)
{
#function to calculate Random sum (Tweedie) densities.
#y is the value of the r.v. Can be a vector
#mu.N is the mean of the Poisson summing r.v. Can be a vector of length(y)
#mu.Z is the mean of the Gamma rv Can be a vector of length(y)
#alpha is the `other' parameter of the gamma distribution s.t. var = ( mu.Z^2)/alpha Can be a vector of length(y)
#If mu.N, mu.Z or alpha are scalare but y isn't then they will be used for all y. If lengths mis-match then error
#LOG=TRUE gives the density on the log scale
#do.checks=TRUE checks the input vectors for compatability and gives errors / changes them as appropriate.
#do.checks=FALSE doesn't check and relies on the user to have things right. If not right then catastrophic failure may occur.
# if( any( is.null( c( y, mu.N, mu.Z, alpha)))){
# print( "Error: null input values -- please check. Null values are:")
# tmp <- double( is.null( c( y, mu.N, mu.Z, alpha)))
# names( tmp) <- c( "y", "mu.N","mu.Z","alpha")
# print( tmp)
# print( "Exitting")
# return()
# }
mu.N <- lambda
if( !all( is.element( c( length( mu.N), length( mu.Z), length( alpha)), c( length( y), 1)))){
print( "Error: length of parameter vectors does not match length of random variable vector")
print( "Exitting")
return()
}
if( length( mu.N) != length( y))
mu.N <- rep( mu.N, length( y))
if( length( mu.Z) != length( y))
mu.Z <- rep( mu.Z, length( y))
if( length( alpha) != length( y))
alpha <- rep( alpha, length( y))
res <- .Call( "dTweedie", as.numeric( y), as.numeric( mu.N), as.numeric( mu.Z), as.numeric( alpha), as.integer( LOG),PACKAGE="SpeciesMix")
return( res)
}
"dPoisGamDerivs" <-
function ( y=NULL, lambda=NULL, mu.Z=NULL, alpha=NULL, do.checks=TRUE)
{
#function to calculate Random sum (Tweedie) densities.
#y is the value of the r.v. Can be a vector
#mu.N is the mean of the Poisson summing r.v. Can be a vector of length(y)
#mu.Z is the mean of the Gamma rv Can be a vector of length(y)
#alpha is the `other' parameter of the gamma distribution s.t. var = ( mu.Z^2)/alpha Can be a vector of length(y)
#If mu.N, mu.Z or alpha are scalare but y isn't then they will be used for all y. If lengths mis-match then error
#LOG=TRUE gives the density on the log scale
#do.checks=TRUE checks the input vectors for compatability and gives errors / changes them as appropriate.
#do.checks=FALSE doesn't check and relies on the user to have things right. If not right then catastrophic failure may occur.
mu.N <- lambda
if( do.checks){
if( any( is.null( c( y, mu.N, mu.Z, alpha)))){
print( "Error: null input values -- please check. Null values are:")
tmp <- double( is.null( c( y, mu.N, mu.Z, alpha)))
names( tmp) <- c( "y", "mu.N","mu.Z","alpha")
print( tmp)
print( "Exitting")
return()
}
if( !all( is.element( c( length( mu.N), length( mu.Z), length( alpha)), c( length( y), 1)))){
print( "Error: length of parameter vectors does not match length of random variable vector")
print( "Exitting")
}
if( length( mu.N) != length( y))
mu.N <- rep( mu.N, length( y))
if( length( mu.Z) != length( y))
mu.Z <- rep( mu.Z, length( y))
if( length( alpha) != length( y))
alpha <- rep( alpha, length( y))
}
res <- .Call( "dTweedieDeriv", as.numeric( y), as.numeric( mu.N), as.numeric( mu.Z), as.numeric( alpha),PACKAGE="SpeciesMix")
colnames( res) <- c("lambda","mu.Z","alpha")
return( res)
}
"dTweedie" <-
function ( y, mu, phi, p, LOG=TRUE)
{
lambda <- ( mu^( 2-p)) / ( phi*(2-p))
alpha <- ( 2-p) / ( p-1)
tau <- phi*(p-1)*mu^(p-1)
mu.Z <- alpha * tau
dens <- dPoisGam( y, lambda, mu.Z, alpha, LOG)
return( dens)
}
"estimate.pi" <-
function (j,sp,spname,datsp,fmM,pi,G,first.fit)
{
tmp.like <- rep(0,G)
tau <- rep(0,G)
sel.sp <- which(sp==spname[j])
link.fun <- make.link("logit")
for(i in 1:G) {
if(length(fmM[[i]]$coef)==1){lpre <- link.fun$linkinv(first.fit$x[sel.sp,]*fmM[[i]]$coef)
}else{ lpre <- link.fun$linkinv(first.fit$x[sel.sp,]%*%fmM[[i]]$coef)}
##lpre <- link.fun$linkinv(first.fit$x[sel.sp,]%*%fmM[[i]]$coef)
## tmp.like[i] <- sum(dbinom(datsp$obs[sel.sp],1,fmM[[i]]$fitted[sel.sp],log=T))
obs <- datsp$obs[sel.sp]
lpre[obs==0]<- 1- lpre[obs==0]
tmp.like[i] <- sum(log(lpre))
}
eps <- max(tmp.like)
sum.like <- (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
for(i in 1:G) {
tau[i] <- exp((log(pi[i]) + tmp.like[i]) - sum.like )
}
return(list(tau=tau,sum.like=sum.like))
}
"estimate.pi.gaussian" <-
function (j,sp,spname,datsp,fmM,pi,G,first.fit)
{
tmp.like <- rep(0,G)
tau <- rep(0,G)
sel.sp <- which(sp==spname[j])
for(i in 1:G) {
##lpre <- dnorm(first.fit$y[sel.sp],mean=fmM[[i]]$fitted[sel.sp],sd=fmM[[i]]$theta,log=TRUE)
lpre <- dnorm(first.fit$y[sel.sp],mean=fmM[[i]]$fitted[sel.sp],sd=1,log=TRUE)
##lpre <- dnbinom(first.fit$y[sel.sp],mu=exp(cbind(1,first.fit$x[sel.sp,])%*%c(fmM[[i]]$sp.intercept[j],fmM[[i]]$coef)),size=fmM[[i]]$theta,log=TRUE)
##lpre <- dpois(first.fit$y[sel.sp],exp(cbind(1,first.fit$x[sel.sp,])%*%c(fmM[[i]]$sp.intercept[j],fmM[[i]]$coef)),log=TRUE)
#lpre <- dpois(first.fit$y[sel.sp],exp(first.fit$x[sel.sp,]%*%fmM[[i]]$coef),log=TRUE)
tmp.like[i] <- sum(lpre)
}
eps <- max(tmp.like)
sum.like <- (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
for(i in 1:G) {
tau[i] <- exp((log(pi[i]) + tmp.like[i]) - sum.like )
}
return(list(tau=tau,sum.like=sum.like))
}
"estimate.pi.nbinom" <-
function (j,sp,spname,datsp,fmM,pi,G,first.fit)
{
tmp.like <- rep(0,G)
tau <- rep(0,G)
sel.sp <- which(sp==spname[j])
for(i in 1:G) {
lpre <- dnbinom(first.fit$y[sel.sp],mu=fmM[[i]]$fitted[sel.sp],size=fmM[[i]]$theta,log=TRUE)
##lpre <- dnbinom(first.fit$y[sel.sp],mu=exp(cbind(1,first.fit$x[sel.sp,])%*%c(fmM[[i]]$sp.intercept[j],fmM[[i]]$coef)),size=fmM[[i]]$theta,log=TRUE)
##lpre <- dpois(first.fit$y[sel.sp],exp(cbind(1,first.fit$x[sel.sp,])%*%c(fmM[[i]]$sp.intercept[j],fmM[[i]]$coef)),log=TRUE)
#lpre <- dpois(first.fit$y[sel.sp],exp(first.fit$x[sel.sp,]%*%fmM[[i]]$coef),log=TRUE)
tmp.like[i] <- sum(lpre)
}
eps <- max(tmp.like)
sum.like <- (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
for(i in 1:G) {
tau[i] <- exp((log(pi[i]) + tmp.like[i]) - sum.like )
}
return(list(tau=tau,sum.like=sum.like))
}
"estimate.pi.tweedie" <-
function (j, sp, spname, datsp, fmM, pi, G, first.fit)
{
tmp.like <- rep(0, G)
tau <- rep(0, G)
sel.sp <- which(sp == spname[j])
for (i in 1:G) {
lpre <- dTweedie(first.fit$y[sel.sp], mu = fmM[[i]]$fitted[sel.sp],
phi = fmM[[i]]$phi, p = fmM[[i]]$p, LOG = TRUE)
tmp.like[i] <- sum(lpre)
}
eps <- max(tmp.like)
sum.like <- (log(sum(pi * exp((tmp.like) - (eps)))) + (eps))
for (i in 1:G) {
tau[i] <- exp((log(pi[i]) + tmp.like[i]) - sum.like)
}
return(list(tau = tau, sum.like = sum.like))
}
"fitMix" <-
function (form,datsp,sp,G=2,ite.max=500,trace=TRUE,full.model=FALSE,r1=FALSE)
{
## dat2 has colums obs,sp
##
temp.warn <- getOption( "warn")
options( warn=-1)
sp.name <- unique(sp)
S <- length(unique(sp))
n <- length(which(sp==sp.name[1]))
cat("Fitting Group",G,"\n")
if(trace) cat("Iteration | LogL \n")
##dat.tau <- data.frame(matrix(0,dim(datsp)[1],G))
dat.tau <- 0
##dat <- data.frame(datsp,dat.tau)
pi <- rep(0,G)
ite <- 1
logL <- -99999999
old.logL <- -88888888
## set up initial GLM
t1 <- create.starting.values(S,G,n,form,datsp)
pi <- t1$pi;fmM <- t1$fmM;tau <- t1$tau;dat.tau <- t1$dat.tau;first.fit <- t1$first.fit
while(abs(logL-old.logL) > 0.0001 & ite<=ite.max){
old.logL <- logL
for(i in 1:G){
##dat.tau[,i] <- rep(tau[,i],each=n)
pi[i] <- sum(tau[,i])/S
}
if(any(pi==0)) { ## occasionally with complicated models the random starting values result in a pi[i]==0; so restart with new random starts
cat("pi has gone to zero - restarting fitting \n")
t1 <- create.starting.values(S,G,n,form,datsp)
pi <- t1$pi;fmM <- t1$fmM;tau <- t1$tau;dat.tau <- t1$dat.tau;first.fit <- t1$first.fit
ite <- 1
}
fmM <- lapply(1:G,weighted.glm,first.fit,tau,n,fmM,sp)
logL <- 0
tmp.like <- matrix(0,S,G)
est.tau <- lapply(1:S,estimate.pi,sp,sp.name,datsp,fmM,pi,G,first.fit)
for(j in 1:S){
if(is.atomic(est.tau[[j]])){ print (est.tau[[j]])} else
{
tau[j,] <- est.tau[[j]]$tau
logL <- logL+est.tau[[j]]$sum.like
}
}
if(trace) cat(ite," | ",logL,"\n")
ite <- ite+1
}
fm.out <- data.frame(matrix(0,G,length(fmM[[1]]$coef)))
names(fm.out) <- names(fmM[[1]]$coef)
tau <- data.frame(tau)
names(tau) <- paste("grp.",1:G,sep="")
EN <- -sum(unlist(tau)*log(unlist(tau)))
d <- length(unlist(fm.out)) + length(tau)-1
for(i in 1:G) {
fm.out[i,] <- fmM[[i]]$coef
##dat.tau[,i] <- rep(tau[,i],each=n)
}
names(pi) <- paste("G",1:G,sep=".")
t.pi <- additive.logistic(pi,TRUE)
parms <- c(t.pi[1:(G-1)],unlist(fm.out))
logL.full <- logL
logL <- logLmix(parms,first.fit,G,S,sp,sp.name)
options(warn=temp.warn)
if(full.model) return(list(logl=logL,aic = -2*logL + 2*d,tau=round(tau,4),pi=pi,bic=-2*logL + log(S)*d,ICL= -2*logL + log(S)*d +2*EN,coef=fm.out,fmM=fmM,model.tau=dat.tau,covar=0,aic.full= -2*logL.full + 2*d,bic.full= -2*logL.full + log(S)*d))
return(list(logl=logL,aic = -2*logL + 2*d,tau=round(tau,4),pi=pi,bic=-2*logL + log(S)*d,ICL= -2*logL + log(S)*d +2*EN,coef=fm.out,covar=0,aic.full= -2*logL.full + 2*d,bic.full= -2*logL.full + log(S)*d))
}
"fitmix.cpp" <-
function (form,datsp,sp,G=2,pars=NA,trace=TRUE,calc.hes=FALSE,r1)
{
if(!is.numeric(sp)){
sp <- as.integer(factor(sp))
}
sp.name <- unique(sp)
S <- length(unique(sp))
n <- length(which(sp==sp.name[1]))
X <- model.matrix(form, data = datsp[sp==sp.name[1],])
##X <- model.frame(form, data = datsp)
# X <- X[sp==sp.name[1],]
#y <- model.response(form, data = datsp)
y <- datsp$obs
if(is.na(pars[1])) {
##pars <- rep(0.01,G-1+(ncol(X)*G))
pars <- runif(G-1+(ncol(X)*G),-2,2)
pars[1:(G-1)] <- runif(G-1)
}
offset<-rep(0,n)
## hes <- rep(0,length(pars)^2)
gradient <- rep(0,length(pars))
tau <- matrix(0,S,G) ##must leave this in as defines S & G
model.type<- as.integer(1)
loglike <- try(.Call("SpeciesMix",pars,y,X,sp,tau,gradient,offset,model.type,PACKAGE="SpeciesMix"))
calc.deriv <- function(p){
gradient <- rep(0,length(pars))
ll <- .Call("Calculate_Gradient",p,y,X,sp,tau,gradient,offset,as.integer(1),PACKAGE="SpeciesMix")
return(gradient)
}
hes <- 0
if(calc.hes){
hes <- nd2(pars,calc.deriv)
dim(hes) <- rep(length(pars),2)
dim(hes) <- rep(length(pars),2)
rownames(hes) <- colnames(hes) <- c(paste("G.",1:(G-1),sep=""),paste("G",1:G,rep(colnames(X),each=G),sep="."))
}
if(!is.numeric(loglike)) loglike <- 0
pi <- pars[1:(G-1)]
coef <- pars[ (G):length(pars)]
r.logl <- logLmix(pars,list(y=y,x=model.matrix(form, data = datsp)),G,S,sp,sp.name,out.tau=TRUE)
pi <- additive.logistic(pi)
names(pi) <- paste("G.",1:G,sep="")
coef <- matrix(coef,G,ncol(X))
rownames(coef) <- paste("G.",1:G,sep="")
colnames(coef) <- colnames(X)
AIC <- 2*loglike + 2*length(pars)
BIC <- 2*loglike + log(S)*length(pars)
##list(logl=loglike,r.logl=r.logl$logl,pi=pi,coef=coef,tau=round(exp(r.logl$tau),4),aic=AIC,bic=BIC,hessian=hes,gradient=gradient)
list(logl=loglike,pi=pi,coef=coef,tau=round(exp(r.logl$tau),4),aic=AIC,bic=BIC,hessian=hes,gradient=gradient)
}
"fitmix.gaussian" <-
function (form,datsp,sp,G=2,ite.max=500,trace=TRUE,full.model=FALSE)
{
## fitting abundance data with a negative binomial distribution
## dat2 has colums obs,sp
##
temp.warn <- getOption( "warn")
options( warn=-1)
sp.name <- unique(sp)
S <- length(unique(sp))
n <- length(which(sp==sp.name[1]))
cat("Fitting Group",G,"\n")
if(trace) cat("Iteration | LogL \n")
##dat.tau <- data.frame(matrix(0,dim(datsp)[1],G))
dat.tau <- 0
##dat <- data.frame(datsp,dat.tau)
pi <- rep(0,G)
ite <- 1
logL <- -99999999
old.logL <- -88888888
## set up initial GLM
t1 <- create.starting.values.gaussian(S,G,n,form,datsp)
pi <- t1$pi;fmM <- t1$fmM;tau <- t1$tau;dat.tau <- t1$dat.tau;first.fit <- t1$first.fit
while(abs(logL-old.logL) > 0.0001 & ite<=ite.max){
old.logL <- logL
for(i in 1:G){
##dat.tau[,i] <- rep(tau[,i],each=n)
pi[i] <- sum(tau[,i])/S
}
if(any(pi==0)) { ## occasionally with complicated models the random starting values result in a pi[i]==0; so restart with new random starts
cat("pi has gone to zero - restarting fitting \n")
t1 <- create.starting.values.gaussian(S,G,n,form,datsp)
pi <- t1$pi;fmM <- t1$fmM;tau <- t1$tau;dat.tau <- t1$dat.tau;first.fit <- t1$first.fit
ite <- 1
}
fmM <- lapply(1:G,weighted.glm.gaussian,first.fit,tau,n,fmM,sp)
logL <- 0
tmp.like <- matrix(0,S,G)
est.tau <- lapply(1:S,estimate.pi.gaussian,sp,sp.name,datsp,fmM,pi,G,first.fit)
for(j in 1:S){
if(is.atomic(est.tau[[j]])){ print (est.tau[[j]])} else
{
tau[j,] <- est.tau[[j]]$tau
logL <- logL+est.tau[[j]]$sum.like
}
}
if(trace) cat(ite," | ",logL,"\n")
ite <- ite+1
}
fm.out <- data.frame(matrix(0,G,length(fmM[[1]]$coef)))
int.out <- data.frame(matrix(0,S,G))
names(int.out) <- paste("grp.",1:G,sep="")
names(fm.out) <- names(fmM[[1]]$coef)
tau <- data.frame(tau)
names(tau) <- paste("grp.",1:G,sep="")
EN <- -sum(unlist(tau)*log(unlist(tau)))
d <- length(unlist(fm.out)) + length(tau)-1
fm.theta <- rep(0,G)
for(i in 1:G) {
fm.out[i,] <- fmM[[i]]$coef
int.out[,i] <- fmM[[i]]$sp.intercept
##dat.tau[,i] <- rep(tau[,i],each=n)
fm.theta[i] <- fmM[[i]]$theta
}
names(pi) <- paste("G",1:G,sep=".")
t.pi <- additive.logistic(pi,TRUE)
parms <- c(t.pi[1:(G-1)],fm.theta,unlist(fm.out),unlist(int.out))
logL.full <- logL
logL <- logLmix.gaussian(parms,first.fit,G,S,sp,sp.name)
options(warn=temp.warn)
if(full.model) return(list(logl=logL,aic = -2*logL + 2*d,tau=round(tau,4),pi=pi,bic=-2*logL + log(S)*d,ICL= -2*logL + log(S)*d +2*EN,coef=fm.out,sp.intercept=int.out,theta=fm.theta,fmM=fmM,model.tau=dat.tau,covar=0,aic.full= -2*logL.full + 2*d,bic.full= -2*logL.full + log(S)*d))
return(list(logl=logL,aic = -2*logL + 2*d,tau=round(tau,4),pi=pi,bic=-2*logL + log(S)*d,ICL= -2*logL + log(S)*d +2*EN,coef=fm.out,sp.intercept=int.out,theta=fm.theta,covar=0,aic.full= -2*logL.full + 2*d,bic.full= -2*logL.full + log(S)*d))
}
"fitmix.nbinom" <-
function (form, datsp, sp, G = 2, ite.max = 500, trace = TRUE,
full.model = FALSE)
{
temp.warn <- getOption("warn")
options(warn = -1)
sp.name <- unique(sp)
S <- length(unique(sp))
n <- length(which(sp == sp.name[1]))
cat("Fitting Group", G, "\n")
if (trace)
cat("Iteration | LogL \n")
dat.tau <- 0
pi <- rep(0, G)
ite <- 1
logL <- -99999999
old.logL <- -88888888
if (ite != 1)
t1 <- create.starting.values.nbinom(S, G, n, form, datsp)
t1 <- create.starting.values.nbinom.kmeans(S, G, n, form,
datsp)
pi <- t1$pi
fmM <- t1$fmM
tau <- t1$tau
dat.tau <- t1$dat.tau
first.fit <- t1$first.fit
while (abs(logL - old.logL) > 1e-04 & ite <= ite.max) {
old.logL <- logL
for (i in 1:G) {
pi[i] <- sum(tau[, i])/S
}
if (any(pi == 0)) {
cat("pi has gone to zero - restarting fitting \n")
t1 <- create.starting.values.nbinom(S, G, n, form,
datsp)
pi <- t1$pi
fmM <- t1$fmM
tau <- t1$tau
dat.tau <- t1$dat.tau
first.fit <- t1$first.fit
ite <- 1
}
fmM <- lapply(1:G, weighted.glm.nbinom, first.fit, tau,
n, fmM, sp)
for (j in 1:S) {
tmp <- rep(0, G)
for (g in 1:G) tmp[g] <- fmM[[g]]$sp.intercept[j]
tmp <- sum(tmp * tau[j, ])
for (g in 1:G) fmM[[g]]$sp.intercept[j] <- tmp
}
logL <- 0
tmp.like <- matrix(0, S, G)
est.tau <- lapply(1:S, estimate.pi.nbinom, sp, sp.name,
datsp, fmM, pi, G, first.fit)
for (j in 1:S) {
if (is.atomic(est.tau[[j]])) {
print(est.tau[[j]])
}
else {
tau[j, ] <- est.tau[[j]]$tau
logL <- logL + est.tau[[j]]$sum.like
}
}
if (trace)
cat(ite, " | ", logL, "\n")
ite <- ite + 1
}
fm.out <- data.frame(matrix(0, G, length(fmM[[1]]$coef)))
int.out <- rep(0, S)
names(fm.out) <- names(fmM[[1]]$coef)
tau <- data.frame(tau)
names(tau) <- paste("grp.", 1:G, sep = "")
EN <- -sum(unlist(tau) * log(unlist(tau)))
d <- length(unlist(fm.out)) + length(tau) - 1
fm.theta <- rep(0, G)
int.out <- fmM[[1]]$sp.intercept
for (i in 1:G) {
fm.out[i, ] <- fmM[[i]]$coef
}
names(pi) <- paste("G", 1:G, sep = ".")
t.pi <- additive.logistic(pi, TRUE)
parms <- c(t.pi[1:(G - 1)], unlist(fm.out), int.out, rep(1,
S))
logL.full <- logL
logL <- logLmix.nbinom(parms, first.fit, G, S, sp, sp.name)
options(warn = temp.warn)
if (full.model)
return(list(logl = logL, aic = -2 * logL + 2 * d, tau = round(tau,
4), pi = pi, bic = -2 * logL + log(S) * d, ICL = -2 *
logL + log(S) * d + 2 * EN, coef = fm.out, sp.intercept = int.out,
theta = fm.theta, fmM = fmM, model.tau = dat.tau,
covar = 0, aic.full = -2 * logL.full + 2 * d, bic.full = -2 *
logL.full + log(S) * d, pars = parms))
return(list(logl = logL, aic = -2 * logL + 2 * d, tau = round(tau,
4), pi = pi, bic = -2 * logL + log(S) * d, ICL = -2 *
logL + log(S) * d + 2 * EN, coef = fm.out, sp.intercept = int.out,
theta = fm.theta, covar = 0, aic.full = -2 * logL.full +
2 * d, bic.full = -2 * logL.full + log(S) * d, pars = parms))
}
"fitmix.nbinom.cpp" <-
function (form,datsp,sp,G=2,pars=NA,trace=TRUE,calc.hes=FALSE)
{
if(!is.numeric(sp)){
sp <- as.integer(factor(sp))
}
sp.name <- unique(sp)
S <- length(unique(sp))
n <- length(which(sp==sp.name[1]))
X <- model.matrix(form, data = datsp[sp==sp.name[1],])
## offset <- model.offset(form, data = datsp[sp==sp.name[1],])
offset <- model.frame(form, data = datsp[sp==sp.name[1],])
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,n)
##offset <- rep(0,n)
##X <- model.frame(form, data = datsp)
# X <- X[sp==sp.name[1],]
#y <- model.response(form, data = datsp)
y <- datsp$obs
if(is.na(pars[1])) {
##pars <- rep(0.01,G-1+(ncol(X)*G))
sp.int <- rep(0.5,S)
sp.dispersion <- rep(1,S)
fm <- matrix(runif(ncol(X)*G,-1,1),G,ncol(X))
pars <- c(runif(G-1),unlist(fm),sp.int,sp.dispersion)
}
## hes <- rep(0,length(pars)^2)
gradient <- rep(0,length(pars))
tau <- matrix(0,S,G) ##must leave this in as defines S & G
loglike <- try(.Call("SpeciesMix",pars,y,X,sp,tau,gradient,offset,as.integer(2),PACKAGE="SpeciesMix"))
calc.deriv <- function(p){
gradient <- rep(0,length(pars))
ll <- .Call("Calculate_Gradient",p,y,X,sp,tau,gradient,offset,as.integer(2),PACKAGE="SpeciesMix")
return(gradient)
}
r.deriv <- function(p){ logLmix.nbinom(p,list(y=y,x=model.matrix(form, data = datsp)),G,S,sp,sp.name,out.tau=FALSE)}
#r.grad <- nd2(pars,r.deriv)
##print(r.grad)
hes <- 0
covar <- 0
if(calc.hes){
hes <- nd2(pars,calc.deriv)
dim(hes) <- rep(length(pars),2)
dim(hes) <- rep(length(pars),2)
covar <- try(solve(hes))
#rownames(hes) <- colnames(hes) <- c(paste("G.",1:(G-1),sep=""),paste("G",1:G,rep(colnames(X),each=G),sep="."))
}
if(!is.numeric(loglike)) loglike <- 0
pi <- pars[1:(G-1)]
#coef <- pars[ (G):length(pars)]
coef <- pars[-1*(1:((G-1)))] ## remove pi
sp.int <- coef[(length(coef)-(2*S-1)):(length(coef)-S)]
sp.dispersion <- coef[(length(coef)-(S-1)):length(coef)]
fm <- coef[-1*((length(coef)-(2*S-1)):length(coef))]
offset <- model.frame(form, data = datsp)
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,nrow(datsp))
r.logl <- logLmix.nbinom(pars,list(y=y,x=model.matrix(form, data = datsp),offset=offset),G,S,sp,sp.name,out.tau=TRUE)
print(r.logl$logl)
pi <- additive.logistic(pi)
names(pi) <- paste("G.",1:G,sep="")
coef <- matrix(fm,G,ncol(X))
rownames(coef) <- paste("G.",1:G,sep="")
colnames(coef) <- colnames(X)
AIC <- 2*loglike + 2*length(pars)
BIC <- 2*loglike + log(S)*length(pars)
##list(logl=loglike,r.logl=r.logl$logl,pi=pi,coef=coef,tau=round(exp(r.logl$tau),4),aic=AIC,bic=BIC,hessian=hes,gradient=gradient)
list(logl=loglike,pi=pi,coef=coef,sp.intercept=sp.int,sp.dispersion=sp.dispersion,tau=round(exp(r.logl$tau),4),aic=AIC,bic=BIC,hessian=hes,gradient=gradient,covar=covar)#,r.grad=r.grad)
}
"fitmix.tweedie" <-
function (form, datsp, sp, G = 2, ite.max = 500, trace = TRUE,
full.model = FALSE)
{
temp.warn <- getOption("warn")
options(warn = -1)
sp.name <- unique(sp)
S <- length(unique(sp))
n <- length(which(sp == sp.name[1]))
cat("Fitting Group", G, "\n")
if (trace)
cat("Iteration | LogL \n")
dat.tau <- 0
pi <- rep(0, G)
ite <- 1
logL <- -99999999
old.logL <- -88888888
if (ite != 1)
t1 <- create.starting.values.tweedie(S, G, n, form, datsp)
t1 <- create.starting.values.tweedie.kmeans(S, G, n, form,
datsp)
pi <- t1$pi
fmM <- t1$fmM
tau <- t1$tau
dat.tau <- t1$dat.tau
first.fit <- t1$first.fit
while (abs(logL - old.logL) > 1e-04 & ite <= ite.max) {
old.logL <- logL
for (i in 1:G) {
pi[i] <- sum(tau[, i])/S
}
if (any(pi == 0)) {
cat("pi has gone to zero - restarting fitting \n")
t1 <- create.starting.values.tweedie(S, G, n, form,
datsp)
pi <- t1$pi
fmM <- t1$fmM
tau <- t1$tau
dat.tau <- t1$dat.tau
first.fit <- t1$first.fit
ite <- 1
}
fmM <- lapply(1:G, weighted.glm.tweedie, first.fit, tau,
n, fmM, sp)
for (j in 1:S) {
tmp <- rep(0, G)
for (g in 1:G) tmp[g] <- fmM[[g]]$sp.intercept[j]
tmp <- sum(tmp * tau[j, ])
for (g in 1:G) fmM[[g]]$sp.intercept[j] <- tmp
}
logL <- 0
tmp.like <- matrix(0, S, G)
est.tau <- lapply(1:S, estimate.pi.tweedie, sp, sp.name,
datsp, fmM, pi, G, first.fit)
for (j in 1:S) {
if (is.atomic(est.tau[[j]])) {
print(est.tau[[j]])
}
else {
tau[j, ] <- est.tau[[j]]$tau
logL <- logL + est.tau[[j]]$sum.like
}
}
if (trace)
cat(ite, " | ", logL, "\n")
ite <- ite + 1
}
fm.out <- data.frame(matrix(0, G, length(fmM[[1]]$coef)))
int.out <- fmM[[1]]$sp.intercept
names(fm.out) <- names(fmM[[1]]$coef)
tau <- data.frame(tau)
names(tau) <- paste("grp.", 1:G, sep = "")
EN <- -sum(unlist(tau) * log(unlist(tau)))
d <- length(unlist(fm.out)) + length(tau) - 1
fm.phi <- fm.p <- rep(0, G)
for (i in 1:G) {
fm.out[i, ] <- fmM[[i]]$coef
fm.phi[i] <- fmM[[i]]$phi
fm.p[i] <- fmM[[i]]$p
}
names(pi) <- paste("G", 1:G, sep = ".")
t.pi <- additive.logistic(pi, TRUE)
parms <- c(t.pi[1:(G - 1)], unlist(fm.out), int.out, rep(2,
S))
names(parms) <- c(rep("pi", G - 1), rep("coef", length(unlist(fm.out))),
rep("int", length(int.out)), rep("phi", S))
logL.full <- logL
logL <- logLmix.tweedie(parms, first.fit, G, S, sp, sp.name)
print(logL)
options(warn = temp.warn)
if (full.model)
return(list(logl = logL, aic = -2 * logL + 2 * d, tau = round(tau,
4), pi = pi, bic = -2 * logL + log(S) * d, ICL = -2 *
logL + log(S) * d + 2 * EN, coef = fm.out, sp.intercept = int.out,
phi = fm.phi, p = fm.p, fmM = fmM, model.tau = dat.tau,
covar = 0, aic.full = -2 * logL.full + 2 * d, bic.full = -2 *
logL.full + log(S) * d, pars = parms))
return(list(logl = logL, aic = -2 * logL + 2 * d, tau = round(tau,
4), pi = pi, bic = -2 * logL + log(S) * d, ICL = -2 *
logL + log(S) * d + 2 * EN, coef = fm.out, sp.intercept = int.out,
phi = fm.phi, p = fm.p, covar = 0, aic.full = -2 * logL.full +
2 * d, bic.full = -2 * logL.full + log(S) * d, pars = parms))
}
"fitmix.tweedie.cpp" <-
function (form, datsp, sp, G = 2, pars = NA, trace = TRUE, calc.hes = FALSE)
{
if (!is.numeric(sp)) {
sp <- as.integer(factor(sp))
}
sp.name <- unique(sp)
S <- length(unique(sp))
n <- length(which(sp == sp.name[1]))
X <- model.matrix(form, data = datsp[sp == sp.name[1], ])
offset <- model.frame(form, data = datsp[sp==sp.name[1],])
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,n)
y <- datsp$obs
if (is.na(pars[1])) {
sp.int <- rep(0.5, S)
sp.phi <- rep(1, S)
sp.p <- rep(1.6, S)
fm <- matrix(rep(0.5, ncol(X) * G), G, ncol(X))
pars <- c(rep(0.5, G - 1), unlist(fm), sp.int, sp.phi)
}
pars.og <- pars
pars.og[1] <- pars.og[1] * 0.99
gradient <- rep(0, length(pars))
tau <- matrix(0, S, G)
loglike <- try(.Call("SpeciesMix", pars, y, X, sp, tau, gradient,
offset, as.integer(3),PACKAGE="SpeciesMix"))
calc.deriv <- function(p) {
gradient <- rep(0, length(pars))
ll <- .Call("Calculate_Gradient", p, y, X, sp, tau, gradient,
offset, as.integer(3),PACKAGE="SpeciesMix")
return(gradient)
}
hes <- 0
if (calc.hes) {
hes <- jacobian(calc.deriv, pars, method = "simple")
dim(hes) <- rep(length(pars), 2)
dim(hes) <- rep(length(pars), 2)
}
if (!is.numeric(loglike))
loglike <- 0
pi <- pars[1:(G - 1)]
coef <- pars[-1 * (1:((G - 1)))]
sp.int <- coef[(length(coef) - (2 * S - 1)):(length(coef) -
S)]
sp.dispersion <- coef[(length(coef) - (S - 1)):length(coef)]
sp.phi <- sp.dispersion[1:S]
sp.p <- rep(1.6, S)
fm <- coef[-1 * ((length(coef) - (2 * S - 1)):length(coef))]
offset <- model.frame(form, data = datsp)
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,nrow(datsp))
names(gradient) <- names(pars) <- c(rep("pi", G - 1), rep("coef",
length(unlist(fm))), rep("int", length(sp.int)), rep("phi",
S))
r.deriv <- function(p) {
logLmix.tweedie(p, list(y = y, x = model.matrix(form,
data = datsp)), G, S, sp, sp.name, out.tau = FALSE)
}
r.logl <- logLmix.tweedie(pars, list(y = y, x = model.matrix(form,
data = datsp),offset=offset), G, S, sp, sp.name, out.tau = TRUE)
print(r.logl$logl)
pi <- additive.logistic(pi)
names(pi) <- paste("G.", 1:G, sep = "")
coef <- matrix(fm, G, ncol(X))
rownames(coef) <- paste("G.", 1:G, sep = "")
colnames(coef) <- colnames(X)
AIC <- 2 * loglike + 2 * length(pars)
BIC <- 2 * loglike + log(S) * length(pars)
list(logl = loglike, pi = pi, coef = coef, sp.intercept = sp.int,
phi = sp.phi, p = sp.p, tau = round(exp(r.logl$tau),
4), aic = AIC, bic = BIC, hessian = hes, gradient = gradient)
}
"glm.fit.nbinom" <-
function (x,y,offset=NULL,weights=NULL,mustart=NULL,est.var=FALSE)
{
X <- x
if(is.null(offset)) offset <- 0
if(is.null(weights)) weights <- rep(1,length(y))
gradient <- rep(0,ncol(X)+1)
if(is.null(mustart)){ pars <- gradient+1} else{pars <- mustart}
fitted.values <- rep(0,length(y))
logl <- .Call("Neg_Bin",pars,X,y,weights,offset,gradient,fitted.values,PACKAGE="SpeciesMix")
vcov <- 0
se <- rep(0,length(pars))
if(est.var) {
calc.deriv <- function(p){
gradient <- rep(0,length(pars))
ll <- .Call("Neg_Bin_Gradient",p,X,y,weights,offset,gradient,PACKAGE="SpeciesMix")
return(gradient)
}
hes <- nd2(pars,calc.deriv)
dim(hes) <- rep(length(pars),2)
vcov <- try(solve(hes))
se <- try(sqrt(diag(vcov)))
colnames(vcov) <- rownames(vcov) <- c("theta",colnames(X))
}
names(pars) <- names(se) <- names(gradient) <- c("theta",colnames(X))
return(list(logl=logl,coef=pars[-1],theta=pars[1],se=se[-1],se.theta=se[1],fitted=fitted.values,gradient=gradient,vcov=vcov))
}
"glm.nbinom" <-
function (form,data,weights=NULL,mustart=NULL,est.var=FALSE)
{
X <- model.matrix(form,data)
t1 <- model.frame(form,data)
y <- model.response(t1)
offset <- model.offset(t1)
if(is.null(offset)) offset <- 0
if(is.null(weights)) weights <- rep(1,length(y))
gradient <- rep(0,ncol(X)+1)
if(is.null(mustart)){ pars <- gradient+1}else{pars <- mustart}
fitted.values <- rep(0,length(y))
logl <- .Call("Neg_Bin",pars,X,y,weights,offset,gradient,fitted.values,PACKAGE="SpeciesMix")
vcov <- 0
se <- rep(0,length(pars))
if(est.var) {
calc.deriv <- function(p){
gradient <- rep(0,length(pars))
ll <- .Call("Neg_Bin_Gradient",p,X,y,weights,offset,gradient,PACKAGE="SpeciesMix")
return(gradient)
}
hes <- nd2(pars,calc.deriv)
dim(hes) <- rep(length(pars),2)
vcov <- try(solve(hes))
se <- try(sqrt(diag(vcov)))
colnames(vcov) <- rownames(vcov) <- c("theta",colnames(X))
}
names(pars) <- names(se) <- names(gradient) <- c("theta",colnames(X))
return(list(logl=logl,coef=pars[-1],theta=pars[1],se=se[-1],se.theta=se[1],fitted=fitted.values,gradient=gradient,vcov=vcov))
}
"ldTweedie.lp" <-
function ( parms, y, X.p, offsetty, phi, p, wts=rep( 1, length( y)))
{
mu <- exp( X.p %*% parms[1:ncol( X.p)] + offsetty)
if( is.null( phi) & is.null( p)){
phi <- parms[ncol( X.p) + 1]
p <- parms[ncol( X.p) + 2]
}
if( is.null( phi) & !is.null( p))
phi <- parms[ncol( X.p)+1]
if( !is.null( phi) & is.null( p))
p <- parms[ncol( X.p)+1]
lambda <- ( mu^( 2-p)) / ( phi*(2-p))
alpha <- ( 2-p) / ( p-1)
tau <- phi*(p-1)*mu^(p-1)
mu.Z <- alpha * tau
return( -sum( wts * dPoisGam( y, lambda=lambda, mu.Z=mu.Z, alpha=alpha, LOG=TRUE)))
}
"ldTweedie.lp.deriv" <-
function ( parms, y, X.p, offsetty, phi, p, wts=rep( 1, length( y)))
{
mu <- exp( X.p %*% parms[1:ncol( X.p)] + offsetty)
p.flag <- phi.flag <- FALSE
if( is.null( phi) & is.null( p)){
p.flag <- phi.flag <- TRUE
phi <- parms[ncol( X.p) + 1]
p <- parms[ncol( X.p) + 2]
}
if( is.null( phi) & !is.null( p)){
phi <- parms[ncol( X.p)+1]
phi.flag <- TRUE
}
if( !is.null( phi) & is.null( p)){
p <- parms[ncol( X.p)+1]
p.flag <- TRUE
}
lambda <- ( mu^( 2-p)) / ( phi*(2-p))
alpha <- ( 2-p) / ( p-1)
tau <- phi*(p-1)*mu^(p-1)
mu.Z <- alpha * tau
dTweedparms <- -wts * dPoisGamDerivs( y, lambda=lambda, mu.Z=mu.Z, alpha=alpha)
DTweedparmsDmu <- matrix( c( ( mu^(1-p)) / phi, alpha*phi*( ( p-1)^2)*( mu^(p-2)), rep( 0, length( mu))), nrow=3, byrow=T)
tmp <- rowSums( dTweedparms * t( DTweedparmsDmu))
tmp <- tmp * mu
tmp <- apply( X.p, 2, function( x) x*tmp)
derivs <- colSums( tmp)
if( phi.flag){
DTweedparmsDphi <- matrix( c( -( ( mu^(2-p)) / ( ( phi^2)*(2-p))), alpha*( p-1)*( mu^( p-1)), rep( 0, length( mu))), nrow=3, byrow=T)
tmpPhi <- rowSums( dTweedparms * t( DTweedparmsDphi))#vectorised way of doing odd calculation
derivs <- c( derivs, sum( tmpPhi))
names( derivs)[length( derivs)] <- "phi"
}
if( p.flag){
dalphadp <- -( 1+alpha) / ( p-1)
DTweedparmsDp <- matrix( c( lambda*( 1/(2-p) - log( mu)), mu.Z*( dalphadp/alpha + 1/( p-1) + log( mu)), rep( dalphadp, length( y))), nrow=3, byrow=T)
tmpP <- rowSums( dTweedparms * t( DTweedparmsDp))
derivs <- c( derivs, sum( tmpP))
names( derivs)[length( derivs)] <- "p"
}
return( derivs)
}
"Lmix" <-
function (pars,y,x,G)
{
fm <- pars[-1*(1:(G-1))]
pi <- pars[(1:(G-1))]
dim(fm) <- c(G,(length(pars)-(G-1))/G)
pi <- additive.logistic(pi)
S <- dim(y)[2]
link.fun <- make.link("logit")
like <- 1
lf <- matrix(0,G,S)
lg <- rep(0,S)
dBi <- array(0,dim=c(S,G,dim(fm)[2]))
for(s in 1:S){
tmp.like <- rep(0,G)
for(g in 1:G){
lpre <- x%*%fm[g,]
for(j in 1:dim(fm)[2]){
dBi[s,g,j] <- sum((y[,s]-link.fun$linkinv(lpre))*x[,j])
}
tmp.like[g] <- pi[g]*prod(dbinom(y[,s],1,link.fun$linkinv(lpre),log=F))
lf[g,s] <- prod(dbinom(y[,s],1,link.fun$linkinv(lpre),log=F))
}
lg[s] <- sum(tmp.like)
like <- like*sum(tmp.like)
}
#print(dBi)
# print(log(lg))
# print(log(lf))
dl.dpi <- rep(0,G)
for(g in 1:G) dl.dpi[g] <- ( sum( exp(-log(lg) + log(lf[g,]))))
der <- matrix(0,dim(fm)[1],dim(fm)[2])
for(j in 1:dim(fm)[2]){
for(g in 1:G){
for(s in 1:S){
der[g,j] <- der[g,j]+ exp( -log(lg[s]) + log(pi[g]) + log(lf[g,s])) * dBi[s,g,j]
#if(g==1 & j == 1) print(c(der[g,j],-log(lg[s]),pi[g],log(lf[g,s]),dBi[s,g,j]))
}
}
}
dpi.deta <- matrix(0,G-1,G)
ad.trans <- 1+sum(exp(pars[1:(G-1)]))
for(i in 1:(G-1))
for(g in 1:G-1){
if(i==g) {dpi.deta[i,g] <- exp(pars[i])*exp(pars[g])/ad.trans^2}
else{ dpi.deta[i,g] <- exp(pars[i])/ad.trans - exp(2*pars[g])/ad.trans^2}
}
dpi.deta[,G] <- -rowSums(dpi.deta)
print(dpi.deta)
print(dl.dpi)
d1 <- dpi.deta%*%dl.dpi
print(d1)
list(like,0-der)
}
"logLike.pars" <-
function (pi,coef,sp.form,sp.data,covar.data)
{
G <- length(pi)
pars <- c(additive.logistic(pi,T)[1:(G-1)],unlist(coef))
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=unlist(sp.data),covar.data)
names(data)[1] <- as.character(sp.form)[2]
first.fit <- list(y=data[,1],x=model.matrix(sp.form,data=data))
logl <- -logLmix(pars,first.fit,G,S,sp,sp.name)
logl
}
"logLmix" <-
function (pars,first.fit,G,S,sp,spname,out.tau=FALSE)
{
tau <- matrix(0,S,G)
##tau,out.tau=FALSE
if(G>1){
fm <- pars[-1*(1:(G-1))]
pi <- pars[(1:(G-1))]
## fm <- tau[-1*((length(tau)-(G-2)):length(tau))]
#pi <- tau[((length(tau)-(G-2)):length(tau))]
dim(fm) <- c(G,(length(pars)-(G-1))/G)
##pi[G] <- 1-sum(pi)
pi <- additive.logistic(pi)
} else{
fm <- tau[1:(length(pars)-1)]
dim(fm) <- c(1,length(fm))
pi <- 1
}
link.fun <- make.link("logit")
log.like <- 0
for(j in 1:S){
sel.sp <- which(sp==spname[j])
tmp.like <- rep(0,G)
for(i in 1:G){
if(length(fm[i,])==1){lpre <- first.fit$x[sel.sp,]*fm[i,]
}else{ lpre <- first.fit$x[sel.sp,]%*%fm[i,]}
tmp.like[i] <- sum(dbinom(first.fit$y[sel.sp],1,link.fun$linkinv(lpre),log=T))
}
eps <- max(tmp.like)
log.like <- log.like + (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
tau[j,] <- log(pi) + tmp.like - (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
}
if(out.tau)return(list(logl=log.like,tau=tau))
log.like
}
"logLmix.gaussian" <-
function (pars,first.fit,G,S,sp,spname,out.tau=FALSE)
{
tau <- matrix(0,S,G)
##tau,out.tau=FALSE
if(G>1){
fm <- pars[-1*(1:((G-1)+G))] ## remove pi
sp.int <- fm[(length(fm)-(S*G-1)):length(fm)]
fm <- fm[-1*((length(fm)-(S*G-1)):length(fm))]
pi <- pars[(1:(G-1))]
theta <- pars[G:(G+G-1)]
## fm <- tau[-1*((length(tau)-(G-2)):length(tau))]
#pi <- tau[((length(tau)-(G-2)):length(tau))]
dim(sp.int) <- c(S,G)
dim(fm) <- c(G,length(fm)/G)
##pi[G] <- 1-sum(pi)
pi <- additive.logistic(pi)
} else{
return(0)
fm <- tau[1:(length(pars)-1)]
dim(fm) <- c(1,length(fm))
pi <- 1
}
log.like <- 0
for(j in 1:S){
sel.sp <- which(sp==spname[j])
tmp.like <- rep(0,G)
for(i in 1:G){
lpre <- cbind(1,first.fit$x[sel.sp,])%*%c(sp.int[j,i],fm[i,])
##tmp.like[i] <- sum(dpois(first.fit$y[sel.sp],exp(lpre),log=T))
tmp.like[i] <- sum(dnorm(first.fit$y[sel.sp],mean=lpre,sd=theta[i],log=T))
eps <- max(tmp.like)
log.like <- log.like + (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
tau[j,] <- log(pi) + tmp.like - (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
}
}
if(out.tau)return(list(logl=log.like,tau=tau))
log.like
}
"logLmix.nbinom" <-
function (pars,first.fit,G,S,sp,spname,out.tau=FALSE)
{
tau <- matrix(0,S,G)
##tau,out.tau=FALSE
if(G>1){
fm <- pars[-1*(1:((G-1)))] ## remove pi
##sp.int <- fm[(length(fm)-(S*G-1)):length(fm)]
sp.int <- fm[(length(fm)-(2*S-1)):(length(fm)-S)]
sp.dispersion <- fm[(length(fm)-(S-1)):length(fm)]
##fm <- fm[-1*((length(fm)-(S*G-1)):length(fm))]
fm <- fm[-1*((length(fm)-(2*S-1)):length(fm))]
pi <- pars[(1:(G-1))]
theta <- pars[G:(G+G-1)]
## fm <- tau[-1*((length(tau)-(G-2)):length(tau))]
#pi <- tau[((length(tau)-(G-2)):length(tau))]
##dim(sp.int) <- c(S,G)
dim(fm) <- c(G,length(fm)/G)
##pi[G] <- 1-sum(pi)
pi <- additive.logistic(pi)
} else{
return(0)
fm <- tau[1:(length(pars)-1)]
dim(fm) <- c(1,length(fm))
pi <- 1
}
log.like <- 0
for(j in 1:S){
sel.sp <- which(sp==spname[j])
tmp.like <- rep(0,G)
for(i in 1:G){
lpre <- cbind(1,first.fit$x[sel.sp,])%*%c(sp.int[j],fm[i,])+first.fit$offset[sel.sp]##,i],fm[i,])
##tmp.like[i] <- sum(dpois(first.fit$y[sel.sp],exp(lpre),log=T))
tmp.like[i] <- sum(dnbinom(first.fit$y[sel.sp],mu=exp(lpre),size=sp.dispersion[j],log=T))
}
# print(tmp.like)
eps <- max(tmp.like)
log.like <- log.like + (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
tau[j,] <- log(pi) + tmp.like - (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
}
if(out.tau)return(list(logl=log.like,tau=tau))
log.like
}
"logLmix.tweedie" <-
function (pars,first.fit,G,S,sp,spname,out.tau=FALSE)
{
tau <- matrix(0,S,G)
##tau,out.tau=FALSE
if(G>1){
pi <- pars[which(names(pars)=="pi")]
phi <- pars[which(names(pars)=="phi")]
p <- rep(1.6,S)
fm <- pars[which(names(pars)=="coef")]
sp.int <- pars[which(names(pars)=="int")]
##fm <- pars[-1*(1:((G-1)+G))] ## remove pi
##sp.int <- fm[(length(fm)-(S*G-1)):length(fm)]
##fm <- fm[-1*((length(fm)-(S*G-1)):length(fm))]
##pi <- pars[(1:(G-1))]
##theta <- pars[G:(G+G-1)]
## fm <- tau[-1*((length(tau)-(G-2)):length(tau))]
#pi <- tau[((length(tau)-(G-2)):length(tau))]
##dim(sp.int) <- c(S,G)
dim(fm) <- c(G,length(fm)/G)
##pi[G] <- 1-sum(pi)
pi <- additive.logistic(pi)
} else{
return(0)
fm <- tau[1:(length(pars)-1)]
dim(fm) <- c(1,length(fm))
pi <- 1
}
log.like <- 0
for(j in 1:S){
sel.sp <- which(sp==spname[j])
tmp.like <- rep(0,G)
for(i in 1:G){
lpre <- cbind(1,first.fit$x[sel.sp,])%*%c(sp.int[j],fm[i,])+first.fit$offset[sel.sp]
##tmp.like[i] <- sum(dpois(first.fit$y[sel.sp],exp(lpre),log=T))
tmp.like[i] <- sum(dTweedie(y=first.fit$y[sel.sp],mu=exp(lpre),phi=phi[j],p=p[j],LOG=T))
}
##print(tmp.like)
eps <- max(tmp.like)
log.like <- log.like + (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
tau[j,] <- log(pi) + tmp.like - (log(sum(pi*exp((tmp.like)-(eps))))+(eps))
}
if(out.tau)return(list(logl=log.like,tau=tau))
log.like
}
"mix.residuals" <-
function (fmM,form,datsp,sp)
{
cat("calculating residuals \n")
link.fun <- make.link("logit")
x <- model.matrix(form,data=datsp)
spname <- unique(sp)
S <- length(spname)
G <- ncol(fmM$tau)
PIT <- matrix(NA,S,G)
for(g in 1:G){
for(s in 1:length(spname)){
sel.sp <- which(sp==spname[s])
t.obs <- sum(datsp$obs[sel.sp])
pre <- link.fun$linkinv(x[sel.sp,]%*%fmM$coef[g,])
obs <- datsp$obs[sel.sp]
dis <- distr.binom(pre)
# PIT[s,g] <- qnorm(sum(dis[1:(length(dis)-1)],dis[length(dis)]/2))
nSucc <- sum( obs)
transfo <- sum( dis[1:nSucc],dis[nSucc+1]/2)
transfo <- min( transfo, 1)
##transfor <- max( transfo, 0)
PIT[s,g] <- qnorm( transfo)
}
}
PIT
}
"nd2" <-
function( x0, f, m=NULL, D.accur=4, ...) {
# A function to compute highly accurate first-order derivatives
# From Fornberg and Sloan (Acta Numerica, 1994, p. 203-267; Table 1, page 213)
# Adapted by Scott Foster from code nicked off the net 2007
D.n <- length( x0)
if ( is.null( m)) {
D.f0 <- f(x0, ...)
m <- length( D.f0)
}
if ( D.accur == 2) {
D.w <- tcrossprod( rep( 1, m),c( -1/2, 1/2))
D.co <- c( -1, 1)
}
else {
D.w <- tcrossprod( rep( 1, m),c( 1/12, -2/3, 2/3, -1/12))
D.co <- c( -2, -1, 1, 2)
}
D.n.c <- length( D.co)
macheps <- .Machine$double.eps
D.h <- macheps^( 1/3)*abs( x0)
D.deriv <- matrix( NA, nrow=m, ncol=D.n)
for ( ii in 1:D.n) {
D.temp.f <- matrix( 0, m, D.n.c)
for ( jj in 1:D.n.c) {
D.xd <- x0+D.h[ii]*D.co[jj]*( 1:D.n == ii)
D.temp.f[,jj] <- f( D.xd, ...)
}
D.deriv[,ii] <- rowSums( D.w*D.temp.f)/D.h[ii]
}
return( as.double( D.deriv))
}
"nH2" <-
function( pt, fun, accur=c(4,4), type="H.Diag", ...) {
# A function to compute highly accurate second order Hessian diags and other off-diags
# partially from Fornberg and Sloan (Acta Numerica, 1994, p. 203-267; Table 1, page 213)
# Adapted by Scott Foster from code nicked off the net 2007
H.n <- length( pt)
derivs <- function( d.x0, ...) { nd2( x0=d.x0, f=fun, m=1, D.accur=accur[2], ...) }
Hes <- nd2( x0=pt, f=derivs, D.accur=accur[2], ...)
Hes <- matrix( Hes, nrow=length( pt))
Hes <- ( Hes+t( Hes))/2
if ( type == "H.Diag") {
macheps <- .Machine$double.eps
H.h <- macheps^(1/4)*abs(pt)
H.f0 <- fun( pt, ...)
H.m <- length( H.f0)
if ( accur[1] == 2) {
H.w <- tcrossprod( rep( 1, H.m), c( 1, -2, 1))
H.co <- c( -1, 0, 1)
}
else {
H.w <- tcrossprod( rep( 1, H.m), c( -1/12, 4/3, -5/2, 4/3, -1/12))
H.co <- c( -2, -1, 0, 1, 2)
}
H.n.c <- length( H.co)
Hes.diag <- double( length=H.n)
for ( ii in 1:H.n) {
H.temp.f <- matrix( 0, H.m, H.n.c)
for ( jj in 1:H.n.c) {
if ( H.co[jj] != 0) {
H.xd <- pt+H.h[ii]*H.co[jj]*( 1:H.n == ii)
H.temp.f[,jj] <- fun( H.xd, ...)
}
else
H.temp.f[,jj] <- H.f0
}
Hes.diag[ii] <- rowSums( H.w*H.temp.f)/( H.h[ii]^2)
}
diag( Hes) <- Hes.diag
}
return( Hes)
}
".onLoad" <-
function (libname, pkgname)
{
##require(MASS)
# Generic DLL loader
dll.path <- file.path( libname, pkgname, 'libs')
if( nzchar( subarch <- .Platform$r_arch))
dll.path <- file.path( dll.path, subarch)
this.ext <- paste( sub( '.', '[.]', .Platform$dynlib.ext,
fixed=TRUE), '$', sep='')
dlls <- dir( dll.path, pattern=this.ext, full.names=FALSE)
names( dlls) <- dlls
if( length( dlls))
lapply( dlls, function( x) library.dynam( sub( this.ext, '', x),
package=pkgname, lib.loc=libname))
}
"predict.archetype" <-
function (object,new.obs,...)
{
mixture.model <- object
if(class(mixture.model)[2]=="bernoulli"){
G <- length(mixture.model$pi)
covar <- mixture.model$covar[-(1:(G-1)),-(1:(G-1))]
coef <- mixture.model$coef
model.fm <- as.formula(mixture.model$formula)
model.fm[[2]] <- NULL
X <- model.matrix(model.fm,new.obs)
link.fun <- make.link("logit")
outvar <- matrix(NA,dim(X)[1],G)
outpred <- matrix(NA,dim(X)[1],G)
colnames(outvar) <- colnames(outpred) <- paste("G",1:G,sep=".")
for(g in 1:G){
lp <- as.numeric(X%*%coef[g,])
outpred[,g] <- link.fun$linkinv(lp)
dhdB <- (exp(lp)/(1+exp(lp)))*X - exp(lp)^2/((1+exp(lp))^2)*X
c2 <- covar[seq(g,dim(covar)[1],G),seq(g,dim(covar)[1],G)]
for(k in 1:dim(X)[1]){
outvar[k,g] <- (dhdB[k,]%*%c2)%*%(dhdB[k,])
}
}
}
if(class(mixture.model)[2]=="negbin"|class(mixture.model)[2]=="tweedie"){
G <- length(mixture.model$pi)
##covar <- mixture.model$covar[-1*c(1:(G-1),(dim(mixture.model$covar)[1]-2*length(mixture.model$sp.intercept)+1):dim(mixture.model$covar)[1]),-1*c(1:(G-1),(dim(mixture.model$covar)[1]-2*length(mixture.model$sp.intercept)+1):dim(mixture.model$covar)[1])] # remove pi,sp.int,sp.disp from matrix
covar <- mixture.model$covar[-1*c(1:(G-1),(dim(mixture.model$covar)[1]-length(mixture.model$sp.intercept)+1):dim(mixture.model$covar)[1]),-1*c(1:(G-1),(dim(mixture.model$covar)[1]-length(mixture.model$sp.intercept)+1):dim(mixture.model$covar)[1])] # remove pi,sp.disp from matrix
sp.int <- mixture.model$sp.intercept
coef <- mixture.model$coef
model.fm <- as.formula(mixture.model$formula)
model.fm[[2]] <- NULL
X <- cbind(model.matrix(model.fm,new.obs),1)
offset <- model.frame(model.fm, data = new.obs)
offset <- model.offset(offset)
if(is.null(offset)) offset <- rep(0,nrow(X))
outvar <- matrix(NA,dim(X)[1],G)
outpred <- matrix(NA,dim(X)[1],G)
colnames(outvar) <- colnames(outpred) <- paste("G",1:G,sep=".")
for(g in 1:G){
s.outvar <- matrix(NA,dim(X)[1],length(sp.int))
s.outpred <- matrix(NA,dim(X)[1],length(sp.int))
for(s in 1:length(sp.int)){
lp <- as.numeric(X%*%c(coef[g,],sp.int[s])+offset)
s.outpred[,s] <- exp(lp)
dhdB <- exp(lp)*X
c2 <- covar[c(seq(g,G*(dim(X)[2]-1),G),G*(dim(X)[2]-1)+s),c(seq(g,G*(dim(X)[2]-1),G),G*(dim(X)[2]-1)+s)]
for(k in 1:dim(X)[1]){
s.outvar[k,s] <- (dhdB[k,]%*%c2)%*%(dhdB[k,])
}
}
outpred[,g] <- apply(s.outpred*rep(mixture.model$tau[,g],each=dim(X)[1]),1,mean)/sum(mixture.model$tau[,g])
outvar[,g]<- apply(s.outvar*rep(mixture.model$tau[,g],each=dim(X)[1]),1,mean)/sum(mixture.model$tau[,g])
}
}
list(fit=outpred,se.fit=sqrt(outvar))
}
"print.archetype" <-
function (x,...)
{
cat("\nMixing probabilities\n")
print(x$pi)
cat("\nCoefficents\n")
print(x$coef)
if(!is.na(x$se[1])){
cat("\nStandard Errors of coefficents\n")
print(x$se)
}
cat("\nPosterior Probabilities\n")
print(x$tau)
}
"pTweedie" <-
function ( quant, mu, phi, p)
{
lambda <- ( mu^( 2-p)) / ( phi*(2-p))
alpha <- ( 2-p) / ( p-1)
tau <- phi*(p-1)*mu^(p-1)
mu.Z <- alpha * tau
ps <- 0
# ps <- pPoisGam( quant, lambda, mu.Z, alpha)
# require( tweedie)
# ps <- ptweedie( as.numeric( q), mu=as.numeric( mu), phi=as.numeric( phi), power=as.numeric( p))
return( ps)
}
"rPoisGam" <-
function ( n, lambda, mu.Z, alpha)
{
mu.N <- lambda
#simulate n random variables from the same compound poisson distribution
my.fun <- function (parms)
return( rgamma( n=1, scale=parms[3], shape=parms[1]*parms[2]))
# return( sum( rgamma( n=parms[1], scale=parms[3], shape=parms[2])))
tmp <- c( n, length( mu.N), length(mu.Z), length( alpha))
names( tmp) <- c( "n","mu.N","mu.Z","alpha")
if( !all( is.element( tmp[-1], c( 1, tmp[1])))) {
print( "rPoisGam: error -- length of arguments are not compatible")
return( tmp)
}
if( tmp["mu.N"]==1)
mu.N <- rep( mu.N, tmp["n"])
if( tmp["mu.Z"]==1)
mu.Z <- rep( mu.Z, tmp["n"])
if( tmp["alpha"]==1)
alpha <- rep( alpha, tmp["n"])
np <- matrix( rpois( n, mu.N), ncol=1)
beta <- mu.Z / alpha
y <- apply( cbind( np, alpha, beta), 1, my.fun)
return( y)
}
"rTweedie" <-
function ( n, mu, phi, p)
{
lambda <- ( mu^( 2-p)) / ( phi*(2-p))
alpha <- ( 2-p) / ( p-1)
tau <- phi*(p-1)*mu^(p-1)
mu.Z <- alpha * tau
rans <- rPoisGam( n, lambda, mu.Z, alpha)
return( rans)
}
"SpeciesMix" <-
function (sp.form,sp.data,covar.data,G=2, pars=NA, em.prefit=TRUE,em.steps=3, em.refit = 1 , dist="bernoulli",est.var = FALSE,residuals=FALSE,trace=TRUE)
{
if(dist=="bernoulli") return(SpeciesMix.bernoulli(sp.form,sp.data,covar.data,G, pars, em.prefit,em.steps, em.refit ,est.var,residuals,trace,r1=FALSE))
## if(dist=="bernoulli-random") return(SpeciesMix.bernoulli(sp.form,sp.data,covar.data,G, pars, em.prefit,em.steps, em.refit ,est.var,residuals,trace,r1=TRUE))
if(dist=="negbin") return(SpeciesMix.nbinom(sp.form,sp.data,covar.data,G, pars, em.prefit,em.steps, em.refit ,est.var,residuals,trace))
if(dist=="tweedie") return(SpeciesMix.tweedie(sp.form,sp.data,covar.data,G, pars, em.prefit,em.steps, em.refit ,est.var,residuals,trace))
##if(dist=="gaussian") return(SpeciesMix.gaussian(sp.form,sp.data,covar.data,G, pars, em.prefit,em.steps, em.refit ,est.var,residuals,trace))
print("incorrect distribution type, options are bernoulli, negbin or tweedie")
return(0)
}
"SpeciesMix.bernoulli" <-
function (sp.form,sp.data,covar.data,G=2, pars=NA, em.prefit=TRUE,em.steps=4, em.refit = 1 , est.var = FALSE,residuals=FALSE,trace=TRUE,r1=FALSE)
{
t.covar.data <- covar.data
t.sp.data <- sp.data
sp.form <- update.formula(sp.form,obs~1+.)
if(em.prefit | G==1){
prefit <- SpeciesMix.em(sp.form,sp.data,covar.data,G,ite.max=em.steps,em.refit=em.refit,r1)
pars <- c(additive.logistic(prefit$pi,T)[1:(G-1)],unlist(prefit$coef))
}
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
if(G==1) return(prefit)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=as.numeric(unlist(sp.data)),covar.data)
names(data)[1] <- as.character(sp.form)[2]
fmM.out <- fitmix.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var,r1)
while(fmM.out$logl==0) {
prefit <- SpeciesMix.em(sp.form,t.sp.data,t.covar.data,G,ite.max=em.steps,em.refit=1,r1)
pars <- c(additive.logistic(prefit$pi,T)[1:(G-1)],unlist(prefit$coef))
fmM.out <- fitmix.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var,r1)
}
rownames(fmM.out$tau) <- sp.name ## add names to taus
fmM.out$se <- NA
if(est.var){
fmM.out$covar <- try(solve(fmM.out$hessian))
if(class(fmM.out$covar)!="try-error"){
##colnames(fmM.out$covar) <- rownames(fmM.out$covar) <- names(fmM.out$gradient)
tmp <- sqrt(diag(fmM.out$covar))
tmp <- tmp[(G):length(tmp)]
fmM.out$se <- matrix(tmp,G,ncol(fmM.out$coef))
colnames(fmM.out$se) <- colnames(fmM.out$coef)
rownames(fmM.out$se) <- rownames(fmM.out$coef)
}
}
if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
fmM.out$formula <- sp.form
class(fmM.out) <- c("archetype","bernoulli")
fmM.out
}
"SpeciesMix.em" <-
function (sp.form,sp.data,covar.data,G=2,em.refit=1,ite.max=500, est.var = FALSE,residuals=FALSE,trace=TRUE,r1=FALSE)
{
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=unlist(sp.data),covar.data)
fmM.out <- fitMix(sp.form,data,sp,G,ite.max,trace,r1)
if(em.refit>1)
for(i in 2:em.refit){
fmM <- fitMix(sp.form,data,sp,G,ite.max,trace,r1)
if(fmM$logl>fmM.out$logl) fmM.out <- fmM
}
if(est.var){
var <- 0
t.pi <- additive.logistic(fmM.out$pi,TRUE)
parms <- c(t.pi[1:(G-1)],unlist(fmM.out$coef))
first.fit <- list(y=data[,1],x=model.matrix(sp.form,data=data))
fun.est.var <- function(x){-logLmix(x,first.fit,G,S,sp,sp.name)}
var <- solve( nH2( pt=parms, fun=fun.est.var))
colnames( var) <- rownames( var) <- names( parms)
fmM.out$covar <- var
}
if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
fmM.out$formula <- sp.form
fmM.out
}
"SpeciesMix.em.gaussian" <-
function (sp.form,sp.data,covar.data,G=2,em.refit=1,ite.max=500, est.var = FALSE,residuals=FALSE,trace=TRUE)
{
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=unlist(sp.data),covar.data)
fmM.out <- fitmix.gaussian(sp.form,data,sp,G,ite.max,trace)
if(em.refit>1)
for(i in 2:em.refit){
fmM <- fitmix.gaussian(sp.form,data,sp,G,ite.max,trace)
if(fmM$logl>fmM.out$logl) fmM.out <- fmM
}
##est.var <- FALSE
if(est.var){
var <- 0
t.pi <- additive.logistic(fmM.out$pi,TRUE)
##parms <- c(t.pi[1:(G-1)],unlist(fmM.out$coef))
parms <- c(t.pi[1:(G-1)],fmM.out$theta,unlist(fmM.out$coef),unlist(fmM.out$sp.intercept))
first.fit <- list(y=data[,1],x=model.matrix(sp.form,data=data)[,-1])
fun.est.var <- function(x){-logLmix.gaussian(x,first.fit,G,S,sp,sp.name)}
var <- solve( nH2( pt=parms, fun=fun.est.var))
colnames( var) <- rownames( var) <- names( parms)
fmM.out$covar <- var
}
residuals <- FALSE
if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
fmM.out$formula <- sp.form
fmM.out
}
"SpeciesMix.em.nbinom" <-
function (sp.form,sp.data,covar.data,G=2,em.refit=1,ite.max=500, est.var = FALSE,residuals=FALSE,trace=TRUE)
{
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=unlist(sp.data),covar.data)
fmM.out <- fitmix.nbinom(sp.form,data,sp,G,ite.max,trace)
if(em.refit>1)
for(i in 2:em.refit){
fmM <- fitmix.nbinom(sp.form,data,sp,G,ite.max,trace)
if(fmM$logl>fmM.out$logl) fmM.out <- fmM
}
##est.var <- FALSE
if(est.var){
var <- 0
t.pi <- additive.logistic(fmM.out$pi,TRUE)
##parms <- c(t.pi[1:(G-1)],unlist(fmM.out$coef))
parms <- c(t.pi[1:(G-1)],unlist(fmM.out$coef),fmM.out$sp.intercept,rep(1,S))##unlist(fmM.out$sp.intercept))
first.fit <- list(y=data[,1],x=model.matrix(sp.form,data=data)[,-1])
# first.fit$x <- cbind(sp.mat,first.fit$x)
fun.est.var <- function(x){-logLmix.nbinom(x,first.fit,G,S,sp,sp.name)}
#
deriv <- nd2(parms,fun.est.var)
var <- solve( nH2( pt=parms, fun=fun.est.var))
colnames( var) <- rownames( var) <- names( parms)
fmM.out$covar <- var
}
residuals <- FALSE
if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
fmM.out$formula <- sp.form
fmM.out
}
"SpeciesMix.em.tweedie" <-
function (sp.form,sp.data,covar.data,G=2,em.refit=1,ite.max=500, est.var = FALSE,residuals=FALSE,trace=TRUE)
{
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=unlist(sp.data),covar.data)
fmM.out <- fitmix.tweedie(sp.form,data,sp,G,ite.max,trace)
if(em.refit>1)
for(i in 2:em.refit){
fmM <- fitmix.tweedie(sp.form,data,sp,G,ite.max,trace)
if(fmM$logl>fmM.out$logl) fmM.out <- fmM
}
est.var <- FALSE
if(est.var){
var <- 0
t.pi <- additive.logistic(fmM.out$pi,TRUE)
##parms <- c(t.pi[1:(G-1)],unlist(fmM.out$coef))
parms <- c(t.pi[1:(G-1)],unlist(fmM.out$coef),unlist(fmM.out$sp.intercept),fmM.out$theta)
first.fit <- list(y=data[,1],x=model.matrix(sp.form,data=data)[,-1])
fun.est.var <- function(x){-logLmix.tweedie(x,first.fit,G,S,sp,sp.name)}
var <- solve( nH2( pt=parms, fun=fun.est.var))
colnames( var) <- rownames( var) <- names( parms)
fmM.out$covar <- var
}
residuals <- FALSE
if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
fmM.out$formula <- sp.form
fmM.out
}
"SpeciesMix.gaussian" <-
function (sp.form,sp.data,covar.data,G=2, pars=NA, em.prefit=TRUE,em.steps=4, em.refit = 1 , est.var = FALSE,residuals=FALSE,trace=TRUE)
{
t.covar.data <- covar.data
t.sp.data <- sp.data
sp.form <- update.formula(sp.form,obs~1+.)
em.steps=100
##if(em.prefit | G==1){
prefit <- SpeciesMix.em.gaussian(sp.form,sp.data,covar.data,G,ite.max=em.steps,est.var=est.var,em.refit=em.refit)
return(prefit)
##pars <- c(additive.logistic(prefit$pi,T)[1:(G-1)],unlist(prefit$coef))
##}
##S <- dim(sp.data)[2]
##if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
##n <- dim(sp.data)[1]
##sp <- rep(sp.name,each=n)
##if(G==1) return(prefit)
##var.names <- colnames(covar.data)
##covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
##names(covar.data) <- var.names
##sp.data <- as.data.frame(sp.data)
##data <- data.frame(obs=as.numeric(unlist(sp.data)),covar.data)
##names(data)[1] <- as.character(sp.form)[2]
##fmM.out <- fitmix.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var)
##while(fmM.out$logl==0) {
## prefit <- SpeciesMix.em(sp.form,t.sp.data,t.covar.data,G,ite.max=em.steps,em.refit=1)
## pars <- c(additive.logistic(prefit$pi,T)[1:(G-1)],unlist(prefit$coef))
##fmM.out <- fitmix.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var)
##}
##rownames(fmM.out$tau) <- sp.name ## add names to taus
##if(est.var){
## fmM.out$covar <- try(solve(fmM.out$hessian))
##}
## if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
##fmM.out$formula <- sp.form
##class(fmM.out) <- "archetype"
##fmM.out
}
"SpeciesMix.nbinom" <-
function (sp.form,sp.data,covar.data,G=2, pars=NA, em.prefit=TRUE,em.steps=3, em.refit = 1 , est.var = FALSE,residuals=FALSE,trace=TRUE)
{
t.covar.data <- covar.data
t.sp.data <- sp.data
sp.form <- update.formula(sp.form,obs~1+.)
#em.steps=100
if(em.prefit | G==1){
prefit <- SpeciesMix.em.nbinom(sp.form,sp.data,covar.data,G,ite.max=em.steps,est.var=FALSE,em.refit=em.refit)
##return(prefit)
pars <- prefit$pars#c(additive.logistic(prefit$pi,T)[1:(G-1)],unlist(prefit$coef))
}
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
if(G==1) return(prefit)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=as.numeric(unlist(sp.data)),covar.data)
names(data)[1] <- as.character(sp.form)[2]
sp.form <- update.formula(sp.form,obs~-1+.)
fmM.out <- fitmix.nbinom.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var)
#while(fmM.out$logl==0) {
# prefit <- SpeciesMix.em.nbinom(sp.form,sp.data,covar.data,G,ite.max=em.steps,est.var=est.var,em.refit=em.refit)
# pars <- prefit$pars
# fmM.out <- fitmix.nbinom.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var)
#}
rownames(fmM.out$tau) <- sp.name ## add names to taus
fmM.out$se <- NA
if(est.var){
fmM.out$covar <- try(solve(fmM.out$hessian))
if(class(fmM.out$covar)!="try-error"){
#colnames(fmM.out$covar) <- rownames(fmM.out$covar) <- names(fmM.out$gradient)
tmp <- sqrt(diag(fmM.out$covar))
tmp <- tmp[-1*(1:((G-1)))]
tmp <- tmp[-1*((length(tmp)-(2*S-1)):length(tmp))]
fmM.out$se <- matrix(tmp,G,ncol(fmM.out$coef))
colnames(fmM.out$se) <- colnames(fmM.out$coef)
rownames(fmM.out$se) <- rownames(fmM.out$coef)
}
##t1 <- sqrt(diag(fmM.out$covar))
##t1 <- t1[-(1:G-1)]
##fmM.out$se.coef <- t1[1:length(fmM.out$coef)]
##t1 <- t1[-(1:length(fmM.out$coef))]
##dim(fmM.out$se.coef) <- dim(fmM.out$coef)
##fmM.out$se.int <- t1[1:S]
##fmM.out$se.disp <- t1[-(1:S)]
}
if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
fmM.out$formula <- sp.form
class(fmM.out) <- c("archetype","negbin")
fmM.out
}
"SpeciesMix.tweedie" <-
function (sp.form,sp.data,covar.data,G=2, pars=NA, em.prefit=TRUE,em.steps=3, em.refit = 1 , est.var = FALSE,residuals=FALSE,trace=TRUE)
{
t.covar.data <- covar.data
t.sp.data <- sp.data
sp.form <- update.formula(sp.form,obs~1+.)
pars <- NA
if(em.prefit | G==1){
prefit <- SpeciesMix.em.tweedie(sp.form,sp.data,covar.data,G,ite.max=em.steps,est.var=est.var,em.refit=em.refit)
##return(prefit)
pars <- prefit$pars##pars <- c(additive.logistic(prefit$pi,T)[1:(G-1)],unlist(prefit$coef))
}
S <- dim(sp.data)[2]
if(is.null(colnames(sp.data))){sp.name <- 1:S} else {sp.name <- colnames(sp.data)}
n <- dim(sp.data)[1]
sp <- rep(sp.name,each=n)
if(G==1) return(prefit)
var.names <- colnames(covar.data)
covar.data <- data.frame(kronecker( rep( 1, S), as.matrix(covar.data)))
names(covar.data) <- var.names
sp.data <- as.data.frame(sp.data)
data <- data.frame(obs=as.numeric(unlist(sp.data)),covar.data)
names(data)[1] <- as.character(sp.form)[2]
sp.form <- update.formula(sp.form,obs~-1+.)
fmM.out <- fitmix.tweedie.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var)
##while(fmM.out$logl==0) {
## prefit <- SpeciesMix.em(sp.form,t.sp.data,t.covar.data,G,ite.max=em.steps,em.refit=1)
## pars <- c(additive.logistic(prefit$pi,T)[1:(G-1)],unlist(prefit$coef))
##fmM.out <- fitmix.cpp(sp.form,data,sp,G,pars=pars,calc.hes=est.var)
##}
rownames(fmM.out$tau) <- sp.name ## add names to taus
fmM.out$se <- NA
if(est.var){
##fmM.out$covar <- try(solve(fmM.out$hessian))
fmM.out$covar <- try(solve(fmM.out$hessian))
if(class(fmM.out$covar)!="try-error"){
colnames(fmM.out$covar) <- rownames(fmM.out$covar) <- names(fmM.out$gradient)
tmp <- sqrt(diag(fmM.out$covar))
tmp <- tmp[-1*(1:((G-1)))]
tmp <- tmp[-1*((length(tmp)-(2*S-1)):length(tmp))]
fmM.out$se <- matrix(tmp,G,ncol(fmM.out$coef))
colnames(fmM.out$se) <- colnames(fmM.out$coef)
rownames(fmM.out$se) <- rownames(fmM.out$coef)
}
}
if(residuals) fmM.out$residuals <- mix.residuals(fmM.out,sp.form,data,sp)
fmM.out$formula <- sp.form
class(fmM.out) <- c("archetype","tweedie")
fmM.out
}
"tglm" <-
function ( mean.form, data, wts=NULL, phi=NULL, p=NULL, inits=NULL, vcov=TRUE, residuals=TRUE, trace=1, iter.max=150)
{
if( is.null( wts))
wts <- rep( 1, nrow( data))
e<-new.env()
e$wts <- wts
environment(mean.form) <- e
temp.p <- model.frame( mean.form, data=as.data.frame( data), weights=wts)
y <- model.response( temp.p)
names( y) <- NULL
X.p <- model.matrix( mean.form, data)
offset.p <- model.offset( temp.p)
if ( is.null( offset.p))
offset.p <- rep( 0, nrow( temp.p))
wts1 <- as.vector( model.weights( temp.p))
fm1 <- NULL
if( is.null( inits)){
inits <- rep( 0, ncol( X.p))
if( trace!=0)
print( "Obtaining initial mean values from log-linear Poisson model -- this might be stupid")
abit <- 1
ystar <- round( y, digits=0)
fm1 <- glm( ystar~-1+X.p, family=poisson( link="log"), weights=wts1)
inits <- fm1$coef
}
if( is.null( phi) & length( inits)==ncol( X.p)){
if( is.null( fm1))
inits <- c( inits, 1)
else{
if( trace!=0)
print( "Obtaining initial dispersion from the smaller of the Pearson or Deviance estimator (or 25)")
if( is.null( p) & length( inits)==ncol( X.p))
ptemp <- 1.9
else
ptemp <- p
disDev <- fm1$deviance/( length( y) - length( fm1$coef))
disPear <- sum( ( wts1 * ( y-fm1$fitted)^2) / ( fm1$fitted^ptemp)) / ( length( y) - length( fm1$coef))
dis <- min( disDev, disPear, 25)
inits <- c( inits, dis)
}
}
if( is.null( p) & length( inits)==ncol( X.p) + is.null( phi))
inits <- c( inits, 1.6)
if( length( inits) != ncol( X.p) + is.null( phi) + is.null( p)) {
print( "Initial values supplied are of the wrong length -- please check")
tmp <- c( length( inits), ncol( X.p) + is.null( phi) + is.null( p))
names( tmp) <- c( "inits", "nParams")
return( tmp)
}
fmTGLM <- tglm.fit( x=X.p, y=y, wts=wts1, offset=offset.p, inits=inits, phi=phi, p=p, vcov=vcov, residuals=residuals, trace=trace, iter.max=iter.max)
return( fmTGLM)
}
"tglm.fit" <-
function ( x, y, wts=NULL, offset=rep( 0, length( y)), inits, phi=NULL, p=NULL, vcov=TRUE, residuals=TRUE, trace=1, iter.max=150)
{
if( trace!=0){
print( "Estimating parameters")
if( is.null( phi) & is.null( p))
cat("iter:", "-logl", colnames(x), "phi", "p", "\n", sep = "\t")
if( is.null( phi) & !is.null( p))
cat("iter:", "-logl", colnames(x), "phi", "\n", sep = "\t")
if( !is.null( phi) & is.null( p))
cat("iter:", "-logl", colnames(x), "p", "\n", sep = "\t")
if( !is.null( phi) & !is.null( p))
cat("iter:", "-logl", colnames(x), "\n", sep = "\t")
}
eps <- 1e-5
my.lower <- rep( -Inf, ncol( x))
my.upper <- rep( Inf, ncol( x))
if( is.null( phi))
{my.lower <- c( my.lower, eps); my.upper <- c( my.upper, Inf)}
if( is.null( p))
{my.lower <- c( my.lower, 1+eps); my.upper <- c( my.upper, 2-eps)}
fm <- nlminb( start=inits, objective=ldTweedie.lp, gradient=ldTweedie.lp.deriv, hessian=NULL, lower=my.lower, upper=my.upper, y=y, X.p=x, offsetty=offset, phi=phi, p=p, control=list(trace=trace, iter.max=iter.max), wts=wts)
parms <- fm$par
tmp <- colnames( x)
if( !is.null( phi) & !is.null( p))
tmp <- tmp
if( !is.null( phi) & is.null( p))
tmp <- c( tmp, "p")
if( is.null( phi) & !is.null( p))
tmp <- c( tmp, "phi")
if( is.null( phi) & is.null( p))
tmp <- c( tmp, "phi", "p")
names( parms) <- tmp
if( vcov){
if( trace!=0)
print( "Calculating variance matrix of estimates")
vcovar <- nd2(x0=parms, f=ldTweedie.lp.deriv, y=y, X.p=x, offsetty=offset, phi=phi, p=p, wts=wts)
vcovar <- 0.5 * ( vcovar + t( vcovar))
vcovar <- solve( vcovar)
rownames( vcovar) <- colnames( vcovar) <- names( parms)
}
else{
if( trace!=0)
print( "Not calculating variance matrix of estimates")
vcovar <- NULL
}
scores <- -ldTweedie.lp.deriv( parms=parms, y=y, X.p=x, offsetty=offset, phi=phi, p=p, wts=wts)
if( trace !=0)
print( "Calculating means")
mu <- exp( x %*% parms[1:ncol( x)] + offset)
if( residuals){
if( trace!=0)
print( "Calculating quantile residuals")
if( is.null( phi)){
phi1 <- parms[ncol( x)+1]
if( is.null( p))
p1 <- parms[ncol(x)+2]
else
p1 <- p
}
else{
phi1 <- phi
if( is.null( p))
p1 <- parms[ncol( x)+1]
else
p1 <- p
}
resids <- matrix( rep( pTweedie( y, mu, phi1, p1), 2), ncol=2)
resids[y==0,1] <- 0.5 * dTweedie( y[y==0], mu[y==0], phi1, p1, LOG=FALSE)
nzero <- sum( y==0)
resids[y==0,2] <- runif( nzero, min=rep( 0, nzero), max=2*resids[y==0,1])
resids <- qnorm( resids)
colnames( resids) <- c("expect","random")
}
else{
if( trace!=0)
print( "Not calculating quantile residuals")
resids <- NULL
}
if( trace!=0)
print( "Done")
ICadj <- 0
if( !is.null( phi))
ICadj <- ICadj+1
if( !is.null( p))
ICadj <- ICadj+1
AIC <- -2*(-fm$objective) + 2*(length( parms)-ICadj)
BIC <- -2*(-fm$objective) + log( nrow( x))*(length( parms)-ICadj)
res <- list( coef=parms, logl=-fm$objective, scores=scores, vcov=vcovar, conv=fm$convergence, message=fm$message, niter=fm$iterations, evals=fm$evaluations, call=match.call(), fitted=mu, residuals=resids, AIC=AIC, BIC=BIC)
class( res) <- "tglm"
return( res)
}
"weighted.glm" <-
function (g,first.fit,tau,n,fmM,sp)
{
dat.tau <- rep(tau[,g],each=n)
##lpre <- first.fit$x%*%fmM[[g]]$coef
##f.mix <- glm.fit(x=first.fit$x,y=first.fit$y,weights=dat.tau[,g],family=binomial(),etastart=fmM[[g]]$linear.predictors)
f.mix <- glm.fit(x=first.fit$x,y=first.fit$y,weights=dat.tau,family=binomial(),start=fmM[[g]]$coef)
return(list(coef=f.mix$coef))
##list(residuals=f.mix$residuals,fitted=f.mix$fitted,linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
##list(linear.predictors=f.mix$linear.predictors,coef=f.mix$coef)
}
"weighted.glm.gaussian" <-
function (g,first.fit,tau,n,fmM,sp)
{
dat.tau <- rep(tau[,g],each=n)
sp.name <- unique(sp)
X <- first.fit$x
sp.mat <- matrix(0,dim(X)[1],length(sp.name))
for(i in 1:length(sp.name)){
sp.mat[sp==sp.name[i],i] <- 1
}
X <- cbind(sp.mat,X)
f.mix <- glm.fit(x=X,y=first.fit$y,weights=dat.tau,family=gaussian())
#f.mix$theta <- 1
## f.mix <- glm.fit.nbinom(x=X,y=first.fit$y,weights=dat.tau)
sp.intercept <- f.mix$coef[1:length(sp.name)]
sp.intercept[is.na(sp.intercept)] <- 0
return(list(coef=f.mix$coef[-1:-length(sp.name)],theta=sqrt(f.mix$deviance/f.mix$df.residual),sp.intercept=sp.intercept,fitted=f.mix$fitted))#,lpre=f.mix$linear.predictors))
}
"weighted.glm.nbinom" <-
function (g,first.fit,tau,n,fmM,sp)
{
dat.tau <- rep(tau[,g],each=n)
sp.name <- unique(sp)
X <- first.fit$x
sp.mat <- matrix(0,dim(X)[1],length(sp.name))
for(i in 1:length(sp.name)){
sp.mat[sp==sp.name[i],i] <- 1
}
X <- cbind(sp.mat,X)
#f.mix <- glm.fit(x=X,y=first.fit$y,weights=dat.tau,family=poisson())
#f.mix$theta <- 1
f.mix <- glm.fit.nbinom(x=X,y=first.fit$y,offset=first.fit$offset,weights=dat.tau)
sp.intercept <- f.mix$coef[1:length(sp.name)]
sp.intercept[is.na(sp.intercept)] <- 0
return(list(coef=f.mix$coef[-1:-length(sp.name)],theta=f.mix$theta,sp.intercept=sp.intercept,fitted=f.mix$fitted))#,lpre=f.mix$linear.predictors))
}
"weighted.glm.tweedie" <-
function (g, first.fit, tau, n, fmM, sp)
{
dat.tau <- rep(tau[, g], each = n)
sp.int.offset <- rep(fmM[[g]]$sp.intercept, each = n)+first.fit$offset
sp.name <- unique(sp)
X <- first.fit$x
f.mix <- tglm.fit(x = X, y = first.fit$y, wts = dat.tau,
offset = sp.int.offset, vcov = FALSE, residuals = FALSE,
trace = 0, inits = fmM[[g]]$coef, phi = fmM[[g]]$phi,
p = 1.6)
return(list(coef = f.mix$coef, phi = fmM[[g]]$phi, p = 1.6,
sp.intercept = fmM[[g]]$sp.intercept, fitted = f.mix$fitted))
}
Try the SpeciesMix package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.