# Edps/Soc 584 and Psych 594 example on how to
#
# Read in data set into a matiex
# Write functions to compute
# (i) n, xbar
# (ii) S
# (iii) R
# (iv) n, xbar, S, and R
#
#
#
# First get a sample data set
#
scores <- read.table('scores_data.txt', header=TRUE)
attach(scores)
names(scores) # look at the variable names
head(scores) # see the first 6 lines of data
tail(scores) # see the last 6 lines of data
X <- as.matrix(scores[,2:5]) # put the variables for the test scores
# into matrix X
p <- ncol(X) # Check to make sure things worked OK
p
n <- nrow(X)
head(X)
tail(X)
#
# Define function to compute mean vector that prints out the
# n = sample size
# xbar = column vector with means in it
#
meanvector <- function(X){
n <- nrow(X)
one <- array(1, dim=c(n))
xbar <- t(X) %*% one/n
list(n=n, xbar=xbar)
}
#
# Examples of running the function:
#
# (a) Print n and xbar to console
meanvector(X)
# (b) Print only xbar to console
meanvector(X)$xbar
# (c) Print only n to console
meanvector(X)$n
# (d) Set xbartest equal to the mean vector
xbartest <-meanvector(X)$xbar
#
# Function that computes S = sample covariance matrix
#
Smatrix <- function(X) {
n <- nrow(X)
one <- array(1, dim=c(n))
xbar <- t(X) %*% one/n
W <- t(X - one %*% t(xbar)) %*% (X - one %*% t(xbar))
S <- W / (n-1)
list(S=S)
}
#
# Function that computes R = sample covariance matrix
#
# Note that this one calls the function Smatrix
#
Rmatrix <- function(X){
S <- Smatrix(X)$S
Dsqrt <- diag(diag(sqrt(S)) )
invD <- solve(Dsqrt)
R <- invD %*% S %*% invD
list (R=R)
}
#
# Function that computes
# n = sample size
# xbar= column vector of means
# S = sample covariance matrix
# R = sampe correlation matirx
#
samplestats <- function(X) {
n <- nrow(X)
one <- array(1, dim=c(n))
xbar <- t(X) %*% one/n
W <- t(X - one %*% t(xbar)) %*% (X - one %*% t(xbar))
S <- W / (n-1)
Dsqrt <- diag(diag(sqrt(S)) )
invD <- solve(Dsqrt)
R <- invD %*% S %*% invD
list (n = n, xbar=xbar, S = S, R=R )
}
#
# Read male and female data into different matrices
# and compute mean vectors and the difference between them
#
male <- as.matrix(scores[1:32,2:5])
female <- as.matrix(scores[33:64,2:5])
maleall <- samplestats(male)
mbar <- maleall$xbar
femaleall <- samplestats(female)
fbar <- femaleall$xbar
difference <- mbar - fbar