4.10 The Cox model fitting
Hening Cui
4/4/2023
Preparing Data
loading library
library(Seurat)## Attaching SeuratObjectlibrary(spacexr)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2
## ──## ✔ ggplot2 3.4.2 ✔ purrr 1.0.1
## ✔ tibble 3.2.1 ✔ dplyr 1.1.1
## ✔ tidyr 1.3.0 ✔ stringr 1.5.0
## ✔ readr 2.1.4 ✔ forcats 1.0.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(spatstat)## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 3.0-6
## Loading required package: spatstat.random
## spatstat.random 3.1-3
## Loading required package: spatstat.explore
## Loading required package: nlme
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapse
##
## spatstat.explore 3.0-6
## Loading required package: spatstat.model
## Loading required package: rpart
## spatstat.model 3.2-1
## Loading required package: spatstat.linnet
## spatstat.linnet 3.0-6
##
## spatstat 3.0-3
## For an introduction to spatstat, type 'beginner'library(nlme)loading dataset
load(file= "/Users/cuihening/Desktop/shikun/seurat/robjs/chicken_visium.4.Robj")TM <- FetchData(object=chicken_visium,vars=c('TMSB4X','ident', 'orig.ident'))
tmsb4x=cbind(cell = rownames(TM), TM)Try first with D4 data
Prepare data
D4 coordinate
coords_d4 <- GetTissueCoordinates(chicken_visium, image= 'slice1')
colnames(coords_d4) <- c("x", "y")
coords_d4 <- cbind(cell = rownames(coords_d4),coords_d4)D4 dataset
D4 = tmsb4x %>%
filter(orig.ident=='D4') %>%
merge(coords_d4, by="cell") EDA
Build D4 data ppp project
X <- D4$x
Y <- D4$y
xrange <- range(X, na.rm=T)
yrange <- range(Y, na.rm=T)
D4_tm = ppp(X,Y,xrange,yrange,marks = D4$TMSB4X)
par(mfrow=c(2,2))
plot(unmark(D4_tm),main="D4 only point pattern")
plot(D4_tm,main="TMSB4X marker spatial plot")
hist(marks(D4_tm), main="D4-Distribution of TMSB4X exprssion")
D4_ty = ppp(X,Y,xrange,yrange,marks = D4$ident)
plot(D4_ty, main="D4 cell type distribution")
Different cell type density plot
D4cell_type= split(D4_ty)
plot(D4cell_type[0:9], main = "D4 cell location by type")
plot(density(D4cell_type[0:9]), main = "D4 Densities of cell type")
For type 2 cell fitting model
Type 2 cell data
D4_ty2 = D4cell_type$'2'Fitting LMM model for the TMSB4
D4_ty2df = D4 %>% filter(ident==2)
LMM <- lme (TMSB4X ~ x+y, random = ~1 | cell, data = D4_ty2df, method='REML')
summary(LMM)## Linear mixed-effects model fit by REML
## Data: D4_ty2df
## AIC BIC logLik
## -139.9427 -134.72 74.97133
##
## Random effects:
## Formula: ~1 | cell
## (Intercept) Residual
## StdDev: 0.003814857 0.001430571
##
## Fixed effects: TMSB4X ~ x + y
## Value Std.Error DF t-value p-value
## (Intercept) 0.04710706 0.012752401 21 3.693976 0.0013
## x -0.00007802 0.000051935 21 -1.502276 0.1479
## y 0.00011638 0.000060047 21 1.938094 0.0662
## Correlation:
## (Intr) x
## x -0.436
## y -0.111 -0.845
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -0.62281873 -0.20537346 0.04306407 0.14733176 0.57358082
##
## Number of Observations: 24
## Number of Groups: 24Fitting Cox spatial model
cf=coef(LMM)
TMfunction <- function(x,y){ cf$`(Intercept)` -0.00007802*x + 0.00011638*y}
D4_fit = kppm(unmark(D4_ty2) ~TMSB4X, "LGCP", data = D4_ty2, covariates=list(TMSB4X=TMfunction))## Warning in cf$`(Intercept)` - 7.802e-05 * x: longer object length is not a
## multiple of shorter object length
## Warning in cf$`(Intercept)` - 7.802e-05 * x: longer object length is not a
## multiple of shorter object lengthsummary(D4_fit)## Inhomogeneous Cox point process model
## Fitted to point pattern dataset 'unmark(D4_ty2)'
## Fitted by minimum contrast
## Summary statistic: inhomogeneous K-function
## Minimum contrast fit (object of class "minconfit")
## Model: Log-Gaussian Cox process
## Covariance model: exponential
## Fitted by matching theoretical K function to unmark(D4_ty2)
##
## Internal parameters fitted by minimum contrast ($par):
## sigma2 alpha
## 4.100185 91.651261
##
## Fitted covariance parameters:
## var scale
## 4.100185 91.651261
## Fitted mean of log of random intensity: [pixel image]
##
## Converged successfully after 83 function evaluations
##
## Starting values of parameters:
## sigma2 alpha
## 1.00000 17.08318
## Domain of integration: [ 0 , 79.67 ]
## Exponents: p= 2, q= 0.25
##
## ----------- TREND -----
## Point process model
## Fitted to data: X
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = X, trend = trend, rename.intercept = FALSE, covariates = covariates,
## covfunargs = covfunargs, use.gam = use.gam, forcefit = TRUE,
## improve.type = ppm.improve.type, improve.args = ppm.improve.args,
## nd = nd, eps = eps)
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Planar point pattern: 24 points
## Average intensity 0.000189 points per square unit
## Window: rectangle = [83.6455, 482.0155] x [66.8408, 385.5368] units
## (398.4 x 318.7 units)
## Window area = 126959 square units
##
## Dummy quadrature points:
## 32 x 32 grid of dummy points, plus 4 corner points
## dummy spacing: 12.449062 x 9.959249 units
##
## Original dummy parameters: =
## Planar point pattern: 1028 points
## Average intensity 0.0081 points per square unit
## Window: rectangle = [83.6455, 482.0155] x [66.8408, 385.5368] units
## (398.4 x 318.7 units)
## Window area = 126959 square units
## Quadrature weights:
## (counting weights based on 32 x 32 array of rectangular tiles)
## All weights:
## range: [31, 124] total: 127000
## Weights on data points:
## range: [31, 62] total: 940
## Weights on dummy points:
## range: [31, 124] total: 126000
## --------------------------------------------------------------------------------
## FITTED :
##
## Nonstationary Poisson process
##
## ---- Intensity: ----
##
## Log intensity: ~TMSB4X
## Model depends on external covariate 'TMSB4X'
## Covariates provided:
## TMSB4X: function(x, y)
##
## Fitted trend coefficients:
## (Intercept) TMSB4X
## -8.996264 8.095000
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -8.996264 0.7768394 -10.51884 -7.473687 *** -11.5805966
## TMSB4X 8.095000 14.1250853 -19.58966 35.779659 0.5730939
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept) TMSB4X
## -8.996264 8.095000
##
## Fitted exp(theta):
## (Intercept) TMSB4X
## 1.238717e-04 3.278038e+03
##
## ----------- COX -----------
## Model: log-Gaussian Cox process
##
## Covariance model: exponential
## Fitted covariance parameters:
## var scale
## 4.100185 91.651261
## Fitted mean of log of random intensity: [pixel image]
##
## Final standard error and CI
## (allowing for correlation of Cox process):
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -8.996264 4.65421 -18.11835 0.1258199 -1.93293042
## TMSB4X 8.095000 84.79243 -158.09511 174.2851136 0.09546843plot(D4_fit )## Warning in cf$`(Intercept)` - 7.802e-05 * x: longer object length is not a
## multiple of shorter object length## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
plot(D4_fit )## Warning in cf$`(Intercept)` - 7.802e-05 * x: longer object length is not a
## multiple of shorter object length## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
Pseudotime analysis
Take type2 as an example.
Type 2 all time
Type2 = tmsb4x %>%
filter(ident==2) %>%
rename(time = orig.ident)The function to get coordinate
getcoor=function(slice, time){
coords <- GetTissueCoordinates(chicken_visium, image= slice)
colnames(coords) <- c("x", "y")
coords <- cbind(cell = rownames(coords),coords)
coords$time = time
return(coords)
}All the time point coordinates
coords_d4 = getcoor('slice1','D4')
coords_d10 = getcoor('slice1_D10-C1','D10')
coords_d7 = getcoor('slice1_D7-B1','D7')
coords_d14 = getcoor('slice1_D14-D1','D14')Integrate to the type 2 dataset
ty2_coor= Type2 %>%
merge(coords_d4,by="cell") %>%
rbind(merge(Type2, coords_d7,by="cell")) %>%
rbind(merge(Type2, coords_d10,by="cell")) %>%
rbind(merge(Type2, coords_d14,by="cell")) %>%
mutate(time = factor(time.x))Build the ppp object
buildppp=function(dataset,marker){
X <- dataset$x
Y <- dataset$y
xrange <- range(X, na.rm=T)
yrange <- range(Y, na.rm=T)
ppppro = ppp(X,Y,xrange,yrange,marks = dataset[marker])
return(ppppro)
}EDA
type2ppp = buildppp(ty2_coor, 'TMSB4X')
hist(marks(type2ppp), main="Type 2Distribution of TMSB4X exprssion")
type2_alltime = buildppp(ty2_coor, 'time')
ty2_time= split(type2_alltime)
plot(ty2_time[0:4], main = "Type 2 cell location by type")
plot(density(ty2_time[0:4]), main = "Type 2 cell Densities across time")
Model fit
fitcox = function(data, tim, ppplist){
dataset= data %>% filter(time == tim)
LMM = lme (TMSB4X ~ x+y, random = ~1 | cell, data = dataset, method='REML')
ppppro = ppplist[[tim]]
cf=coef(LMM)
TMfunction <- function(x,y){ cf$`(Intercept)` -cf$x*x + cf$y*y}
cox_fit = kppm(unmark(ppppro) ~TMSB4X, "LGCP", data = ppppro, covariates=list(TMSB4X=TMfunction))
return(cox_fit)}D4 = fitcox(ty2_coor,'D4',ty2_time)## Warning in cf$x * x: longer object length is not a multiple of shorter object
## length## Warning in cf$`(Intercept)` - cf$x * x: longer object length is not a multiple
## of shorter object length## Warning in cf$y * y: longer object length is not a multiple of shorter object
## length## Warning in cf$x * x: longer object length is not a multiple of shorter object
## length## Warning in cf$`(Intercept)` - cf$x * x: longer object length is not a multiple
## of shorter object length## Warning in cf$y * y: longer object length is not a multiple of shorter object
## lengthsummary(D4)## Inhomogeneous Cox point process model
## Fitted to point pattern dataset 'unmark(ppppro)'
## Fitted by minimum contrast
## Summary statistic: inhomogeneous K-function
## Minimum contrast fit (object of class "minconfit")
## Model: Log-Gaussian Cox process
## Covariance model: exponential
## Fitted by matching theoretical K function to unmark(ppppro)
##
## Internal parameters fitted by minimum contrast ($par):
## sigma2 alpha
## 3.834707 42.610325
##
## Fitted covariance parameters:
## var scale
## 3.834707 42.610325
## Fitted mean of log of random intensity: [pixel image]
##
## Converged successfully after 59 function evaluations
##
## Starting values of parameters:
## sigma2 alpha
## 1.00000 17.08318
## Domain of integration: [ 0 , 105.9 ]
## Exponents: p= 2, q= 0.25
##
## ----------- TREND -----
## Point process model
## Fitted to data: X
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = X, trend = trend, rename.intercept = FALSE, covariates = covariates,
## covfunargs = covfunargs, use.gam = use.gam, forcefit = TRUE,
## improve.type = ppm.improve.type, improve.args = ppm.improve.args,
## nd = nd, eps = eps)
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Planar point pattern: 24 points
## Average intensity 0.000124 points per square unit
## Window: rectangle = [51.0357, 506.8714] x [106.7723, 530.5386] units
## (455.8 x 423.8 units)
## Window area = 193168 square units
##
## Dummy quadrature points:
## 32 x 32 grid of dummy points, plus 4 corner points
## dummy spacing: 14.24487 x 13.24269 units
##
## Original dummy parameters: =
## Planar point pattern: 1028 points
## Average intensity 0.00532 points per square unit
## Window: rectangle = [51.0357, 506.8714] x [106.7723, 530.5386] units
## (455.8 x 423.8 units)
## Window area = 193168 square units
## Quadrature weights:
## (counting weights based on 32 x 32 array of rectangular tiles)
## All weights:
## range: [47.2, 189] total: 193000
## Weights on data points:
## range: [47.2, 94.3] total: 1590
## Weights on dummy points:
## range: [47.2, 189] total: 192000
## --------------------------------------------------------------------------------
## FITTED :
##
## Nonstationary Poisson process
##
## ---- Intensity: ----
##
## Log intensity: ~TMSB4X
## Model depends on external covariate 'TMSB4X'
## Covariates provided:
## TMSB4X: function(x, y)
##
## Fitted trend coefficients:
## (Intercept) TMSB4X
## -16.61408 65.90236
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -16.61408 1.744459 -20.03316 -13.19500 *** -9.523918
## TMSB4X 65.90236 13.953855 38.55331 93.25142 *** 4.722879
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept) TMSB4X
## -16.61408 65.90236
##
## Fitted exp(theta):
## (Intercept) TMSB4X
## 6.089708e-08 4.178618e+28
##
## ----------- COX -----------
## Model: log-Gaussian Cox process
##
## Covariance model: exponential
## Fitted covariance parameters:
## var scale
## 3.834707 42.610325
## Fitted mean of log of random intensity: [pixel image]
##
## Final standard error and CI
## (allowing for correlation of Cox process):
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -16.61408 5.994865 -28.36380 -4.864361 ** -2.771385
## TMSB4X 65.90236 52.189034 -36.38626 168.190990 1.262763D7 = fitcox(ty2_coor,'D7',ty2_time)## Warning in cf$x * x: longer object length is not a multiple of shorter object
## length## Warning in cf$`(Intercept)` - cf$x * x: longer object length is not a multiple
## of shorter object length## Warning in cf$y * y: longer object length is not a multiple of shorter object
## length## Warning in cf$x * x: longer object length is not a multiple of shorter object
## length## Warning in cf$`(Intercept)` - cf$x * x: longer object length is not a multiple
## of shorter object length## Warning in cf$y * y: longer object length is not a multiple of shorter object
## lengthsummary(D7)## Inhomogeneous Cox point process model
## Fitted to point pattern dataset 'unmark(ppppro)'
## Fitted by minimum contrast
## Summary statistic: inhomogeneous K-function
## Minimum contrast fit (object of class "minconfit")
## Model: Log-Gaussian Cox process
## Covariance model: exponential
## Fitted by matching theoretical K function to unmark(ppppro)
##
## Internal parameters fitted by minimum contrast ($par):
## sigma2 alpha
## 1.929279 61.741071
##
## Fitted covariance parameters:
## var scale
## 1.929279 61.741071
## Fitted mean of log of random intensity: [pixel image]
##
## Converged successfully after 75 function evaluations
##
## Starting values of parameters:
## sigma2 alpha
## 1.00000 15.30785
## Domain of integration: [ 0 , 105.9 ]
## Exponents: p= 2, q= 0.25
##
## ----------- TREND -----
## Point process model
## Fitted to data: X
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = X, trend = trend, rename.intercept = FALSE, covariates = covariates,
## covfunargs = covfunargs, use.gam = use.gam, forcefit = TRUE,
## improve.type = ppm.improve.type, improve.args = ppm.improve.args,
## nd = nd, eps = eps)
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Planar point pattern: 243 points
## Average intensity 0.00126 points per square unit
## Window: rectangle = [51.0357, 506.8714] x [106.7723, 530.5386] units
## (455.8 x 423.8 units)
## Window area = 193168 square units
##
## Dummy quadrature points:
## 40 x 40 grid of dummy points, plus 4 corner points
## dummy spacing: 11.39589 x 10.59416 units
##
## Original dummy parameters: =
## Planar point pattern: 1604 points
## Average intensity 0.0083 points per square unit
## Window: rectangle = [51.0357, 506.8714] x [106.7723, 530.5386] units
## (455.8 x 423.8 units)
## Window area = 193168 square units
## Quadrature weights:
## (counting weights based on 40 x 40 array of rectangular tiles)
## All weights:
## range: [30.2, 121] total: 193000
## Weights on data points:
## range: [30.2, 60.4] total: 10900
## Weights on dummy points:
## range: [30.2, 121] total: 182000
## --------------------------------------------------------------------------------
## FITTED :
##
## Nonstationary Poisson process
##
## ---- Intensity: ----
##
## Log intensity: ~TMSB4X
## Model depends on external covariate 'TMSB4X'
## Covariates provided:
## TMSB4X: function(x, y)
##
## Fitted trend coefficients:
## (Intercept) TMSB4X
## -7.699397 19.687661
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -7.699397 0.9128054 -9.488463 -5.910331 *** -8.434872
## TMSB4X 19.687661 17.5106222 -14.632528 54.007850 1.124327
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept) TMSB4X
## -7.699397 19.687661
##
## Fitted exp(theta):
## (Intercept) TMSB4X
## 4.531003e-04 3.550115e+08
##
## ----------- COX -----------
## Model: log-Gaussian Cox process
##
## Covariance model: exponential
## Fitted covariance parameters:
## var scale
## 1.929279 61.741071
## Fitted mean of log of random intensity: [pixel image]
##
## Final standard error and CI
## (allowing for correlation of Cox process):
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -7.699397 1.825011 -11.27635 -4.122442 *** -4.2188231
## TMSB4X 19.687661 34.104925 -47.15676 86.532085 0.5772674D10 = fitcox(ty2_coor,'D10',ty2_time)## Warning in cf$x * x: longer object length is not a multiple of shorter object
## length## Warning in cf$`(Intercept)` - cf$x * x: longer object length is not a multiple
## of shorter object length## Warning in cf$y * y: longer object length is not a multiple of shorter object
## length## Warning in cf$x * x: longer object length is not a multiple of shorter object
## length## Warning in cf$`(Intercept)` - cf$x * x: longer object length is not a multiple
## of shorter object length## Warning in cf$y * y: longer object length is not a multiple of shorter object
## lengthsummary(D10)## Inhomogeneous Cox point process model
## Fitted to point pattern dataset 'unmark(ppppro)'
## Fitted by minimum contrast
## Summary statistic: inhomogeneous K-function
## Minimum contrast fit (object of class "minconfit")
## Model: Log-Gaussian Cox process
## Covariance model: exponential
## Fitted by matching theoretical K function to unmark(ppppro)
##
## Internal parameters fitted by minimum contrast ($par):
## sigma2 alpha
## 2.036141 91.074313
##
## Fitted covariance parameters:
## var scale
## 2.036141 91.074313
## Fitted mean of log of random intensity: [pixel image]
##
## Converged successfully after 63 function evaluations
##
## Starting values of parameters:
## sigma2 alpha
## 1.00000 16.49134
## Domain of integration: [ 0 , 105.9 ]
## Exponents: p= 2, q= 0.25
##
## ----------- TREND -----
## Point process model
## Fitted to data: X
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = X, trend = trend, rename.intercept = FALSE, covariates = covariates,
## covfunargs = covfunargs, use.gam = use.gam, forcefit = TRUE,
## improve.type = ppm.improve.type, improve.args = ppm.improve.args,
## nd = nd, eps = eps)
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Planar point pattern: 283 points
## Average intensity 0.00147 points per square unit
## Window: rectangle = [51.0357, 506.8714] x [106.7723, 530.5386] units
## (455.8 x 423.8 units)
## Window area = 193168 square units
##
## Dummy quadrature points:
## 40 x 40 grid of dummy points, plus 4 corner points
## dummy spacing: 11.39589 x 10.59416 units
##
## Original dummy parameters: =
## Planar point pattern: 1604 points
## Average intensity 0.0083 points per square unit
## Window: rectangle = [51.0357, 506.8714] x [106.7723, 530.5386] units
## (455.8 x 423.8 units)
## Window area = 193168 square units
## Quadrature weights:
## (counting weights based on 40 x 40 array of rectangular tiles)
## All weights:
## range: [30.2, 121] total: 193000
## Weights on data points:
## range: [30.2, 60.4] total: 11800
## Weights on dummy points:
## range: [30.2, 121] total: 181000
## --------------------------------------------------------------------------------
## FITTED :
##
## Nonstationary Poisson process
##
## ---- Intensity: ----
##
## Log intensity: ~TMSB4X
## Model depends on external covariate 'TMSB4X'
## Covariates provided:
## TMSB4X: function(x, y)
##
## Fitted trend coefficients:
## (Intercept) TMSB4X
## -5.600218 -17.303762
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -5.600218 0.7246593 -7.020524 -4.179912 *** -7.728071
## TMSB4X -17.303762 13.5430663 -43.847684 9.240160 -1.277684
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept) TMSB4X
## -5.600218 -17.303762
##
## Fitted exp(theta):
## (Intercept) TMSB4X
## 3.697058e-03 3.055426e-08
##
## ----------- COX -----------
## Model: log-Gaussian Cox process
##
## Covariance model: exponential
## Fitted covariance parameters:
## var scale
## 2.036141 91.074313
## Fitted mean of log of random intensity: [pixel image]
##
## Final standard error and CI
## (allowing for correlation of Cox process):
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -5.600218 2.020185 -9.559707 -1.640729 ** -2.7721318
## TMSB4X -17.303762 35.721408 -87.316435 52.708912 -0.4844087D14 = fitcox(ty2_coor,'D14',ty2_time)## Warning in cf$x * x: longer object length is not a multiple of shorter object
## length## Warning in cf$`(Intercept)` - cf$x * x: longer object length is not a multiple
## of shorter object length## Warning in cf$y * y: longer object length is not a multiple of shorter object
## length## Warning in cf$x * x: longer object length is not a multiple of shorter object
## length## Warning in cf$`(Intercept)` - cf$x * x: longer object length is not a multiple
## of shorter object length## Warning in cf$y * y: longer object length is not a multiple of shorter object
## lengthsummary(D14)## Inhomogeneous Cox point process model
## Fitted to point pattern dataset 'unmark(ppppro)'
## Fitted by minimum contrast
## Summary statistic: inhomogeneous K-function
## Minimum contrast fit (object of class "minconfit")
## Model: Log-Gaussian Cox process
## Covariance model: exponential
## Fitted by matching theoretical K function to unmark(ppppro)
##
## Internal parameters fitted by minimum contrast ($par):
## sigma2 alpha
## 1.935681 164.157719
##
## Fitted covariance parameters:
## var scale
## 1.935681 164.157719
## Fitted mean of log of random intensity: [pixel image]
##
## Converged successfully after 79 function evaluations
##
## Starting values of parameters:
## sigma2 alpha
## 1.0000 15.2501
## Domain of integration: [ 0 , 105.9 ]
## Exponents: p= 2, q= 0.25
##
## ----------- TREND -----
## Point process model
## Fitted to data: X
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = X, trend = trend, rename.intercept = FALSE, covariates = covariates,
## covfunargs = covfunargs, use.gam = use.gam, forcefit = TRUE,
## improve.type = ppm.improve.type, improve.args = ppm.improve.args,
## nd = nd, eps = eps)
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Planar point pattern: 324 points
## Average intensity 0.00168 points per square unit
## Window: rectangle = [51.0357, 506.8714] x [106.7723, 530.5386] units
## (455.8 x 423.8 units)
## Window area = 193168 square units
##
## Dummy quadrature points:
## 40 x 40 grid of dummy points, plus 4 corner points
## dummy spacing: 11.39589 x 10.59416 units
##
## Original dummy parameters: =
## Planar point pattern: 1604 points
## Average intensity 0.0083 points per square unit
## Window: rectangle = [51.0357, 506.8714] x [106.7723, 530.5386] units
## (455.8 x 423.8 units)
## Window area = 193168 square units
## Quadrature weights:
## (counting weights based on 40 x 40 array of rectangular tiles)
## All weights:
## range: [30.2, 121] total: 193000
## Weights on data points:
## range: [30.2, 60.4] total: 12700
## Weights on dummy points:
## range: [30.2, 121] total: 180000
## --------------------------------------------------------------------------------
## FITTED :
##
## Nonstationary Poisson process
##
## ---- Intensity: ----
##
## Log intensity: ~TMSB4X
## Model depends on external covariate 'TMSB4X'
## Covariates provided:
## TMSB4X: function(x, y)
##
## Fitted trend coefficients:
## (Intercept) TMSB4X
## -2.730574 -95.114968
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -2.730574 0.3458162 -3.408362 -2.052787 *** -7.896027
## TMSB4X -95.114968 9.2593981 -113.263055 -76.966881 *** -10.272263
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept) TMSB4X
## -2.730574 -95.114968
##
## Fitted exp(theta):
## (Intercept) TMSB4X
## 6.518183e-02 4.921463e-42
##
## ----------- COX -----------
## Model: log-Gaussian Cox process
##
## Covariance model: exponential
## Fitted covariance parameters:
## var scale
## 1.935681 164.157719
## Fitted mean of log of random intensity: [pixel image]
##
## Final standard error and CI
## (allowing for correlation of Cox process):
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -2.730574 2.588156 -7.803267 2.342118 -1.055027
## TMSB4X -95.114968 61.670000 -215.985946 25.756010 -1.542322par(mfrow=c(4,2))
plot(D4, main="D4 type 2")## Warning in cf$x * x: longer object length is not a multiple of shorter object
## length## Warning in cf$`(Intercept)` - cf$x * x: longer object length is not a multiple
## of shorter object length## Warning in cf$y * y: longer object length is not a multiple of shorter object
## length## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function nameplot(D7, main="D7 type 2")## Warning in cf$x * x: longer object length is not a multiple of shorter object
## length## Warning in cf$`(Intercept)` - cf$x * x: longer object length is not a multiple
## of shorter object length## Warning in cf$y * y: longer object length is not a multiple of shorter object
## length## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function nameplot(D10, main="D10 type 2")## Warning in cf$x * x: longer object length is not a multiple of shorter object
## length## Warning in cf$`(Intercept)` - cf$x * x: longer object length is not a multiple
## of shorter object length## Warning in cf$y * y: longer object length is not a multiple of shorter object
## length## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function nameplot(D14, main="D14 type 2")## Warning in cf$x * x: longer object length is not a multiple of shorter object
## length## Warning in cf$`(Intercept)` - cf$x * x: longer object length is not a multiple
## of shorter object length## Warning in cf$y * y: longer object length is not a multiple of shorter object
## length## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name