###################################################################################### # Edps/Psych/Stat 587 # Fall 2016 # Carolyn J. Anderson # # Last revised # # Function to do same constrasts as SAS (used standard multivariate concepts # regarding linear combination of normaly distributed random variables) # # Ho: L*Gamma = 0 Ha: L*Gamma <> 0 # # Where L is an (c X p) matirx of c contrast on the p fixed effects # Gamma is (p X 1) vector of fixed effects # # # Input: model # matrix (vector) with number of columns=number of fixed effects # and elements are constrast coefficients (sum to 0) # # Output: # F statistics # Guess as df # p-value for F # X2 (chi-square test statistic) # df # p-value for X2 # # Requirements: Model should be fit by lmer with lmerTest (the latter is used # as a guess at den df for F). # ##################################################################################### contrast <- function(model,L) { gamma <- fixef(model) nfixed <- length(gamma) # number of fixed effects ftable <- summary(model)[10] # use this for df (satterthwaite) ftable <- matrix(unlist(ftable),nrow=nfixed) # so I can put out df for fixed effects cov.gamma <- vcov(model) # covariance matrix of parameter estimates # check that number of columms of L = number of fixed effects nL <- ncol(L) if (nfixed != nL) { print('The number of columns of L must equal number of fixed effects.') return( ) } # check for linearly dependent rows of L ck <- L %*% t(L) det.ck <- det(ck) if (det.ck <= .0000001) { print('Check for linearly dependent rows of L') return() } estimate<- L %*% gamma # value of constrast(s) # F test statistic F <- as.numeric((t(estimate) %*% solve(L %*% cov.gamma %*% t(L)) %*% estimate)/nrow(L)) # Wald test statistic X2 <- as.numeric((t(estimate) %*% solve(L %*% cov.gamma %*% t(L)) %*% estimate)) which.df <- which(L !=0, arr.ind=TRUE) # Find elements of L that are non-zero ddf <- ftable[which.df[1,2],3] # I don't know what these should be ndf <- nrow(L) # these are fine pvalue <- pf(F,ndf,ddf,lower.tail=FALSE) # p-value of F pchi <- pchisq(X2,ndf,lower.tail=FALSE) # p-value for Wald # output results in nice table/format result <- c(F, ndf, ddf, pvalue, X2, ndf, pchi) names(result) <- c('F', 'num df', 'den df guess', 'p-value', 'X2', 'df', 'p-value') return(result) }