# How to compute an accuracy measure based on RMSE? Is my large dataset normally distributed?

I have several datasets on the order of thousands of points. The values in each dataset are X,Y,Z referring to a coordinate in space. The Z-value represents a difference in elevation at coordinate pair (x,y).

Typically in my field of GIS, elevation error is referenced in RMSE by subtracting the ground-truth point to a measure point (LiDAR data point). Usually a minimum of 20 ground-truthing check points are used. Using this RMSE value, according to NDEP (National Digital Elevation Guidelines) and FEMA guidelines, a measure of accuracy can be computed: Accuracy = 1.96*RMSE.

This Accuracy is stated as: "The fundamental vertical accuracy is the value by which vertical accuracy can be equitably assessed and compared among datasets. Fundamental accuracy is calculated at the 95-percent confidence level as a function of vertical RMSE."

I understand that 95% of the area under a normal distribution curve lies within 1.96*std.deviation, however that does not relate to RMSE.

Generally I am asking this question: Using RMSE computed from 2-datasets, how can I relate RMSE to some sort of accuracy (i.e. 95-percent of my data points are within +/- X cm)? Also, how can I determine if my dataset is normally distributed using a test that works well with such a large dataset? What is "good enough" for a normal distribution? Should p<0.05 for all tests, or should it match the shape of a normal distribution?

I found some very good information on this topic in the following paper:

http://paulzandbergen.com/PUBLICATIONS_files/Zandbergen_TGIS_2008.pdf

-
Watch out! Your use of ks.test is incorrect. According to the help page, you need to use 'pnorm' instead of 'dnorm'. Moreover, setting the comparison distribution's parameters to the mean and SD of the sample itself will substantially inflate the p-value: "If a single-sample test is used, the parameters specified in ... must be pre-specified and not estimated from the data." –  whuber Jan 31 '12 at 20:29
Well, actually, that formula will not give you a confidence interval: it will be far too large for that. It really is a crude (but standard) way to estimate a tolerance interval, which is the middle 95% of the whole population of differences. There are good reasons to suppose the differences will not have a normal distribution: larger absolute differences tend to be associated with larger topographic slopes. Assuming your 4000 points are a random sample of those differences, why don't you just report their 2.5 and 97.5 percentiles? –  whuber Jan 31 '12 at 20:54
Your data form a statistical sample of the elevations that could be measured. When you talk about "accuracy" you are making claims about how closely your DEMs represent the entire population of elevations. In your case, it is impossible to assess accuracy by comparing datasets: you have to "field-truth" your data. Thus, the guidelines are really talking about relative agreement of two datasets. Finally, their use of "confidence level" is mistaken, as I explained earlier. I accept you have to work within the framework of awful guidance like this, but you deserve to know what's correct. –  whuber Jan 31 '12 at 21:46
That's starting to sound like a useful question for you. Because you haven't yet received any answers, why don't you just completely edit the current question to incorporate the information you have disclosed in these comments? I would suggest broadening it somewhat: after quoting the guidelines (to show what kind of methods are usually employed in your field), you might ask quite generally how to use the distribution of the ordered pairs of differences in elevations to assess accuracy (assuming one of the data sets is the reference). –  whuber Jan 31 '12 at 23:22
All: Updated my main post and question to reflect the changes from the comments. –  Matthew Bilskie Feb 1 '12 at 0:21

This is a partial answer for:

"Is my large dataset normally distributed"?

A good start would be to observe your "z" laser values histograms.

You did not specify the size of your files, but I do that with ".las" files in [R] using the readLAS function created by Michael Sumner.

Here is a reproducible example:

Download ".las" file available here (it is the first example -124 Mb).

Run the readLAS function in R.

#Read the large LiDAR (.las) dataset with readLAS function.
nfile = paste("C:\\...\\lda_4800K_data.las",sep="") #set the file path correctly
lasframe = data.frame(lasfiles)
z = lasframe\$z

[1] 1227.979 1228.462 1228.442 1228.320 1228.209 1228.002
length(z)
[1] 4826890 #very large dataset

#Histogram
library(ggplot2)
ggplot(lasframe,aes(x=z)) +
geom_histogram(binwidth=25,aes(y=..density..), fill="white", color="black") +
ylab("Density") + xlab("Elevation (feet)") +
theme_bw() +
coord_flip()


If you want to analyse small samples over the dataset you can use Shapiro-Wilk Normality Test.
In R, shapiro.test only supports 5000 observations at most.

z_sample = sample(z,40,replace=T)
shapiro.test(z_sample) #high p-value indicates the data is normal (null hypothesis)

Shapiro-Wilk normality test

data:  z_sample
W = 0.9323, p-value = 0.01911 #not normal


However, normality tests on large datasets are not so useful as it is explained in this answer provided by Greg Snow.

"On the other hand, with really large datasets the central limit theorem comes into play and for common analyses (regression, t-tests, ...) you really don't care if the population is normally distributed or not."

#qq-plot (quantiles from empirical distribution - quantiles from theoretical distribution)