R/dynrMi.R

Defines functions dynr.mi

Documented in dynr.mi

##' Multiple Imputation of dynrModel objects
##' 
##' @param model dynrModel object
##' @param m number of multiple imputations
##' @param aux.variable names of auxiliary variables used in imputation
##' @param imp.obs logical. whether to impute the observed variables
##' @param imp.exo logical. whether to impute the exogenous variables
##' @param lag numeric. the number of lags to use
##' 
##' @details
##' This function is in alpha-testing form.  Please do not use or rely on it for now. A full implementation is in progress.
dynr.mi <- function(model, m=5, aux.variable, imp.obs=FALSE, imp.exo=FALSE, lag){    #multiple lag; #factor  #get variable names
	
	data <- model$data$original.data
	k <- length(coef(model)) #- length([email protected]@paramnames)             #number of parameters estimated
	nolag <- TRUE
	
	
	ynames <- model$data$observed.names
	xnames <- model$data$covariate.names
	y <- data[, colnames(data)==ynames]
	x <- data[, colnames(data)==xnames]
	ID <- model$data$id
	time <- model$data$time
	
	au <- data[, colnames(data)==aux.variable]  
	
	dataforlag <- cbind(ID,y,x)
	
	#dataforlag <- ts(dataforlag)
	
	#head(stats::lag(as.xts(dataforlag), 2))
	#head(lag(dataforlag, 2))
	
	datalag <- dataforlag
	#  dataforlag %>%
	#  dplyr::group_by(ID) %>%
	#  dplyr::mutate_all(lag) #Linying, what does this "lag" do?
	
	dataformice <- cbind(dataforlag[,-1], datalag[,-1], au)
	dataformice <- data.frame(dataformice)
	
	colnames(dataformice)=c()
	m <- m
	imp <- mice::mice(dataformice, m=m)
	
	
	pmcarqhat <- matrix(NA, nrow=m, ncol=k) #parameter estimates from each imputation
	pmcaru <- array(NA, dim=c(k,k,m)) #vcov of par estimates from each imputation
	
	for(j in 1:m){
		
		completedata <- mice::complete(imp, action=j)
		colnames(completedata) <- c(ynames, xnames,
		                      paste0("lag", ynames),
		                      paste0("lag", xnames),
		                      colnames(au))
		
		if(imp.obs==TRUE){
			imp.data.obs <- completedata[, colnames(completedata)==ynames]
		} else{
			imp.data.obs <- y
		}
		
		if(imp.exo==TRUE){
			imp.data.exo <- completedata[, colnames(completedata)==xnames]
		} else{
			imp.data.exo <- x
		}
		
		newdata <- cbind(ID, time, imp.data.obs, imp.data.exo)
		
		#Whoa WFT is this?!?
		save(newdata,file="test.rdata")
		#Whoa WFT is this?!?
		
		colnames(newdata) <- c("ID", "Time", "wp", "hp", "ca", "cn")
		
		data <- dynr.data(newdata, id="ID", time="Time",
		                observed=c("wp", "hp"), covariates=c("ca","cn"))
		
		modelnew <- model
		modelnew$data <- data
		#   model <- dynr.model([email protected], [email protected],
		#                      [email protected], [email protected], data=data, [email protected],
		#                     outfile=paste("trial4",i,".c",sep=""))
		
		
		# model <- dynr.model(dynamics=dynm, measurement=meas,
		#                     noise=mdcov, initial=initial, data=data,#transform=trans,
		#                     outfile=paste("trial2",i,".c",sep=""))
		
		
		trial <- dynr.cook(modelnew)  #names(trial) get names of the params
		#summary(trial)
		
		#getting parameter estimates
		pmcarqhat[j,] <- coef(trial)
		pmcaru[, ,j] <- vcov(trial)
	}
	
	pqbarmcarimpute <- apply(pmcarqhat, 2, mean) 
	pubarmcarimpute <- apply(pmcaru, 1:2, mean)
	#ubar <- apply(u, c(2, 3), mean)
	pe.mcarimpute <- pmcarqhat - matrix(pqbarmcarimpute, nrow = m, ncol = k, byrow = TRUE)
	pb.mcarimpute <- (t(pe.mcarimpute) %*% pe.mcarimpute)/(m - 1)
	pvcovmcarimpute <- pubarmcarimpute + (1 + 1/m) * pb.mcarimpute #vcov for estimates
	psemcarimpute <- sqrt(diag(pvcovmcarimpute))
	
	t <- pqbarmcarimpute/psemcarimpute
	# TODO Don't just use 2 for the CIs !!!
	ci.upper <- pqbarmcarimpute + 2*psemcarimpute
	ci.lower <- pqbarmcarimpute - 2*psemcarimpute
	p <- pt( abs(t), df=nobs(model) - k, lower.tail=FALSE)
	
	result <- cbind(pqbarmcarimpute, psemcarimpute, t, ci.lower, ci.upper,p)
	
	colnames(result) <- c("Estimate", "Std. Error", "t value", "ci.lower", "ci.upper", #"",
	                     "Pr(>|t|)")
	row.names(result) <- names(coef(model))
	
	#TODO the current 'result' should be what is printed by summary() on the thing returned from dynr.mi()
	# This summary should use methods similar to summary.dynrCook and summary.lm
	# TODO return a 'mi' object similar to rest of mice
	return(result)
}

Try the dynr package in your browser

Any scripts or data that you put into this service are public.

dynr documentation built on Sept. 25, 2018, 9:05 a.m.