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
#out<-rglmb(n = 1000, y=y, x=x, mu=mu, P=solve(V1), wt = wt,
# family = binomial(logit), Gridtype = 3)
#summary(out)
famfunc<-glmbfamfunc(binomial(logit))
f1<-famfunc$f1
f2<-famfunc$f2 # Used in optim and glmbsim_cpp
f3<-famfunc$f3 # Used in optim
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
###### 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
Env2
#int n, NumericVector y, NumericMatrix x, NumericMatrix mu, NumericMatrix P,
#NumericVector alpha, NumericVector wt, Function f2,
#Rcpp::List Envelope, Rcpp::CharacterVector family,
#Rcpp::CharacterVector link, int progbar
### Note: getting the types correct here is important but potentially difficult for users
### May be better to call an R function wrapper that checks (and converts when possible)
## to correct types
sim=rnnorm_reg_std(n=as.integer(n),y=as.vector(y),x=as.matrix(x2),mu=as.matrix(mu2,ncol=1),
P=as.matrix(P2),alpha=as.vector(alpha),wt=as.vector(wt2),
f2=f2,Envelope=Env2,family="binomial",link="logit",as.integer(0))
out=L2Inv%*%L3Inv%*%t(sim$out)
for(i in 1:n){
out[,i]=out[,i]+mu
}
summary(t(out))
mean(sim$draws)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.