mh_gibbs_modular: mh_gibbs_modular Function

In nmolanog/bayesuirt:

Description

this functions implements the metropolis-hastings within gibbs algorithm proposed in molano thesis for a given model.

Usage

mh_gibbs_modular(dat, idname, psi0, theta0, B = NULL, b = NULL, grad.fix,
  grad.mix, nlf, inv_link, hprim, Nc = 10000)

Arguments

dat	a data frame containing test info. rows for individuals and columns for items. aditionally, a variable which identify individuals must be supied.
idname	atomic character giving the name of the id variable in dat.
psi0	a vector of initial values for parameter psi0. See nlf_2pl documentation.
theta0	vector of initial values for latent traits. See nlf_2pl documentation.
B	Prior covariance matrix for psi.
b	Prior expected values for psi.
grad.fix	gradient function of fixed effects
grad.mix	gradient function of random effects
nlf	nonlinear function
inv_link	inverse of link function
hprim	derivative of link function
Nc	number of mcmc iterations.

Value

a list with the following objects:

• fixp: a matrix where each row corresponds to an mcmc simulation of psi parameters.

• theta: a matrix where each row corresponds to an mcmc simulation of latent traits.

• Dev: a vector where deviance is calculated for each mcmc iteration.

• B: Prior covariance matrix used for psi.

• b: Prior Expected values used for psi.

Examples

###Simulate a data set with model m1 of Yuli
nitems<-7
nind<-100
sort(alpha_o<-rnorm(nitems))
sort(beta_o<-rnorm(nitems))
summary(theta_o<-rnorm(nind))

temp.b<-matrix(beta_o,ncol=nitems,nrow=nind, byrow =T)
temp.a<-matrix(alpha_o,ncol=nitems,nrow=nind, byrow =T)
temp.t<-matrix(theta_o,ncol=nitems,nrow=nind)
logit_expt_theta_M<-exp(temp.t)/(1+exp(temp.t))
exp_alpha_M<-exp(temp.a)
exp_beta_M<-exp(temp.b)
### m1 model of yuli
pmod<-1-(1-(logit_expt_theta_M)^(exp_alpha_M))^(exp_beta_M)

summary(pmod)
test<-data.frame(ifelse(runif(nitems*nind)<pmod,1,0))
colnames(test)<-paste("t",1:nitems,sep="")
test%>%{apply(., 1,sum)/nitems} %>%summary
test%>%apply( 2,sum)/nind
test<-data.frame(test,id=1:nind)

###non linear function
nlf_m1<-function(X,psi,z,theta,nitems){
 logit_zt<-exp(z%*%theta)/(1+exp(z%*%theta))
 exp_xa<-exp(X%*%psi[(nitems+1):(2*nitems)])
 exp_xb<-exp(X%*%psi[1:nitems])
 return(drop(1-(1-(logit_zt)^(exp_xa))^(exp_xb)))
 }
 ###inverse link
 inv_link_m1<-function(x){x}
 ###fixed effects gradient
 grad.fix_m1<-function(X,psi,z,theta,nitems){
   logit_zt<-exp(z%*%theta)/(1+exp(z%*%theta))
   exp_xa<-exp(X%*%psi[(nitems+1):(2*nitems)])
   exp_xb<-exp(X%*%psi[1:nitems])
   logitzt_high_expa<-logit_zt^(exp_xa)
   a_grad<-drop(exp_xa*exp_xb*(logitzt_high_expa)*log(logit_zt)*((1-logitzt_high_expa)^(exp_xb-1)))*X
   b_grad<- drop(-exp_xb*log(1-(logitzt_high_expa))*((1-(logit_zt^exp_xa))^(exp_xb)))*X
   return(cbind(b_grad,a_grad))
 }
 ###random effects gradient
 grad.mix_m1<-function(X,psi,z,theta,nitems){
   logit_zt<-exp(z%*%theta)/(1+exp(z%*%theta))
   exp_xa<-exp(X%*%psi[(nitems+1):(2*nitems)])
   exp_xb<-exp(X%*%psi[1:nitems])
   logitzt_high_expa<-logit_zt^(exp_xa)
   return(drop(exp_xa*exp_xb*logitzt_high_expa*((1-logitzt_high_expa)^(exp_xb-1))*(1-logit_zt))*z)
 }

 ### derivative of link
 hprim_m1<-function(x){rep(1,length(x))}

 bres<-mh_gibbs_modular(test,"id",psi0=c(beta_o,alpha_o),theta0=theta_o,B=NULL,b=NULL,
grad.fix=grad.fix_m1,grad.mix=grad.mix_m1,nlf=nlf_m1,inv_link=inv_link_m1,
hprim=hprim_m1,Nc=10000)

nmolanog/bayesuirt documentation built on May 23, 2019, 8:52 a.m.