TABLE OF CONTENTS
estimation/updatedisppears [ Functions ]
NAME
updatedisppears --- Pearson dispersion parameter estimators
FUNCTION
Compute dispersion parameter estimates as Pearson-type estimators with a bias correction for the BLUP shrinkage effect.
SYNOPSIS
825 updatedisppears <- function(m, Ji, frailtyoutput, dispparams, corrval)
INPUTS
m number of clusters
Ji cluster sizes
frailtyoutput all output of fupdatefrailties4
dispparams a list of current estimates of dispersion parameters,
with components:
sigma2hat, sigma2dhat, nu2hat, nu2dhat, thetahat
corrval a ``correction factor'' that can be used to implement
Ma's degree-of-freedom adjustment
OUTPUTS
sigma2hat cluster-level dispersion for process 1
sigma2dhat cluster-level dispersion for process 2
nu2hat subject-level dispersion for process 1
nu2dhat subject-level dispersion for process 2
thetahat subject-level covariance
SOURCE
828 { 829 830 Jicum <- c(0, cumsum(Ji)) 831 jimax <- max(Ji) 832 833 ## Extract variables from the parameters passed into the function 834 pij <- frailtyoutput$pqrs$pij; qij <- frailtyoutput$pqrs$qij; 835 rij <- frailtyoutput$pqrs$rij; sij <- frailtyoutput$pqrs$sij; 836 pijprime <- frailtyoutput$pqrs$pijprime; qijprime <- frailtyoutput$pqrs$qijprime; 837 rijprime <- frailtyoutput$pqrs$rijprime; sijprime <- frailtyoutput$pqrs$sijprime; 838 pi <- frailtyoutput$pqrs$pi; qi <- frailtyoutput$pqrs$qi; 839 ri <- frailtyoutput$pqrs$ri; si <- frailtyoutput$pqrs$si; 840 piprime <- frailtyoutput$pqrs$piprime; qiprime <- frailtyoutput$pqrs$qiprime; 841 riprime <- frailtyoutput$pqrs$riprime; siprime <- frailtyoutput$pqrs$siprime; 842 wij <- frailtyoutput$pqrs$wij;zij <- frailtyoutput$pqrs$zij; 843 wi <- frailtyoutput$pqrs$wi 844 Uihat <- frailtyoutput$Uihat;Vihat <- frailtyoutput$Vihat; 845 Uijhat <- frailtyoutput$Uijhat;Vijhat <- frailtyoutput$Vijhat; 846 sigma2hat <- dispparams$sigma2hat; sigma2dhat <- dispparams$sigma2dhat 847 nu2hat <- dispparams$nu2hat; nu2dhat <- dispparams$nu2dhat 848 thetahat <- dispparams$thetahat 849 850 851 ## Initialize mean squared distance vectors 852 ci <- rep(0, m);cid <- rep(0, m); 853 cij <- matrix(0, m, jimax);cijd <- matrix(0, m, jimax); 854 cijprime <- matrix(0, m, jimax) 855 kUU <- matrix(0, m, jimax);kUV <- matrix(0, m, jimax); 856 kVU <- matrix(0, m, jimax);kVV <- matrix(0, m, jimax) 857 858 for(i in 1:m){ 859 860 ## Compute mean squared distances of cluster frailty estimators 861 ci[i] <- sigma2hat * wi[i] * (1 + sigma2dhat * si[i]) 862 cid[i] <- sigma2dhat * wi[i] * (1 + sigma2hat * pi[i]) 863 864 for(j in 1:Ji[i]){ 865 866 ## Compute useful covariance terms for the bias correction 867 kUU[i, j] <- sigma2hat * wi[i] * ((1 + sigma2dhat * si[i]) * 868 (sigma2hat * pi[i] + nu2hat * pij[i, j] + thetahat * qij[i, j]) - 869 sigma2dhat * qi[i] * (sigma2hat * qi[i] + nu2hat * qij[i, j] + 870 thetahat * sij[i, j])) 871 kUV[i, j] <- sigma2hat * wi[i] * ((1 + sigma2dhat * si[i]) * 872 (thetahat * pij[i, j] + nu2dhat * qij[i, j]) - sigma2dhat * 873 qi[i] * (thetahat * qij[i, j] + nu2dhat * sij[i, j] - 1)) 874 kVU[i, j] <- sigma2dhat * wi[i] * ((1 + sigma2hat * pi[i]) * 875 (thetahat * sij[i, j] + nu2hat * qij[i, j]) - sigma2hat * 876 qi[i] * (thetahat * qij[i, j] + nu2hat * pij[i, j] - 1)) 877 kVV[i, j] <- sigma2dhat * wi[i] * ((1 + sigma2hat * pi[i]) * 878 (sigma2dhat * si[i] + nu2dhat * sij[i, j] + thetahat * qij[i, j]) - 879 sigma2hat * qi[i] * (sigma2dhat * qi[i] + nu2dhat * qij[i, j] + 880 thetahat * pij[i, j])) 881 882 ## Compute mean squared distances of individual frailty estimators 883 cij[i, j] <- (nu2hat * pij[i, j] + thetahat * qij[i, j] - 1) * 884 (kUU[i, j] - sigma2hat - nu2hat) + (nu2hat * qij[i, j] + thetahat * 885 sij[i, j]) * (kVU[i, j] - thetahat) 886 cijd[i, j] <- (nu2dhat * sij[i, j] + thetahat * qij[i, j] - 1) * 887 (kVV[i, j] - sigma2dhat - nu2dhat) + (nu2dhat * qij[i, j] + thetahat * 888 pij[i, j]) * (kUV[i, j] - thetahat) 889 cijprime[i, j] <- thetahat - (1 - nu2hat * pij[i, j] - thetahat * 890 qij[i, j]) * kUV[i, j] + (nu2hat * qij[i, j] + thetahat * 891 sij[i, j]) * kVV[i, j] - (sigma2dhat + nu2dhat) * 892 (nu2hat * qij[i, j] + thetahat * sij[i, j]) - thetahat * 893 (nu2hat * pij[i, j] + thetahat * qij[i, j]) 894 895 } 896 } 897 898 ## Compute estimators of dispersion parameters 899 sigma2hatnew <- corrval * 1 / m * sum((Uihat - 1)^2 + ci) 900 sigma2dhatnew <- corrval * 1 / m * sum((Vihat - 1)^2 + cid) 901 nu2hatnew <- 0;nu2dhatnew <- 0; thetahatnew <- 0 902 for(i in 1:m){ 903 nu2hattemp <- 0;nu2dhattemp <- 0;thetahattemp <- 0; 904 for(j in 1:Ji[i]){ 905 nu2hattemp <- nu2hattemp + (Uijhat[Jicum[i] + j] - Uihat[i])^2 + 906 ci[i] + cij[i, j] - 2 * (sigma2hat - kUU[i, j]) 907 nu2dhattemp <- nu2dhattemp + (Vijhat[Jicum[i] + j] - Vihat[i])^2 + 908 cid[i] + cijd[i, j] - 2 * (sigma2dhat - kVV[i, j]) 909 thetahattemp <- thetahattemp + (Uijhat[Jicum[i] + j] - Uihat[i]) * 910 (Vijhat[Jicum[i] + j] - Vihat[i]) + cijprime[i, j] + kUV[i, j] + 911 kVU[i, j] - sigma2hat * sigma2dhat * wi[i] * qi[i] 912 } 913 nu2hatnew <- nu2hatnew + nu2hattemp / Ji[i]; 914 nu2dhatnew <- nu2dhatnew + nu2dhattemp / Ji[i]; 915 thetahatnew <- thetahatnew + thetahattemp / Ji[i]; 916 } 917 # Ad - hoc correction adjustment 918 nu2hatnew <- corrval * nu2hatnew / m 919 nu2dhatnew <- corrval * nu2dhatnew / m 920 thetahatnew <- corrval * thetahatnew / m 921 # Format for output 922 return(list(sigma2hat = sigma2hatnew, 923 sigma2dhat = sigma2dhatnew, 924 nu2hat = nu2hatnew, 925 nu2dhat = nu2dhatnew, 926 thetahat = thetahatnew)) 927 }