| ######################################################################### |
| # R (ver 2.14.1) program to compute QLS estimates for the mixture model # |
| ######################################################################### |
| # Function to estimate the correlation matrix |
| # between the repeated measurements |
| correlation.est <- function(residuals, tol=1e-10) |
| { |
| t <- ncol(residuals) |
| Z <- t(residuals)%*%residuals |
| # start the decomposition algorithm with an identity matrix |
| Lambda0 <- diag(t) |
| ev <- eigen(Z) |
| Lambdak <- diag(diag(ev$vec%*%diag(sqrt(ev$val))%*%t(ev$vec))) |
| Diff <- diag(Lambdak -Lambda0) |
| while(sum(Diff∧2) > tol) |
| { |
| Lambda0 <- Lambdak |
| ev <- eigen(sqrt(Lambda0)%*% Z %*%sqrt(Lambda0)) |
| M <- ev$vec%*% diag(sqrt(ev$val)) %*%t(ev$vec) |
| Lambdak <- diag(diag(M)) |
| Diff <- diag(Lambdak - Lambda0) |
| } |
| Rtilde <- solve(sqrt(Lambdak))%*% M %*%solve(sqrt(Lambdak)) |
| Rhat <- Rtilde%*%diag(as.vector(solve(Rtilde*Rtilde)%*%rep(1,t)))%* |
| %Rtilde |
| ev <- eigen(Rhat) |
| if (any(ev$val<0)) |
| Rhat <- solve(sqrt(diag(diag(Z))))%*% Z %*% solve |
| (sqrt(diag(diag(Z)))) return(Rhat) |
| } |
| # Function to calculate the regression parameter beta. |
| regression.est <- function(x, y, t, Rhat) |
| { |
| mt <- nrow(x) |
| Sigma <- solve(kronecker(diag(mt/t), Rhat)) |
| XRinvX <- t(x)%*%Sigma%*%x |
| XRinvY <- t(x)%*%Sigma%*%y |
| betahat <- solve(XRinvX)%*%XRinvY |
| return(betahat) |
| } |
| # The main program starts here |
| d <- read.table("c:/hatch-eggs.txt", header=TRUE) |
| proportion <- d$Hatch/d$Total |
| y <- log(proportion/(1-proportion)) |
| x <- model.matrix(~Salinity+Temperature, data=d) |
| tol <- 1e-10 |
| t <- length(d$ID)/length(unique(d$ID)) |
| mt <- nrow(x) |
| m <- mt/t |
| k <- ncol(x) |
| Rhatinit <- diag(t) |
| betahat <- regression.est(x, y, t, Rhatinit) |
| residuals <- matrix(y-x%*%betahat, ncol=t, byrow=TRUE) |
| Rhat <- correlation.est(residuals) |
| betanew <- regression.est(x, y, t, Rhat) |
| while(sum((betanew-betahat)∧2)>tol) |
| { |
| betahat <- betanew |
| residuals <- matrix(y-x%*%betahat, ncol=t, byrow=TRUE) |
| Rhat <- correlation.est(residuals) |
| betanew <- regression.est(x, y, t, Rhat) |
| } |
| # Calculate the scale parameter |
| residuals <- matrix(y-x%*%betahat, ncol=t, byrow=TRUE) |
| Z <- t(residuals)%*% residuals |
| Rhat <- correlation.est(residuals) |
| scale <- sum(diag(solve(Rhat)%*%Z))/(mt-k) |
| # Calculate model based standard errors and z-scores for betas |
| Sigma <- solve(kronecker(diag(m), Rhat)) |
| Covbeta <- scale*solve(t(x)%*%Sigma%*%x) |
| stderrbeta <- sqrt(diag(Covbeta)) |
| zstat <- betanew/stderrbeta |
| # Prepare and print the output |
| output <- cbind(betanew, stderrbeta, zstat, 2*(1-pnorm(abs(zstat)))) |
| colnames(output) <- c("Estimate", "Std. Error", "z value", "P-value") |
| list(scale = scale, Rhat=Rhat, beta = output) |
