R/us_matrix.R In wgeesel: Weighted Generalized Estimating Equations and Model Selection

```us_matrix=function(x_i,y_i,w_i,beta,rho,phi,corstr,family,z_i,r_i,alpha){
mat_est=list()
non_missing_length=length(na.omit(y_i))
if(family=="binomial"){
mu_i=exp(x_i%*%beta)/(1+exp(x_i%*%beta))
cons_cov=c(mu_i*(1-mu_i)) #n_i*1 matrix#
D_i=x_i*cons_cov #n_i*p#
}
else if (family=="gaussian"){
mu_i=x_i%*%beta
cons_cov=rep(1,length(mu_i)) #n_i*1 matrix#
D_i=x_i #n_i*p#
}
else if (family=="poisson"){
mu_i=exp(x_i%*%beta)
cons_cov=c(mu_i) #n_i*1 matrix#
D_i=x_i*cons_cov #n_i*p#
}
S_i=diag(cons_cov) #n_i*n_i#
W_i=diag(w_i)
if(corstr=="exchangeable"){
R_i=diag(1,length(mu_i))+matrix(rho,ncol=length(mu_i),nrow=length(mu_i))-diag(rho,length(mu_i))
}
if(corstr=="independence"){
R_i=diag(1,length(mu_i))
}
if(corstr=="ar1"){
H <- abs(outer(1:length(mu_i),1:length(mu_i), "-"))
R_i=rho^H
}
if(corstr=="unstructured"){
R_i=rho
}
V_i=phi*(sqrt(S_i)%*%R_i%*%sqrt(S_i))
U_i=(t(D_i)%*%solve(V_i)%*%W_i)[,1:non_missing_length]%*%as.matrix(((y_i-mu_i)[1:non_missing_length]),ncol=1)
lam_it=exp(z_i%*%alpha)/(1+exp(z_i%*%alpha))

logit_S_i_list=lapply(2:length(r_i),function(x){
res=0
if(r_i[x-1]==1)res=(r_i[x-1])*(r_i[x]-lam_it[x])*t(t(z_i[x,]))
return(res)
})
logit_S_i=Reduce("+",logit_S_i_list)
mat_est[[1]]=U_i
mat_est[[2]]=logit_S_i
mat_est[[3]]=U_i%*%t(logit_S_i)
mat_est[[4]]=logit_S_i%*%t(logit_S_i)
mat_est[[5]]=D_i
mat_est[[6]]=W_i
mat_est[[7]]=V_i
return(mat_est)
}
```

Try the wgeesel package in your browser

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

wgeesel documentation built on May 2, 2019, 3:41 a.m.