SEEG-Oxford / movement

R package containing useful functions for the analysis of movement data in disease modelling and mapping
20 stars 16 forks source link

Create a new example for the movement() function documentation #69

Closed KathrinTessella closed 8 years ago

KathrinTessella commented 8 years ago

The old example under the movement() function is not useful anymore. Add a new example which applies the movement() function directly.

KathrinTessella commented 8 years ago

A first attempt for a new example is the following:

data(kenya) kenya10 <- raster::aggregate(kenya, 10, sum) flux <- gravity() predictedMovement <- predict(flux, kenya10) movementMatrix <- predictedMovement$movement_matrix class(movementMatrix) <- c('matrix', 'movementmatrix') stopifnot(is.movementmatrix(movementMatrix))

net <- get.network(kenya10, min = 50000) locationData <- data.frame(location = net$locations, population = net$population, x = net$coordinate[,1], y = net$coordinate[,2]) class(locationData) <- c('data.frame', 'locationdataframe') stopifnot(is.locationdataframe(locationData))

s <- movement(movementMatrix ~ locationData, flux)

Note: need to adjust the call for extracting the arguments from the formula to

movementmatrix <- get(as.character(formula[[2]]), mode = 'numeric') locationdataframe <- get(as.character(formula[[3]]), mode = 'list')

to be fully working. However, this example causes an error in the 'optim' function. Error in optim(predictionModel$flux_model$params, fittingwrapper, method = "BFGS", : initial value in 'vmmin' is not finite

This potentially could be address with a different set of initial values or the implementation of the transformations for the parameters to be fitted. I will come back to this once the task in https://github.com/SEEG-Oxford/movement/issues/46 has been implemented.

KathrinTessella commented 8 years ago

@goldingn : I run various predictions of different theoretical flux models with various values of parameter (e.g. change theta for the original radiation flux model), but none of the resulting movement matrices will work with the locationData I extract from the Kenya data

data(kenya)
kenya10 <- raster::aggregate(kenya, 10, sum)
predictedMovement  <- predict(originalRadiation(theta = 0.1), kenya10, symmetric = TRUE)
movementMatrix <- predictedMovement$movement_matrix
class(movementMatrix) <- c('matrix', 'movement_matrix')

net <- getNetwork(kenya10, min = 50000)
locationData <- data.frame(location = net$locations, population = net$population, x = net$coordinate[,1], y = net$coordinate[,2])
class(locationData) <- c('data.frame', 'location_dataframe')

s <- movement(movementMatrix ~ locationData, originalRadiation(theta = 0.1))

In the beginning of the project, Moritz provided me with the data set for the movement matrix 'twitter_west_africa_1.csv' which (once I extract the first 138 rows & columns) works fine using the population data of the Kenya example. But not sure if these data are public available to be used here?

goldingn commented 8 years ago

Alas those data are not publicly available or shareable. They also cover the wrong countries anyway.

I suspect these data don't line up because the network is generated with a different population threshold in the two models, so they are of differing sizes.

This works for me:

# get location data
data(kenya)
kenya10 <- raster::aggregate(kenya, 10, sum)
net <- getNetwork(kenya10, min = 50000)
locationData <- data.frame(location = net$locations, population = net$population, x = net$coordinate[,1], y = net$coordinate[,2])
class(locationData) <- c('data.frame', 'location_dataframe')

# simulate movements (note the values of movementmatrix must be integer)
predictedMovement  <- predict(originalRadiation(theta = 0.1), locationData, symmetric = TRUE)
movementMatrix <- round(predictedMovement$movement_matrix)
class(movementMatrix) <- c('matrix', 'movement_matrix')

# fit a new model to these data
s <- movement(movementMatrix ~ locationData, radiationWithSelection(theta = 0.5))
Training complete.
goldingn commented 8 years ago

Note the round() bit I added - I suspected this was breaking because the movements were non integers (it's not possible to observe half movements)

Could we add some coercion into as.movement_matrix to make sure these are integer-ish (R can handle rounded numerics in this case, and the integer class is rarely used in applications like this), perhaps issuing a warning if rounding was necessary (using e.g. all.equal(x, round(x)))

I'm not sure if this is on the to-do list, but adding the class membership shouldn't be something the user does, but should be done internally via the as.* methods too...

KathrinTessella commented 8 years ago

Thanks! Yes - suspect the rounding was missing!

I also don't like that the

predictedMovement  <- predict(originalRadiation(theta = 0.1), locationData, symmetric = TRUE)

does not return a movement_matrix of the correct class type. I will see if that is easy to adjust as it can be confusing for the user to have to assign the correct class to the object later on.

KathrinTessella commented 8 years ago

All tasks are completed.