MVAdrug3waysTab (R 2.14.1)

MVAdrug3waysTab presents the coefficient estimates based on the saturated model with first and second order interactions as well as the coefficient estimates for the restricted model with first interactions for the drug data. It computes G-squared deviance and Chi-squared test statistics.

Download File

Awdesch Melzer

Applied Multivariate Statistical Analysis

Fri, May 04 2012 by Dedy Dwi Prastyo

-


Description: -


rm(list=ls(all=TRUE))
graphics.off()
install.packages("MASS")
library(MASS)
zi = rbind(c(1,0,1,0,1,0,0,0,0,21),c(1,0,1,0,0,1,0,0,0,32),c(1,0,1,0,0,0,1,0,0,70),c(1,0,1,0,0,0,0,1,0,43),c(1,0,1,0,0,0,0,0,1,19),c(1,0,0,1,1,0,0,0,0,683),c(1,0,0,1,0,1,0,0,0,596),c(1,0,0,1,0,0,1,0,0,705),c(1,0,0,1,0,0,0,1,0,295),c(1,0,0,1,0,0,0,0,1,99),c(0,1,1,0,1,0,0,0,0,46),c(0,1,1,0,0,1,0,0,0,89),c(0,1,1,0,0,0,1,0,0,169),c(0,1,1,0,0,0,0,1,0,98),c(0,1,1,0,0,0,0,0,1,51),c(0,1,0,1,1,0,0,0,0,738),c(0,1,0,1,0,1,0,0,0,700),c(0,1,0,1,0,0,1,0,0,847),c(0,1,0,1,0,0,0,1,0,336),c(0,1,0,1,0,0,0,0,1,196))
y = zi[,10]
I=2  # sex M - F
J=2  # drug Yes - No
K=5  # age category 16-29, 30-44, 45-64, 65-74, 75++
X = rbind(c(1,1,1,0,0,0),c(1,1,0,1,0,0),c(1,1,0,0,1,0),c(1,1,0,0,0,1),c(1,1,-1,-1,-1,-1),c(1,-1,1,0,0,0),c(1,-1,0,1,0,0),c(1,-1,0,0,1,0),c(1,-1,0,0,0,1),c(1,-1,-1,-1,-1,-1),c(-1,1,1,0,0,0),c(-1,1,0,1,0,0),c(-1,1,0,0,1,0),c(-1,1,0,0,0,1),c(-1,1,-1,-1,-1,-1),c(-1,-1,1,0 ,0,0),c(-1,-1,0,1,0,0),c(-1,-1,0,0,1,0),c(-1,-1,0,0,0,1),c(-1,-1,-1,-1,-1,-1))
n = dim(X)
n1 = n[1]
n2 = n[2]
X1 = cbind(X,X[,1]*X[,2],X[,1]*X[,3],X[,1]*X[,4],X[,1]*X[,5],X[,1]*X[,6],X[,2]*X[,3],X[,2]*X[,4],X[,2]*X[,5],X[,2]*X[,6])
X2 = cbind(X1,X[,1]*X[,2]*X[,3],X[,1]*X[,2]*X[,4],X[,1]*X[,2]*X[,5],X[,1]*X[,2]*X[,6])
X3 = cbind(rep(1,n1),X2)
nn = dim(X3)
nn1 = nn[1]
nn2 = nn[2]
df = nn1-nn2
df
b0 = solve(t(X3)%*%X3)%*%t(X3)%*%log(y)
b0
loglik = sum((X3%*%b0)*y)
loglik
nn = dim(X1)
nn1 = nn[1]
nn2 = nn[2]
XX = cbind(rep(1,nn1),X1)
nn = dim(XX)
nn1 = nn[1]
nn2 = nn[2]
df = nn1-nn2
df
fit = glm.nb(y~X1)
b = fit$coefficients
cbind(b)
loglik = sum((XX%*%b)*y)
loglik
lnyfit = XX%*%b
yfit = exp(lnyfit)
e = log(y)-lnyfit
print('degree of freedom')
print(df)
G2 = 2*sum(y*e)
pvalG2 = 1 - pchisq(G2,df)
X2 = sum(((y-yfit)^2)/yfit)
pvalG2 = 1 - pchisq(G2,df)
statstable = cbind(G2, pvalG2, X2, pvalG2, df)
print('Statistics: G2 | pvalue | chisq| pvalue | degrees of freedom')
print(statstable)
print(' ')
print(' observed fitted')
print(' values   values')
print(cbind(y,  yfit))
print(' ')