Maximum likelihood and other parameter estimations

First, a basic definition: A parameter is an unknown, fixed value that describes a characteristic. For example, a mean is a characteristic that describes the average over a population. The true mean is usually not known, but rather estimated.

When we fit a model to our data, we get parameters such as the regression coefficients (β's). In spatial stats, we use a (semi)variogram function to estimate the parameters range, sill, and nugget effect. An introduction to the semivariogram and its parameters may be found here.

A few common methods of parameter estimation used in spatial stats are the least squares (OLS or more commonly WLS) and the likelihood based methods (maximum likelihood MLE or restricted maximum likelihood REML).

Least squares methods fit a model by minimizing the distance between the observed data and the best fit line. Likelihood based methods use the observed data to estimate the population parameters using established distributions. When you code a likelihood estimation, you will input parameters and an underlying distribution. For example, for a spatial stats dataset, you would first investigate the semivariogram to estimate the nugget effect parameter and the distribution (i.e. exponential or linear), and then model the data via MLE.

In R, the likfit command in the package geoR models likelihood based methods. In the same package, variofit models least squares.

Interpretting linear models in R

If you're new to R and stats, check out this awesome post over at the yhatq blog. It walks you through everything, from the code to the analysis, in simple, straight-forward language with code output.

On residuals:

If our residuals are normally distributed, this indicates the mean of the difference between our predictions and the actual values is close to 0 (good) and that when we miss, we're missing both short and long of the actual value, and the likelihood of a miss being far from the actual value gets smaller as the distance from the actual value gets larger.

On variable p-values:

Probability the variable is NOT relevant. You want this number to be as small as possible. If the number is really small, R will display it in scientific notation. In or example 2e-16 means that the odds that parent is meaningless is about ¹⁄_{5000000000000000}

On R-squared:

Metric for evaluating the goodness of fit of your model. Higher is better with 1 being the best. Corresponds with the amount of variability in what you're predicting that is explained by the model.

On the F-test and resulting F-stat:

This takes the parameters of our model (in our case we only have 1) and compares it to a model that has fewer parameters (sic). In theory the model with more parameters should fit better. If the model with more parameters (your model) doesn't perform better than the model with fewer parameters, the F-test will have a high p-value (probability NOT significant boost). If the model with more parameters is better than the model with fewer parameters, you will have a lower p-value.

All quoted text from post "Fitting & Interpreting Linear Models in R" by yhat, published May 18, 2013 at http://blog.yhathq.com/posts/r-lm-summary.html

Fun with GLM links

I can never keep this straight: When to use link(log) vs. link(logit) in generalized linear models:

Use link(log) when you are after relative risks.
Use link(logit) when you are after odds ratios.

Making notes in your code: STATA vs R

I'm used to coding in R at this point, so it took me a second to figure this out:

# This is a note in R code

vs.

* This is a note in STATA code

Other small things that are hard to get use to: in R, if your code spills into the next line, you just keep typing (I usually indent just to keep it looking clean); but in STATA, you have to include "///" at the end of an unfinished line of code, to tell it to continue reading on the next line.

STATA seems to stop running code when it encounters a problem (?) vs. R which keeps going, throwing up error messages when appropriate but following thru to the end.

Simple vector analysis in R

First, create a vector:

ClassGrades<- c(94, 95, 91)
#defines a vector called "ClassGrades" by concatenating individual grades on the exam

Then, run your simple analysis

mean(ClassGrades)
#returns the mean on the exam

hist(GlassGrades)
#produces a histogram of the exam grades

Another R tutorial

http://pages.pomona.edu/~jsh04747/courses/R.pdf

Graphing in R

Frank McCown outlines a quick and simple intro exercise for graphing in R:

http://www.harding.edu/fmccown/r/

A few notes to keep in mind as you run through the exercise:

• If you need a refresher on creating tab delimited files, ehow has you covered:
http://www.ehow.com/how_5079116_create-tab-delimited-files.html
• If you save your tabbed file with .tab per the ehow instructions, make sure to edit Frank's code:
autos_data <- read.table("C:/R/autos.dat", header=T, sep="\t")
#becomes
autos_data <- read.table("C:/R/autos.tab", header=T, sep="\t"
• Make sure you're using forward slashes when you set the file path! This post discusses that pesky character string error.

How the big boys code

Don't just code into the R console, which is the default window that opens when you first open R. Instead, start a new R editor by

File -> New script

Here you can write out your code as well as notes.

To make notes, put # in front of it:

this is code

#this is note

To run a line of code, you can either copy and paste from the console or highlight it and hit control and R.

Setting your working directory

You tried to set your working directory using the location listed in the properties tab of your file of interest by coding something like:

> setwd("C:\Users\Me\Documents\Rfolder\lab1")

But you get the following error:

Error: '\u' used without hex digits in character string starting "C:\u"

Q: What to do???

A: Change all your backslashes \ to forward slashes /

> setwd("C:/Users/Me/Documents/Rfolder/lab1")

> getwd()

[1] "C:/Users/Me/Documents/Rfolder/lab1"

Or you could just shortcut it:

File -> Change dir...

But that may be too easy

Chi squared test of independence in R

1. Create or import table into R
EX: I used the same data from this post, but I had edited the Excel text file to reflect the column and row titles and renamed it "Table2forRhead.txt".

> getwd()
[1] "C:/Documents and Settings/me/My Documents"
>
> setwd("E:")
>
> mydata = read.table("Table2forRhead.txt")
> mydata
Educ Nav
Male 336 355
Female 955 947

2. Run chi square test:
> chisq.test(mydata)

Pearson's Chi-squared test with Yates' continuity correction

data: mydata
X-squared = 0.448, df = 1, p-value = 0.5033

3. Report p-value to test hypothesis of independence
EX: In this problem, p-value = 0.5033, which is > 0.05 (the significance level), therefore we do not reject the null hypothesis that sex is independent of the intervention group of the subjects

BUT... what if we want to test null that invention group is independent of sex (reverse the chi squared test)?

1. Create new txt in Excel, reversing the columns and rows
EX: "Table2forRheadReverse.txt"
2. Repeat above
> getwd()
[1] "E:/"
> mydata2 = read.table("Table2forRheadReverse.txt")
> mydata2
Male Female
Educ 336 955
Nav 355 947
> chisq.test(mydata2)

Pearson's Chi-squared test with Yates' continuity correction

data: mydata2
X-squared = 0.448, df = 1, p-value = 0.5033

THEY ARE THE SAME! So p-value can be interpreted as failing to reject the null hypothesis that sex is independent of the intervention group AND failing to reject the null that intervention group is independent of sex.

Messed up R code?

... And now it won't let you start a new line (aka hitting "enter" starts a continued line, denoted with "+" instead of "<")?

Just hit the red "Stop" button located on the top.

Importing Excel data into R

First save Excel into text file (see earlier post).
EX: I created a 2x2 table with no headers and saved it to E drive (dat path E:\Table2forR.txt)

In R, check working directory
> getwd()

If it's not correct (in this case, not set to E drive), reset it:
> setwd("E:")

Set R to read your table and rename it (renamed from "Table2forR" to "mydata" here)
> mydata = read.table("Table2forR.txt")

Check it:
> mydata

Should get output:
V1 V2
1 336 355
2 955 947

Creating a 2x2 table in R

EX: Want to create a 2x2 table with the following data: Column A value (336 and 955) and Column B value (355 and 947)

Note from R-tutorial: "When we construct a matrix directly with data elements, the matrix content is filled along the column orientation by default"

Code:
> B = matrix (c(336, 955, 355, 947), nrow=2, ncol=2)

Check:
> B

Should see output:
[,1] [,2]
[1,] 336 355
[2,] 955 947

R vocab

From R Tutuorial

A note on data terminology:
"A vector is a sequence of data elements of the same basic type. Members in a vector are officially called components. Nevertheless, we will just call them members in this site."

R links

chi squared test page http://statland.org/R/Rchisq.htm
Note coding differences for goodness of fit vs test of independence

cran's "R data Import/Export" http://cran.r-project.org/doc/manuals/R-data.html#Spreadsheet_002dlike-data

wikiversity's "How to Use R" http://en.wikiversity.org/wiki/How_to_use_R#Tutorials

cran's setting a working directory http://127.0.0.1:23630/library/base/html/getwd.html

R tutorial http://www.r-tutor.com/

Saving Excel as Tab delimited for R

Excel:
Save As -> Other formats (may have to scroll down to see this option) -> Save as type: Text (tab delimited) (scroll down to bottom for this option)