Nothing
#TODO what is this function???
#Is it the empirical likelihood standard error? What assumptions?
#was EL_se
# @export
EL.se<-function(weights,outcome,N=NULL, use.second.order=TRUE, homophily=1){
nas <- is.na(weights)|is.na(outcome)
nas <- nas | is.infinite(weights)|is.infinite(outcome)
if(use.second.order){
wij <- outer(weights[!nas],weights[!nas],"*")
wij <- wij * (1+(homophily-1)*outer(outcome[!nas],outcome[!nas],"!="))
#
wij[wij > 10*stats::median(wij)] <- 10*stats::median(wij)
# Compute the EL standard error (if available)
iLj <- row(wij)>col(wij)
oij <- 0.5*outer(outcome[!nas],outcome[!nas],"+")
estij <- sum(oij*wij*iLj)/sum(wij*iLj)
Gi <- apply((wij*wij*(oij-estij))*iLj,1,sum)
wvi <- apply((wij*wij)*iLj,1,sum)
varest <- 4*sum(Gi^2)/(sum(wvi)^2)
}else{
wi <- weights[!nas]
oi <- outcome[!nas]
#
wi[wi > 10*stats::median(wi)] <- 10*stats::median(wi)
# Compute the EL standard error (if available)
esti <- sum(oi*wi)/sum(wi)
varest <- sum(((oi-esti)*wi*wi)^2)/(sum(wi*wi)^2)
}
if(is.null(N) | !is.numeric(N)){
sqrt( varest )
}else{
sqrt( (N-length(outcome))*varest/(N-1) )
}
}
# @export
EL.se.new<-function(weights,outcome,N=NULL, use.second.order=FALSE, homophily=1){
nas <- is.na(weights)|is.na(outcome)
nas <- nas | is.infinite(weights)|is.infinite(outcome)
if(use.second.order){
wi <- weights[!nas]
wi[wi > 10*stats::median(wi)] <- 10*stats::median(wi)
wij <- outer(weights[!nas],weights[!nas],"*")
# wij <- wij * (1+(homophily-1)*outer(outcome[!nas],outcome[!nas],"!="))
wij[wij > 10*stats::median(wij)] <- 10*stats::median(wij)
# Compute the EL standard error (if available)
iLj <- row(wij) >col(wij)
iLjn <- row(wij)!=col(wij)
oij <- 0.5*outer(outcome[!nas],outcome[!nas],"+")
diag(oij) <- 0
estij <- sum(oij*wij*iLj)/sum(wij*iLj)
#
qij <- 2*(oij-estij)*wij*wij*iLjn
q <- sum(wij*wij*iLjn)
qi <- apply(qij,2,sum)
cross.mat=1-wij/outer(wi,wi,"*")
sgn=sign(sum((cross.mat/wij)[iLjn]))
varest <- (sum(qij*qij)+sgn*sum(cross.mat*outer(qi,qi,"*")*iLjn))/(q*q)
}else{
wij <- outer(weights[!nas],weights[!nas],"*")
# wij <- wij * (1+(homophily-1)*outer(outcome[!nas],outcome[!nas],"!="))
wij[wij > 10*stats::median(wij)] <- 10*stats::median(wij)
# Compute the EL standard error (if available)
iLj <- row(wij) >col(wij)
iLjn <- row(wij)!=col(wij)
oij <- 0.5*outer(outcome[!nas],outcome[!nas],"+")
diag(oij) <- 0
estij <- sum(oij*wij*iLj)/sum(wij*iLj)
qij <- 2*(oij-estij)*wij*wij*iLjn
q <- sum(wij*wij*iLjn)
varest <- sum(qij*qij)/(q*q)
}
# if(is.null(N) | !is.numeric(N)){
sqrt( varest )
# }else{
# sqrt( (N-length(outcome))*varest/(N-1) )
# }
}
EL.est<-function(weights,outcome,N=NULL, use.second.order=TRUE, homophily=1.35){
nas <- is.na(weights)|is.na(outcome)
nas <- nas | is.infinite(weights)|is.infinite(outcome)
if(use.second.order){
wij <- outer(weights[!nas],weights[!nas],"*")
wij <- wij / (1+(homophily-1)*outer(outcome[!nas],outcome[!nas],"=="))
#
# wij[wij > 10*stats::median(wij)] <- 10*stats::median(wij)
# Compute the EL standard error (if available)
iLj <- row(wij)>col(wij)
oij <- 0.5*outer(outcome[!nas],outcome[!nas],"+")
estij <- sum(oij*wij*iLj)/sum(wij*iLj)
Gi <- apply((wij*wij*(oij-estij))*iLj,1,sum)
wvi <- apply((wij*wij)*iLj,1,sum)
varest <- 4*sum(Gi^2)/(sum(wvi)^2)
est <- estij
}else{
wi <- weights[!nas]
oi <- outcome[!nas]
#
wi[wi > 10*stats::median(wi)] <- 10*stats::median(wi)
# Compute the EL standard error (if available)
esti <- sum(oi*wi)/sum(wi)
varest <- sum(((oi-esti)*wi*wi)^2)/(sum(wi*wi)^2)
est <- esti
}
if(is.null(N) | !is.numeric(N)){
attr(est,"EL.se") <- sqrt( varest )
}else{
attr(est,"EL.se") <- sqrt( (N-length(outcome))*varest/(N-1) )
}
est
}
#' Get Horvitz-Thompson estimator assuming inclusion probability proportional
#' to the inverse of network.var (i.e. degree).
#' @param rds.data An rds.data.from
#' @param group.variable The grouping variable.
#' @param network.var The network.size variable.
#' @export
get.h.hat <- function(rds.data,
group.variable,
network.var=attr(rds.data,"network.size")){
if(is(rds.data,"rds.data.frame")){
stopifnot(group.variable %in% names(rds.data))
stopifnot(network.var %in% names(rds.data))
}
else{
stop("rds.data must be of type rds.data.frame")
}
harmonic.mean <- function(x){
x <- x[!is.na(x)]
length(x)/sum(1/x)
}
degrees <- rds.data[[network.var]]
temp <- function(x){
deg <- degrees[as.character(rds.data[[group.variable]]) == x]
return(harmonic.mean(deg))
}
group.names <- sort(unique(as.character(rds.data[[group.variable]])))
group.names <- group.names[!is.na(group.names)]
h.hat <- sapply(group.names,temp)
names(h.hat) <- group.names
return(h.hat)
}
# Compute the horvitz thompson estimate
HT.estimate<-function(weights,outcome){
nas <- is.na(weights)|is.na(outcome)
if(is.factor(outcome)){
num <- rep(0,length=length(levels(outcome[!nas])))
a <- tapply(weights[!nas],as.numeric(outcome[!nas]),sum)
num[as.numeric(names(a))] <- a
}else{
num<-sum(outcome[!nas]*weights[!nas])
}
den<-sum(weights[!nas])
num/den
}
#' Counts the number or recruiter->recruitee transitions between different levels
#' of the grouping variable.
#' @param rds.data An rds.data.frame
#' @param group.variable The name of a categorical variable in rds.data
#' @export
#' @examples
#' data(faux)
#' count.transitions(faux,"X")
count.transitions <- function(rds.data, group.variable)
{
stopifnot(is.rds.data.frame(rds.data))
#extract data
id <- get.id(rds.data)
recruiter.id <- get.rid(rds.data)
grp <- as.factor(rds.data[[group.variable]])
# Make sure the factor labels are alphabetic!
grp=factor(grp,levels=levels(grp)[order(levels(grp))])
#count transitions
ri <- match(recruiter.id,id)
rgrp <- grp[ri]
res <- table(rgrp,grp)
for(i in 1:nrow(res)){
if(sum(res[i,])==0){res[i,] <- 0.1}
}
class(res) <- "matrix"
res
}
# This function assumes that transition counts can be turned into MLEs.
# This will not be possible if any of the rows in the transition counts
# matrix have sum zero.
# IF 7/12/13: not sure what the point of this function is, I simplified the
# code, but why not just use prop.table
#' calculates the mle. i.e. the row proportions of the transition matrix
#' @param transition.counts a matrix or table of transition counts
#' @details depreicated. just use prop.table(transition.counts,1)
transition.counts.to.Markov.mle <- function(transition.counts){
#return(prop.table(tij,1))
tij <- transition.counts + .Machine$double.eps
ti <- margin.table(tij, 1)
stopifnot(ti>0)
sweep(tij, 1, ti, "/", check.margin = FALSE)
}
#' Markov chain statistionary distribution
#' @param mle The transition probabilities
#' @export
#' @return A vector of proportions representing the proportion in each group at the
#' stationary distribution of the Markov chain.
get.stationary.distribution <- function(mle){
stat.dist <- NULL
eigen.analysis <- eigen(t(mle))
# If -1 is an eigenvalue then we have a periodic chain ... which is bad.
values.close.to.minus.one <- abs(eigen(t(mle))$values + 1.0) <
sqrt(.Machine$double.eps)
if(sum(values.close.to.minus.one) == 1){
cat(" Markov chain appears to be periodic.\n")
stat.dist <- rep(NA,nrow(mle))
attr(stat.dist,'status') <- "Periodic Markov MLE"
return(stat.dist)
}
# The stationary distribution will be given by
# the eigenvector with eigenvalue equal to one.
values.close.to.one <- abs(eigen(t(mle))$values - 1.0) <
sqrt(.Machine$double.eps)
# If there is exactly one of these, then we can take it as eigenvalue corresponding to the
# stationary distribution.
if(sum(values.close.to.one) == 1){
stationary.dist.eigen.vector <- eigen.analysis$vectors[,which(values.close.to.one)]
stationary.dist.eigen.vector <- Re(stationary.dist.eigen.vector)
stat.dist <- stationary.dist.eigen.vector/sum(stationary.dist.eigen.vector)
attr(stat.dist,'status') <- "MLE gives unique stationary distribution."
names(stat.dist) <- colnames(mle)
}
return(stat.dist)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.