@@ -233,7 +233,7 @@ init_gam <- function(formula,
233
233
dl <- eval(inp , data , parent.frame())
234
234
if (! control $ keepData ) { rm(data )} # # save space
235
235
names(dl ) <- vars # # list of all variables needed
236
- var.summary <- mgcv ::: variable.summary (gp $ pf ,dl ,nrow(mf )) # # summarize the input data
236
+ var.summary <- variable_summary (gp $ pf ,dl ,nrow(mf )) # # summarize the input data
237
237
rm(dl )
238
238
239
239
G <- gam_setup(gp ,pterms = pterms ,
@@ -244,7 +244,7 @@ init_gam <- function(formula,
244
244
G $ model <- mf ;G $ terms <- terms ;G $ family <- family ;G $ call <- cl
245
245
G $ var.summary <- var.summary
246
246
247
- lambda <- mgcv ::: initial.spg (G $ X ,
247
+ lambda <- initial_spg (G $ X ,
248
248
G $ y ,
249
249
G $ w ,
250
250
family ,
@@ -276,7 +276,7 @@ init_gam <- function(formula,
276
276
return (G )
277
277
}
278
278
279
- # '@importFrom mgcv gam.side smoothCon get.var Rrank interpret.gam
279
+ # '@importFrom mgcv gam.side smoothCon get.var Rrank interpret.gam initial.sp
280
280
# '@importFrom stats .getXlevels model.matrix model.offset na.omit
281
281
# '@importFrom methods cbind2
282
282
# ' @noRd
@@ -824,3 +824,135 @@ clone_smooth_spec <- function(specb, spec){
824
824
}
825
825
specb
826
826
}
827
+
828
+ # ' @noRd
829
+ variable_summary <- function (pf , dl , n ){
830
+ v.n <- length(dl )
831
+ v.name <- v.name1 <- names(dl )
832
+ if (v.n ) {
833
+ k <- 0
834
+ for (i in 1 : v.n ) if (length(dl [[i ]]) > = n ) {
835
+ k <- k + 1
836
+ v.name [k ] <- v.name1 [i ]
837
+ }
838
+ if (k > 0 )
839
+ v.name <- v.name [1 : k ]
840
+ else v.name <- rep(" " , k )
841
+ }
842
+ p.name <- all.vars(pf [- 2 ])
843
+ vs <- list ()
844
+ v.n <- length(v.name )
845
+ if (v.n > 0 )
846
+ for (i in 1 : v.n ) {
847
+ if (v.name [i ] %in% p.name )
848
+ para <- TRUE
849
+ else para <- FALSE
850
+ if (para && is.matrix(dl [[v.name [i ]]]) && ncol(dl [[v.name [i ]]]) >
851
+ 1 ) {
852
+ x <- matrix (apply(dl [[v.name [i ]]], 2 , quantile ,
853
+ probs = 0.5 , type = 3 , na.rm = TRUE ), 1 , ncol(dl [[v.name [i ]]]))
854
+ }
855
+ else {
856
+ x <- dl [[v.name [i ]]]
857
+ if (is.character(x ))
858
+ x <- as.factor(x )
859
+ if (is.factor(x )) {
860
+ x <- x [! is.na(x )]
861
+ lx <- levels(x )
862
+ freq <- tabulate(x )
863
+ ii <- min((1 : length(lx ))[freq == max(freq )])
864
+ x <- factor (lx [ii ], levels = lx )
865
+ }
866
+ else {
867
+ x <- as.numeric(x )
868
+ x <- c(min(x , na.rm = TRUE ), as.numeric(quantile(x ,
869
+ probs = 0.5 , type = 3 , na.rm = TRUE )), max(x ,
870
+ na.rm = TRUE ))
871
+ }
872
+ }
873
+ vs [[v.name [i ]]] <- x
874
+ }
875
+ vs
876
+ }
877
+
878
+ # ' @noRd
879
+ initial_spg <- function (x , y , weights , family , S , rank , off , offset = NULL ,
880
+ L = NULL , lsp0 = NULL , type = 1 , start = NULL , mustart = NULL ,
881
+ etastart = NULL , E = NULL , ... ){
882
+ if (length(S ) == 0 )
883
+ return (rep(0 , 0 ))
884
+ nobs <- nrow(x )
885
+ if (is.null(mustart ))
886
+ mukeep <- NULL
887
+ else mukeep <- mustart
888
+ eval(family $ initialize )
889
+ if (inherits(family , " general.family" )) {
890
+ lbb <- family $ ll(y , x , start , weights , family , offset = offset ,
891
+ deriv = 1 )$ lbb
892
+ pcount <- rep(0 , ncol(lbb ))
893
+ for (i in 1 : length(S )) {
894
+ ind <- off [i ]: (off [i ] + ncol(S [[i ]]) - 1 )
895
+ dlb <- - diag(lbb [ind , ind , drop = FALSE ])
896
+ indp <- rowSums(abs(S [[i ]])) > max(S [[i ]]) * .Machine $ double.eps ^ 0.75 &
897
+ dlb != 0
898
+ ind <- ind [indp ]
899
+ pcount [ind ] <- pcount [ind ] + 1
900
+ }
901
+ lambda <- rep(0 , length(S ))
902
+ for (i in 1 : length(S )) {
903
+ ind <- off [i ]: (off [i ] + ncol(S [[i ]]) - 1 )
904
+ lami <- 1
905
+ dlb <- abs(diag(lbb [ind , ind , drop = FALSE ]))
906
+ dS <- diag(S [[i ]])
907
+ pc <- pcount [ind ]
908
+ ind <- rowSums(abs(S [[i ]])) > max(S [[i ]]) * .Machine $ double.eps ^ 0.75 &
909
+ dlb != 0
910
+ dlb <- dlb [ind ]/ pc [ind ]
911
+ dS <- dS [ind ]
912
+ rm <- max(length(dS )/ rank [i ], 1 )
913
+ while (sqrt(mean(dlb / (dlb + lami * dS * rm )) * mean(dlb )/ mean(dlb +
914
+ lami * dS * rm )) > 0.4 ) lami <- lami * 5
915
+ while (sqrt(mean(dlb / (dlb + lami * dS * rm )) * mean(dlb )/ mean(dlb +
916
+ lami * dS * rm )) < 0.4 ) lami <- lami / 5
917
+ lambda [i ] <- lami
918
+ }
919
+ }
920
+ else {
921
+ if (is.null(mukeep )) {
922
+ if (! is.null(start ))
923
+ etastart <- drop(x %*% start )
924
+ if (! is.null(etastart ))
925
+ mustart <- family $ linkinv(etastart )
926
+ }
927
+ else mustart <- mukeep
928
+ if (inherits(family , " extended.family" )) {
929
+ theta <- family $ getTheta()
930
+ Ddo <- family $ Dd(y , mustart , theta , weights )
931
+ mu.eta2 <- family $ mu.eta(family $ linkfun(mustart ))^ 2
932
+ w <- 0.5 * as.numeric(Ddo $ Dmu2 * mu.eta2 )
933
+ if (any(w < 0 ))
934
+ w <- 0.5 * as.numeric(Ddo $ EDmu2 * mu.eta2 )
935
+ }
936
+ else w <- as.numeric(weights * family $ mu.eta(family $ linkfun(mustart ))^ 2 / family $ variance(mustart ))
937
+ w <- sqrt(w )
938
+ if (type == 1 ) {
939
+ lambda <- mgcv :: initial.sp(w * x , S , off )
940
+ }
941
+ else {
942
+ csX <- colSums((w * x )^ 2 )
943
+ lambda <- rep(0 , length(S ))
944
+ for (i in 1 : length(S )) {
945
+ ind <- off [i ]: (off [i ] + ncol(S [[i ]]) - 1 )
946
+ lambda [i ] <- sum(csX [ind ])/ sqrt(sum(S [[i ]]^ 2 ))
947
+ }
948
+ }
949
+ }
950
+ if (! is.null(L )) {
951
+ lsp <- log(lambda )
952
+ if (is.null(lsp0 ))
953
+ lsp0 <- rep(0 , nrow(L ))
954
+ lsp <- as.numeric(coef(lm(lsp ~ L - 1 + offset(lsp0 ))))
955
+ lambda <- exp(lsp )
956
+ }
957
+ lambda
958
+ }
0 commit comments