hey everyone :), Subject: Re: Error message: object of type 'closure' is not subsettable I am writing a package which should calculate the binary logistic regression. The function itself work perfectly, but if I want to load the function from my package it gives me the above mentioned error. I have tried a lot of things, googled and so on, but I can not figure out what to change so it works. I also found out the mistake lies in the first function (logisticRegression)
logisticRegression <- function(x,y, threshold = 1e-10, maxIter = 100) { calcPi <- function(x,betaCoef) { betaCoef <- as.matrix(betaCoef) return(exp(x%*%betaCoef)/(1 + exp(x%*%betaCoef))) } #initial guess for beta (mostly we start with 0) betaCoef <- rep(0, ncol(x)) #some initial value which is bigger than the threshold for the loop initialValue <- 10000 #count of iteration to make sure its not an infinite loop iterCount <- 0 #iteration process (loop) while(initialValue > threshold) #convergence test { #calculates probabilities using the current estimate of beta p <- as.vector(calcPi(x, betaCoef)) #calculates W which is needed for the Hessian W <- diag(p*(1-p)) #calculates the gradient gradient <- t(x)%*%(y-p) #calculates the hessian hessian <- -(t(x)%*%W)%*%x #calculates beta betaCoef <- betaCoef - solve(hessian) %*% gradient #how much did we change beta? initialValue <- sum((solve(hessian) %*% gradient)^2) #to see if we reached the max # of iteration iterCount <- iterCount + 1 if(iterCount > maxIter) { stop("This is not converging") } } #df Totaldf <- length(y) - 1 Residualdf <- length(y) - length(betaCoef) degreesOfFreedom <- cbind(Totaldf, Residualdf) #fisher information, variance, standard error fisherInformation <- -hessian vCov <- solve(fisherInformation) stdError <- sqrt(diag(vCov)) #wald test statistic waldTestStatistic <- betaCoef/stdError pValue <- 2 * pnorm(-abs(waldTestStatistic), mean = 0, sd = 1) #loglikelihood and deviance logLikelihoodFunction <- y%*%x%*%betaCoef - sum(log(1+exp(x%*%betaCoef))) residualDeviance <- -2*logLikelihoodFunction #AIC akaikeIC <- residualDeviance + 2*ncol(betaCoef) #for R2 cV <- 2/nrow(x) #for C&S R2 #deviance residual sign1 <- replace(y, y == 0, -1) devianceResiduals <- sign1*sqrt(-2*(y*log(p) + (1-y)*log(1-p))) residualQuantiles <- summary(devianceResiduals) #fitted values fittedValues <- log(p/(1-p)) #what to return coefficients <- betaCoef output <- list() output$degreesF <- Residualdf output$nulldf <- Totaldf output$fittedValues <- fittedValues output$devianceResiduals <- devianceResiduals output$residualQuantiles <- residualQuantiles output$coefficients <- coefficients output$vCov <- vCov output$coefMat <- cbind(coefficients, stdError, waldTestStatistic, pValue) colnames(output$coefMat) <- c("Estimate", "Std.Err", "z value", "Pr(>|z|)") output$logLikelihoodFunction <- logLikelihoodFunction output$residualDeviance <- residualDeviance output$akaikeIC <- akaikeIC output$iterCount <- iterCount -1 output$p <- p output$cV <- cV output } logReg <- function(formula, data) { x <- model.matrix(formula, data) regressant <- model.frame(formula, data) y <- regressant[,1] constant <- matrix(rep(1, length(y))) nulllog <- logisticRegression(constant,y) nullDeviance <- nulllog$residualDeviance nullLogLikelihood <- nulllog$logLikelihoodFunction output <- list() output <- logisticRegression(x,y) output$formula <- formula output$nullDeviance <- nullDeviance output$nullLogLikelihood <- nullLogLikelihood output$call <- match.call() #define the S3 class class(output) <- "logReg" print(output) } logReg.default <- function(x, ...) UseMethod("logReg") print.logReg <- function(x, ...) { cat("\nCall:\n") print(x$call) cat("\nCoefficients:\n") colnames(x$coefficients) <- "" print(t(x$coefficients), digits = 5) cat("\nDegrees of Freedom:", paste(x$nulldf, "Total (i.e. Null);", x$degreesF, "Residual")) cat("\nNull Deviance:", paste(round(x$nullDeviance, digits = 2))) cat("\nResidual Deviance:", paste(round(x$residualDeviance, digits = 2))) invisible(x) } Thanks for your help! :D -- Sent from: http://r.789695.n4.nabble.com/R-help-f789696.html [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.