EnvelopeSort: Sorts Envelope function for simulation

Description Usage Arguments Details Value Examples

View source: R/EnvSort.R

Description

Sorts Enveloping function for simulation. how frequently each component of the resulting grid should be sampled during simulation.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
EnvelopeSort(
  l1,
  l2,
  GIndex,
  G3,
  cbars,
  logU,
  logrt,
  loglt,
  logP,
  LLconst,
  PLSD,
  a1
)

Arguments

l1

dimension for model (number of independent variables in X matrix)

l2

dimension for Envelope (number of components)

GIndex

matrix containing information on how each dimension should be sampled (1 means left tail of a restricted normal, 2 center, 3 right tail, and 4 the entire line)

G3

A matrix containing the points of tangencies associated with each component of the grid

cbars

A matrix containing the gradients for the negative log-likelihood at each tangency

logU

A matrix containing the log of the cummulative probability associated with each dimension

logrt

A matrix containing the log of the probability associated with the right tail (i.e. that to the right of the lower bound)

loglt

A matrix containing the log of the probability associated with the left tail (i.e., that to the left of the upper bound)

logP

A matrix containing log-probabilities related to the components of the grid

LLconst

A vector containing constant for each component of the grid used during the accept-reject procedure

PLSD

A vector containing the probability of each component in the Grid

a1

A vector containing the diagonal of the standardized precision matrix

Details

This function sorts the envelope in descending order based on the probability associated with each component in the Grid. Sorting helps speed up simulation once the envelope is constructed.

Value

The function returns a list consisting of the following components (the first six of which are matrics with number of rows equal to the number of components in the Grid and columns equal to the number of parameters):

GridIndex

A matrix containing information on how each dimension should be sampled (1 means left tail of a restricted normal, 2 center, 3 right tail, and 4 the entire line)

thetabars

A matrix containing the points of tangencies associated with each component of the grid

cbars

A matrix containing the gradients for the negative log-likelihood at each tangency

logU

A matrix containing the log of the cummulative probability associated with each dimension

logrt

A matrix containing the log of the probability associated with the right tail (i.e. that to the right of the lower bound)

loglt

A matrix containing the log of the probability associated with the left tail (i.e., that to the left of the upper bound)

LLconst

A vector containing constant for each component of the grid used during the accept-reject procedure

logP

A matrix containing log-probabilities related to the components of the grid

PLSD

A vector containing the probability of each component in the Grid

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
data(menarche2)

summary(menarche2)
plot(Menarche/Total ~ Age, data=menarche2)

Age2=menarche2$Age-13

x<-matrix(as.numeric(1.0),nrow=length(Age2),ncol=2)
x[,2]=Age2

y=menarche2$Menarche/menarche2$Total
wt=menarche2$Total

mu<-matrix(as.numeric(0.0),nrow=2,ncol=1)
mu[2,1]=(log(0.9/0.1)-log(0.5/0.5))/3

V1<-1*diag(as.numeric(2.0))

# 2 standard deviations for prior estimate at age 13 between 0.1 and 0.9
## Specifies uncertainty around the point estimates

V1[1,1]<-((log(0.9/0.1)-log(0.5/0.5))/2)^2 
V1[2,2]=(3*mu[2,1]/2)^2  # Allows slope to be up to 1 times as large as point estimate 

famfunc<-glmbfamfunc(binomial(logit))

f1<-famfunc$f1
f2<-famfunc$f2
f3<-famfunc$f3
f5<-famfunc$f5
f6<-famfunc$f6

dispersion2<-as.numeric(1.0)
start <- mu
offset2=rep(as.numeric(0.0),length(y))
P=solve(V1)
n=1000

## Appears that the type for some of these arguments are important/problematic

#outlist<-rnnorm_reg_cpp(n=as.integer(n),y=as.vector(y),x=as.matrix(x),
#mu=as.vector(mu),P=as.matrix(P),offset2=as.vector(offset2),
#wt=as.vector(wt),dispersion=as.numeric(dispersion2),famfunc=famfunc,
#f1=f1,f2=f2,f3=f3,start=as.vector(start),family="binomial",
#link="logit",Gridtype=as.integer(3))

### This allows use of the rglmb summary function 
### add interface for glmbsim_NGauss_cpp later

#outlist$call<-match.call()
#colnames(outlist$coefficients)<-colnames(x)
#class(outlist)<-c(outlist$class,"rglmb")
#summary(outlist)
#Env1=outlist$Envelope


###### Adjust weight for dispersion

wt2=wt/dispersion2

######################### Shift mean vector to offset so that adjusted model has 0 mean

alpha=x%*%as.vector(mu)+offset2
mu2=0*as.vector(mu)
P2=P
x2=x


#####  Optimization step to find posterior mode and associated Precision

parin=start-mu

opt_out=optim(parin,f2,f3,y=as.vector(y),x=as.matrix(x),mu=as.vector(mu2),
              P=as.matrix(P),alpha=as.vector(alpha),wt=as.vector(wt2),
              method="BFGS",hessian=TRUE
)

bstar=opt_out$par  ## Posterior mode for adjusted model
bstar
bstar+as.vector(mu)  # mode for actual model
A1=opt_out$hessian # Approximate Precision at mode

## Standardize Model

Standard_Mod=glmb_Standardize_Model(y=as.vector(y), x=as.matrix(x),
P=as.matrix(P),bstar=as.matrix(bstar,ncol=1), A1=as.matrix(A1))

bstar2=Standard_Mod$bstar2  
A=Standard_Mod$A
x2=Standard_Mod$x2
mu2=Standard_Mod$mu2
P2=Standard_Mod$P2
L2Inv=Standard_Mod$L2Inv
L3Inv=Standard_Mod$L3Inv

Env2=EnvelopeBuild(as.vector(bstar2), as.matrix(A),y, as.matrix(x2),
as.matrix(mu2,ncol=1),as.matrix(P2),as.vector(alpha),as.vector(wt2),
family="binomial",link="logit",Gridtype=as.integer(3), 
n=as.integer(n),sortgrid=TRUE)

## These now seem to match

#Env1
Env2

knygren/glmbayes documentation built on Sept. 4, 2020, 4:39 p.m.