title: Kalman- Bucy filter: prior mean change

text: I have a question on Kalman-Bucy filter: the prior distribution is $g \sim N(0,σ_g^2 )$, signal is $ds=(μ+g_t )dt+σdZ_t$, posterior distribution becomes $g_t \sim N((\hat{g_t},\hatσ_t^2)$. Here,$σ_g,μ,σ$ are constants. Now, if I change the prior mean from 0 to a, where a is a constant, other things remain the same, may I know how the posterior distribution would change? Thanks a lot!!


title: spatial autocorrelation for time series data

text: I have a 20-yr dataset of an annual count of species abundance for a set of polygons (~200 irregularly shaped, continuous polygons). I have been using regression analysis to infer trends (change in count per year) for each polygon, as well as aggregations of polygon data based on management boundaries. I am sure that there is spatial autocorrelation in the data, which is sure to impact the regression analysis for the aggregated data. My question is - how do I run a SAC test for time series data? Do I need to look at the SAC of residuals from my regression for each year (global Moran's I)? Or can I run one test with all years? Once I've tested that yes there is SAC, is there an easy was to address this? My stats background is minimal and everything I've read on spatio-temporal modeling sounds very complex. I know that R has a distance-weighted autocovariate function - is this at all simple to use? I'm really quite confused on how to assess/addess SAC for this problem and would very much appreciate any suggestions, links, or references. Thanks in advance!


title: Correlation of elevation data and elevation uncertainty as a function of range in a DEM

text: I have a gridded Digital Elevation Model (DEM) which provides terrain elevation ($z$) at cartesian coordinates $(x,y)$. Together with the DEM I know the RMS of the elevations. As far as I understand the RMS can be applied to any point in the DEM and the standard deviation of the elevations is stationary. Now, if I take two points from the DEM, the elevations at these points are not independent, but are auto-correlated if the distance between the points is below a certain range. I can build a variogram from the DEM elevation data which provides the function of elevation variation as a function of distance. For distances beyond the range the elevation variation is independent and just the mean of all elevations in the terrain. However, what I really want to know is the uncertainty of the elevation at different close distances. I would assume that the uncertainty is equal to the RMS at large range, but apporaches zero if the distance apporaches zero. What I need is the function vor values between zero-distance and the range. I would just apply the variogram function of the elevation to the uncertainty and scale it from the mean elevation to the RMS, but I guess this is a pretty much simplified approach.


title: Model evaluation and comparison for selecting the best model

text: When comparing results obtained with different models in R, what should I look for to select the best one? If I use for example the following 4 models applied to the same presence/absence sample taken from a species dataset, with the same variables:   * Generalized linear model   * Generalized additive models Classification   * Regression Tree   * Artificial Neural Networks Should I compare all methods by AIC, Kappa, or cross-validation? Will I ever be certain of selecting the best model? What happens if I compare those 4 models prediction with a Bayes factor? Can I compare them?


title: Relation between causal inference and prediction (classification and regression)

text: I was wondering what relation and differences are between causal inference and prediction (classification and regression)? For example, In prediction, we have predictor/input variables and response/output variables. Does that mean that there is causal relation between input and output variables? So does prediction belong to causal inference? If I understand correctly, causal inference considers to estimate conditional distribution of one random variable given another random variable, and often use graphical models to represent the conditional independence between random variables. So causal inference in this sense isn't prediction, is it? Thanks and regards!


title: How to devise a relationship matrix without diagonal constant when fitting mixed model in R?

text: I have the following small data:               father <- c(NA, NA, NA, 1, 1, 1, 2, 2, 2, 3, 4)     mother <- c(NA, NA, NA, 2, 2, 1, 1, 1, 3, 4, 4)     idd <- c(1:11)     yvar <- rnorm(11, 5, 2)     xvar <- c(1,1,2, 1,2,1, 1,1, 2, 1,2)     mydf <- data.frame (father, mother, idd, yvar, xvar)          require(kinship)     kmat <- kinship(mydf$idd, mydf$father,mydf$mother)                     1     2    3     4     5     6     7     8      9     10    11         1  0.500 0.000 0.00 0.250 0.250 0.500 0.250 0.250 0.0000 0.1250 0.250         2  0.000 0.500 0.00 0.250 0.250 0.000 0.250 0.250 0.2500 0.1250 0.250         3  0.000 0.000 0.50 0.000 0.000 0.000 0.000 0.000 0.2500 0.2500 0.000         4  0.250 0.250 0.00 0.500 0.250 0.250 0.250 0.250 0.1250 0.2500 0.500         5  0.250 0.250 0.00 0.250 0.500 0.250 0.250 0.250 0.1250 0.1250 0.250         6  0.500 0.000 0.00 0.250 0.250 0.750 0.250 0.250 0.0000 0.1250 0.250         7  0.250 0.250 0.00 0.250 0.250 0.250 0.500 0.250 0.1250 0.1250 0.250         8  0.250 0.250 0.00 0.250 0.250 0.250 0.250 0.500 0.1250 0.1250 0.250         9  0.000 0.250 0.25 0.125 0.125 0.000 0.125 0.125 0.5000 0.1875 0.125         10 0.125 0.125 0.25 0.250 0.125 0.125 0.125 0.125 0.1875 0.5000 0.250         11 0.250 0.250 0.00 0.500 0.250 0.250 0.250 0.250 0.1250 0.2500 0.750      This is example of inbred pedigree in plants where anything can be father or mother or vice versa. Even a single individual can be hermaphrodite (both mother and father) under-selfing. Even in this situation the relatedness between individual 1 with individual 1 should be 0.5. Is there any correct-way to correct this problem of kinship matrix estimation? The mixed model fitting using `lmekin` function from kinship package is giving the following error:               require(kinship)      model1 <- lmekin(yvar ~  xvar , random = ~ 1|idd, varlist=list(kmat), data = mydf)         Warning message:     In coxme.varcheck(ncluster, varlist, n, gvars, groups, sparse, rescale,  :       Diagonal of variance matrix is not constant


title: SVM data normalization... what about classifying new (training) data?

text: I've got a big doubt about SVM classification task (and more in general classification task), about **data normalization**. Let's suppose I've a **SVM trained with normalized** data, and **new data to classify**. **1** ) **How do I normalize new data**? Please note that I don't know them when I normalized and trained my SVM. **2** ) **Which is the best/proper normalization method**? Min-max or zero- mean+variance? A possible solution that I thought is: once new data arrives, and as we are working with SVs (that are, part of the training data), we can de-normalize the SVs, re-calculate min-Max/mean-var of the new WHOLE dataset, and normalize the new data and re-normalize the SVs. What about this? Thanks in advance, Ivano


title: Measure of accuracy or measure of agreement for continuous data with true values

text: I have two continuous variables. Variable $X$ represents the true values (basically it is atmospheric pressure measured via a barometer). Variable $Y$ represents a way to approximate the atmospheric pressure based on altitude, temperature, and some other physical parameters. I have $n=17$ values. The question now is, how close the approximations are to the true values, i.e., how close are $X$ and $Y$. It seems to me, that the Bland-Altman-Plot (as descriptive tool) and the Concordence Correlation Coefficient are appropriate to assess the agreement of the two variables and thus the accuracy of $Y$. My only concern is, that one of the variables represents the true value (i.e. $X$) and my understanding is, that the aforementioned methods apply to situations when I'm comparing two measurements without actually knowing the true values, i.e., inter-rater reliability. My question is whether I can use those two methods and if there are other (more) appropriate methods.


title: How to extract bootstrap forecasts from forecast.Arima in the R package forecast?

text: In a forecasting context, I have two different time series: y = {y1, y2, .., yn} and z = {z1, z2, ..., zn}. In R, let's say that the two series are as follows: y <\- arima.sim(list(order = c(1,0,0), ar = 0.7), n = 100) z <\- arima.sim(list(order = c(1,0,0), ar = 0.5), n = 100) Using the auto.arima() function in the R package "forecast", I find the ARIMA models which can be used to describe the temporal dynamics of the two time series: fit.y <\- auto.arima(y) fit.z <\- auto.arima(z) Next, I use the forecast.Arima() function to produce one-step-ahead forecasts and associated forecast intervals for each series using the bootstrap method: for.y <\- forecast(fit.y, h=1, bootstrap=TRUE, npaths=5000) for.z <\- forecast(fit.z, h=1, bootstrap=TRUE, npaths=5000) The forecast function returns a point forecast and a forecast interval. However, the 5000 "bootstrap forecasts" used to generate the forecast intervals are not reported by the forecast function. Moreover, it seems (?) that the forecast function uses the normal distribution quantiles to construct the reported forecast interval, even though the use of bootstrapping would seem to call for the use of the empirical quantiles derived from the distribution of "bootstrap forecasts". My interest in the 5000 "bootstrap forecasts" is motivated by the need to construct a forecast interval for y + z. If I had access to the "bootstrap forecasts" for y and those for z, I could simply add them up to obtain the bootstrap distribution of the one-step ahead forecast of y + z. Is there a way I could extract the 5000 "bootstrap forecasts" from the forecast() function? The other thing that seems to pose challenges in using forecast() on an ARIMA model which includes a drift. It seems to me that the forecast horizon used by forecast() is somehow connected to the xreg parameter. If the ARIMA model includes a drift, the xreg parameter will have the same dimension as the original time series(), so the forecast function will automatically predict up to n steps ahead. If I wanted to predict just one step ahead, how would I bypass this problem? Thank you!


title: POISSON regression with offset variable

text: Our problem is to determine if there is a relationship between the return on equity of firms (ROE) and the presence (number of indicators used) (Y) of a specific type of indicator showing up in the annual report. Initially we had 118 observations and after some adjustments we had 104 firms. In order to verify if there is a linear relationship between Y (Poisson distributed) and the predictor ROE we built 8 clusters based on ROE’s distribution and we computed for each class the mean of ROE (ROEm), the mean of the number of indicators used (ym) and, finally, the total amount of indicators used by each class. Then we plotted in a graph the log_ym and ROEm, their relationship was not linear but quadratic. ![scatterplot](http://i.stack.imgur.com/YzIhq.png) Thus we decided to apply a Poisson regression using as the dependent variable Y and as predictor ROE and ROE^2. This model was not so good: pseudo R^2 was very small (about 0,02) and the hypothesis of over-dispersion was accepted. ![enter image description here](http://i.stack.imgur.com/xq8Rw.png) As we read in some papers, we concluded that there is a tendency for observations to cluster. At this point we decided to apply the Poisson regression on clusterized data using as dependent variable the total amount of indicators used by each class (ym is not discrete) and as predictor ROEm and ROEm^2. As the offset variable we used the number of cases per each cluster. The output of the model is summarized below: ![enter image description here](http://i.stack.imgur.com/7G761.png)   * Is our approach correct?   * Are there simpler methods?   * Are 8 clusters few?    * Results get worst if we use for example 14 clusters instead of 8; does that mean that the model doesn’t fit well?


title: How can I improve the predictive power of this logistic regression model?

text: I am using SPSS to analyze a data set which aims to predict whether individuals have cancer based on five symptoms (a, b, c, d, e). In this data set most individuals have cancer. I ran a Binary Logistic Regression and got the following output: ![Block 0](http://i.stack.imgur.com/aXuqu.png) This tests the model with which only includes the constant, and overall it predicted 91.8% correct. I understand that the fact that I have significant predictors in the "Variables not in the Equation" table means that the addition of one or more of these variables to the model should improve its predictive power. I then looked at the model after all the predictors were included: ![http://i.imgur.com/Y5WdePl.png](http://i.stack.imgur.com/77Fal.png) The prediction is only very slightly different. Now it's predicting that two individuals will not have cancer. The overall percentage correct remains 91.8%.   * Why did no improvement occur, despite the predictors being significant?   * Where should I go from here with this data set? Is it possible to improve the predictive power of the model without including new predictors?   * How should I assess the model? Is the fact that it doesn't improve over the model with only the constant evidence that the model is useless? Full output viewable here. Data set downloadable here as a google doc.


title: Assessing independence of observations using intraclass correlation when some groups have small group sample sizes

text: ### Context I was talking to a researcher in the following situation.   * Participants (n = 500) were sampled from schools.   * Participants came from around 50 different schools.   * The number of participants per school varied with some schools supplying 20 or 30 participants, but a few schools only supplying 3 or 4 or 5 participants.   * The researcher was trying to assess whether there was a substantial violation of the independence of observations assumption. Thus, they were looking at the intra-class correlation of core outcome variables based on the effect of school. ### Question   * Is it problematic to include groups in an intra-class correlation analysis with small numbers per group? If so, what are the implications of this? What might a rule of thumb regarding a minimum number for inclusion be?   * To take it to the extreme, what is some groups only supplied a single participant?   * If there are groups with samples sizes below a given threshold, what is a good subsequent course of action? remove group from analysis? collapse small groups into an "other" group?   * How would any recommendations given relate to the assessment of the independence of observations assumption?


title: Multi-parameter log likelihood of Normal distribution for two separate samples

text: I have been given two sets of independent random variables distributed by two different normal distributions $X_1,...,X_n \sim N(\theta_1,1)$ and $Y_1,...,Y_m \sim N(\theta_2,1)$. And have been asked to find the log likelihood function for the parameter $\theta = (\theta_1,\theta_2)$. I was wondering if a correct approach would be to have $l(\theta) = \sum_{i=1}^n(log(f(\theta_1;x_i)) + \sum_{i=1}^m(log(f(\theta_2;y_i)))$, where $f(\theta_1;x_i)$ is the pdf for $X_i \sim N(\theta_1,1)$ or am I going about this completely the wrong way?


title: Random number generation using t-distribution or laplace distribution

text: I am a newbie in stat. I am completing my thesis in Evolutionary algorithm. I have to generate some random numbers from T-distribution or Laplace distribution. How can I do this? An easy and simple explanation would be appreciated. Thanks in advance.


title: Measuring non-linear dependence

text: Covariance between two random variables defines a measure of how closely are they linearly related to each other. But what if the joint distribution is circlular? Surely there is structure in the distribution. How is this structure extracted?


title: Optim result highly dependent on starting value

text: I want to fit a standardized Student's-t distribution. The log-likelihood is given by: \begin{align*} log \mathcal{L}(\nu | l_1,...,l_n)=\sum_{i=1}^n \left( log \left( (\pi (\nu-2))^{-\frac{1}{2}}\Gamma \left(\frac{\nu}{2} \right)^{-1} \Gamma \left(\frac{\nu+1}{2} \right) \left(1+\frac{l^2}{\nu-2} \right)^{\text{$-\frac{1+\nu}{2}$}}\right)\right) \end{align*} I try do to this with the optim command. My R code is (data):               # log likelihood     pinumber<-3.141592653589793     startvalue<-55          loglikstandardizedt <-function(par){     if(par>0) return(-sum(log((pinumber*(par-2))^(-1/2)*gamma(par/2)^(-1)*gamma((par+1)/2)*(1+standresidalvewma^2/(par-2))^(-(1+par)/2))))     else return(Inf)     }          optim(startvalue,loglikstandardizedt, method="BFGS")     param = optim(startvalue,loglikstandardizedt, method="BFGS")$par      I can control it via looking at the plot:                 denstiystandtresid<-function (x) (pinumber*(param-2))^(-1/2)*gamma(param/2)^(-1)*gamma((param+1)/2)*(1+x^2/(param-2))^(-(1+param)/2)              plot(density(standresidalvewma),ylim=c(0,0.8))         curve(denstiystandtresid,col="red",add=TRUE)      I have the problem, that the output of optim, that means the value for "par" is highly dependend on the starting value. What is the right starting value? And what is the best value for my $\nu$? The variation is very large, I get values from 27 up to 215, so this is a huge difference?


title: Do zero counts need to be adjusted for a likelihood ratio test of poisson/loglinear models?

text: If there are 0's in the contingency table and we're fitting nested poisson/loglinear models (using R's `glm` function) for a likelihood ratio test, do we need to adjust the data prior to fitting the glm models (e.g. add 1/2 to all the counts)? Obviously some parameters cannot be estimated without some adjustment, but how does the adjustment/lack of adjustment effect the LR test?


title: How to eliminate high multicollinearity with a continuous moderating variable, and a categorical independent variable

text: I am looking at whether Corporate Venture Capital-backed firms (1) perform better than Independent Venture capital - backed firms(0) in their POST-IPO performance. My assumption is that this relationship is either strengthened or weakened by the amount of Experience (continuous). After creating a new variable for Corporate VC*Experience, I then put these values in a regression and it gives me a VIF of over 23.


title: Estimating the mean and variance of a stationary time series

text: For a stationary time series can I use the sample mean as an estimate for the time series mean? Also what is the recommended way to estimate the variance? The series is auto- correlated.


title: How to compute the difference between 2 similarity matrices?

text: I had asked this question on Mathoverflow.net and someone suggested that this might be a better forum. So here is the question again Hello, I have two n*n correlation matrices with values ranging between -1,1. (2 correlation matrix because I have the same n terms under 2 different conditions) Taking the absolute of the correlation matrix, I then transformed the correlation into a similarity matrix where values range between [0,1]. The similarity was computed using the formula here: ![enter image description here](http://i.stack.imgur.com/QYxMl.gif) where lij = summation(aiu.auj) and aiu= is correlation value between ith and uth element in the correlation matrix. ki,kj in our case is n-1. This gives me 2 similarity matrices. My question is how can I compute the difference between the similarity matrices (SIM)? Will doing a simple abs(SIM1-SIM2) make statistical/mathematical sense? Is there a robust way to find/quantify difference between two similarity matrices? Any ideas will be truly appreciated. Thanks in advance.


title: How do I plot loglikelihood functions of the Cauchy distribution in R?

text: I need to determine whether or not there is always a unique root of the log- likelihood equation l'(x,0)=0 for different samples sizes of the Cauchy distribution. But I am unsure about how to actually define the log-likelihood function R and about actually being able to plot these graphs. Any advice would be appreciated!


title: What's a good visualization for Poisson regressions?

text: I want to associate code defects with code complexity metrics like closeness. One common model is to view this as a Poisson process, where the duration is how much time is spent coding and the density is a function of the code complexity. I am able to do a regression and get significance values etc. However, it's hard for me to visualize the results (and even harder for my less mathematically inclined colleagues). Is there a good way to view this to look at things like outliers, if it's a linear trend etc.? (Links to R packages are appreciated.) For example, I could just plot `Defects / Time ~ Complexity` but this is fairly noisy and since `Defects` is discretized and very small, it's hard to view the trend. One thing I was thinking of is that I could slice the data into quantiles and then do a regression per quantile and plot the resulting densities - I'm not certain how valid this is though, especially since my data are not normally distributed so it might mislead people about what a quantile is.


title: Practical time series advice

text: Looking for some general time series advice! I have collected performance data at fixed time intervals from a 'shared system' with the aim of investigating the affect of the sharing on the performance of my 'slice' of the system. The performance as claimed should be consistent over time, the actions of the neighbors are unobserved apart from the affect on my performance. Unfortunately I cant post time plots collected (over a two week period) or the plot of differences as I'm new here. There is no reason to believe that the series should have any deterministic trends and my initial belief is that a good model would be               P(t)=constant + random variation      However, the performance series certainly appears non-stationary (and rejects a KPSS test of stationary with no trend) and the correlogram decays very slowly whilst the pacf shows no non-significant correlations after 15 lags. Looking at my time plot and the differenced series the variance appears to be changing. I'm new to time series analysis and would appreciate any pointers and what might be a useful direction to take the analysis in. My understanding (which may well be wrong) is that non constant variance together with serial correlation in a series gives rise to local trends - stochastic trends? Would a ARCH/GARCH model be applicable here or should I try ARIMA models first? I'm less interested in forecasting future performance and more interested in understanding the structure of the data. Any help appreciated!


title: Minimum number of observations per variable for linear regression or MARS model

text: How many cases are required for each variable in order to build a linear regression model and/or a Multivariate Regression Splines Model? Also,   1. Is it a rule of thumb, or does there exist a statistical justification?   2. Any bibliographical references?


title: Logging similarities between vectors with R

text: I'm trying to write a program that automatically groups similarities between vectors. The vectors are comprised of point coordinates. For example (assuming X, Y, and Z are numbers):   Data Set 1: [1, 8, **X1** , **X2** , 6, **X3** , **X4** , 8, 5, 4, 9, 4, **Y1** , **Y2** , **Y3** , **Y4** , **Y5** , 1, 5, 7]   Data Set 2: [ **X1** , **X2** , **X3** , **X4** , 6, 3, 5, 4, 3, 6, **Y1** , **Y2** , **Y3** , 7, **Y4** , **Y5** , 7, 3, 5, 7]   Data Set 3: [ **X3** , 5, 6, 4, 3, 6, **X1** , **Y3** , 3, 5, 6, **X4** , 2, 6, 7, **X2** , **Y2** , 5, 6, **X4** , **Y4** ] X1, X2, X3, and X4, and Y1, Y2, Y3, Y4, and Y5 don't match perfectly with each other between Data Set 1 and 2, but they're within very close proximity, and in the same order, unlike the example in Data Set 3. Is there an algorithm that can do this kind of grouping?


title: Finite sample concentration results for partial correlation

text: Let $X = (X_1, \ldots, X_p) \sim F$ be a multivariate random vector. We have n iid observations from F, $\\{x_i\\}_{i=1}^n$. Let $\hat r_{ab}$ be the empirical estimate of the partial correlation between $X_a$ and $X_b$ given all the other components of F. What results are known for $P[|\hat r_{ab} - r_{ab}| > t] \leq f(t)$ where $r_{ab}$ is the population partial correlation and $f(t)$ is some decreasing function of t? What if $F$ is multivariate Gaussian? References would be welcomed.


title: How do I choose a model/estimator?

text: I've learned many different types of ways to find models with regression, and many different model to find clusters with clustering algorithms however ever since I started with ML, I've always wondered how would one choose to use one algorithm over another? Should we always measure the errors of each model then make our selection based on that? I've learned things like SVM, Perceptron algorithm, Feed forward Neural Networks, Boltzmann Machines however when should I choose one over the other? This has always really confused me. So from what I understand, a data scientist would just,               for i in all_algorithms:         use_algorithm(data)     use_best_algorithm(data)      Side Idea: I was also reading that the algorithm doesn't really matter but it is the amount of data that you have which really matters. Thank you in advance!


title: Mean and SD for distance between occurrences of a given value in sequence drawn at random

text: Given an infinite sequence of values, selected at random (with replacement) from a set of N values, I choose some value X from the set of values. Then, starting from the first occurrence of X in the sequence, I count the number of values between each subsequent occurrence of X, giving a new sequence. What is the mean and SD of this new sequence? For example, if N=2 and the values are 0 and 1, my infinite sequence might start:               0010101101111000101010011001001...      if I choose X to be 1 then the new sequence is:               0010101101111000101010011001001...         1 10 1000   3 1 1  20  2  2...       or the number of 0s that separates each successive 1. Thanks. [I expected this to have a mean and variance of N (vaguely remembering something about Poisson stats from school), but some experimenting seems to suggest the SD is N.]


title: Is a partial F-Test on a model reduced by only one variable valid?

text: For a recent project, I used multiple linear regression to model data. I attempted to choose between my initial full model and a reduced model by performing a partial F-test. The models used were the following: Full model: $\hat{Y}$ = $B_0$ + $B_1$$x_1$ + $B_2$$x_2$ + $B_3$$x_3$ Reduced Model: $\hat{Y}$ = $B_0$ + $B_2$$x_2$ + $B_3$$x_3$ The reduced model has only one covariate removed - $x_1$. The partial F-Test returned a F-value of 2.162 and a P-value greater than 0.1, so I did not reject the null hypothesis and concluded that the addition of $x_1$ in the model may not make any significant linear contribution to the prediction of $Y$. I was later told by my teacher that I should not run a partial F-Test on a reduced model of only one less variable. The teacher said that not only should I not do such a thing, but also that doing so would be wrong. I was told to instead just compare the two models using each model's $R^2$ values, F-values, and T-values. I can see how the partial F-test in this case might not provide helpful information, since the null hypothesis it tests is that $B_1 = 0$, which is the same thing that the T-test on the parameter $x_1$ tests (in the full model). [Note: I'm not sure what the exact relationship between this T-value and the partial F-test's F-value is.] Nonetheless, why would it would be considered wrong to use and interpret the results from the partial F-test?


title: How should one control for group and individual differences in pre-treatment scores in a randomised controlled trial?

text: Andrew Gelman, in the book he wrote with Jennifer Hill, states in Chapter 9, (section 9.3), on page 177: > It is only appropriate to control for pre-treatment predictors, or, more > generally, predictors that would not be affected by the treatment (such as > race or age). This point willl be illustrated more concretely in Section > 9.7... And there (9.7 is entitled "do not control for post treatment variables") he discusses the problem of measuring mediating variables, rather than the pre- post change problem directly. It is important to state here that I think Gelman/Hill is a brilliant text... And I'm thoroughly enjoying understanding it. However, this bit piqued my interest, as it brings to mind Everitt & Pickles's approach to the same problem. Everitt is of the opinion that using a change score (Score B - Score A) will tend to bias your findings in favour of the treatment, whereas including baseline scores in the model is more conservative. They back this up with a simulation - it's pretty persuasive. My understanding up to here has been that what you are controlling for is group differences in baseline scores that might cause the apparent treatment effect to be greater than it is, or to exist, when it does not. It is also my understanding that this is because regression to the mean is at work, so that higher baseline scores will be associated with greater decreases and vice versa, independent of the treatment effect. Everitt is strenuously against "change scores", and Gelman seems to be advising against Including baseline scores in the model. However, Gelman demonstrates this over the next 2-3 pages, including pre-test scores as a predictor. He gives the caveat that you then get a range of plausible treatment effects that are conditional on the pre-test score, not a range of treatment effects representing merely uncertainty in the effects. My opinion is that using "change scores" seems not to really be doing anything _about_ regression to the mean, whereas including the baseline score as a predictor allows baseline group differences to _cancel out_ , essentially introducing a covariance structure. I'm a doctor and I have to make real decisions about which treatments work. So what should I do? Include each person's baseline scores or use "change scores"?


title: Comparing means and variance of one group in two conditions

text: Probably a trivial questions, but I've been wrapping my head around it for a while... I have a small sample (N=7) that performed a test in two conditions. No assumption can be made about the eventual normality of the population. I would like to test if the change in conditions significantly affected the results of the tests across subjects. In this case a t-test on the means for the two conditions across subjects should not be applied. Am I right ? Would a Mann U-test or Wilcoxon signed-rank test be more suitable ? EDIT : (I am aware that these two test would compare the medians of the results) Thanks


title: Removing data above the mean

text: I am working with a data set, for which to remove noisy data I take the average of the sample itself and cut anything above that average. What I'm trying to understand here is does this have a name like I don't think this is mean filtering - but it would be helpful to know what type of filtering process this is since I need to explain it in a paper as for why I did it. Or is this just outlier removal?


title: Strangely imbalanced dataset

text: I'm new to Machine Learning and this forum. I have a beginner's doubt regarding imbalanced dataset. Here it goes: I have a binary classification task, where I'm more interested in accurately classifying the positive class (which is in minority in the target population). Unlike the common problem of **_not_** having enough positive (minority) class instances in the training set, my training dataset contains the positive class in **_majority_**. Here's the target population composition (which I'd expect to find in the environment where my classifier/model would be deployed):   * Positive Class: ~35%   * Negative Class: ~65% Here's my Training Set Composition:   * Positive Class: ~95%   * Negative Class: ~5% As my training set composition drastically differs from the target population composition, will the classification algorithm fail to generalize when classifying instances from the target population? As I mentioned earlier, I'm more interested in accurately classifying the positive class instances, which my training set has in abundance. I read the following description in a publication on imbalanced datasets: "The purpose of machine learning is for the classifier to estimate the probability distribution of the target population. Since that distribution is unknown we try to estimate the population distribution using a sample distribution. Statistics tells us that as long as the sample is drawn randomly, the sample distribution can be used to estimate the population distribution from where it was drawn. Hence, by learning the sample distribution we can learn to approximate the target distribution." Since my training dataset can **_not_** be considered as a random sample of the target distribution, will this affect the generalization power of my classifier? If so, what shall be done to avoid this? Over/under-sampling? Cost Matrices? PS: I tried to search previous posts about issues similar to mine, but all of them dealt with the problem of not having sufficient examples for the minority class (which is the exact opposite of my scenario). Thanks in advance -S


title: Error correction model (to test for asymmetry) with stationary I(0) variables

text: I have price series which are all stationary without taking any difference --> I(0). Can I still perform an ECM model to test for asymmetry? For example: Y= constant X; taking the residuals and separate this term in negative and positive ones (ECT+ and ECT-). Then: D.Y=cons LD.X LD.Y ECT- ECT+ Is this correct for I(0) variables?


title: Dynamic Time Warping - Why do we need the lowest cost path

text: My question is connected with the famous data mining algorithm Dynamic Time Warping also known as only DTW. According to wikipedia, the distance between two sequences is the (i-th,j-th) element of the DTW matrix. But I see in the reseach papers that for classification we need two find the lowest cost path/warping path, which is the sum of all cells in this path (not sure). Can't we use only the value of the DTW[i, j]? Also, I see that many people are dividing the distance by a constant and say that this is a normalization. Why do they do this?


title: Correctness of regression with ARIMA errors model and coefficient interpretation issues

text: I am trying to forecast electricity consumption in GWh for 2 years ahead (from June 2013 ahead), using R (the forecast package). For that purpose, I tried regression with ARIMA errors. I fitted the model using the `auto.arima` function, and I used the following variables in the `xreg` argument in the `forecast.Arima` function: \- Heating and Cooling Degree Days,   \- Dummies for all 12 months and   \- Moving holidays dummies (Easter and Ramadan) I have several questions regarding the model: 1) Is it correct to use all 12 dummies for monthly seasonality, since when I tried to include 11, the function returned error. The `Auto.arima` function returned the model ARIMA(0,1,2) 2)The model returned the following coefficients (I won't specify all of them as there are too many coefficients):               ma1      ma2     HDD     CDD   January  February  March     April     -0.52 -0.16      0.27    0.12  525.84   475.13    472.57    399.01      I am trying to determine the influence of the temperature component over electricity load. In percentages, (interpreting the coefficients just as with the usual regression) the temperature components (`HDD`+`CDD`) account for 11,3% of the electricity consumption. Isn't this too little, considering the fact that the electricity consumption is mostly influenced by the weather component? On the other hand, taking look at the dummies' coefficients, it turns out that the seasonality accounts for the greater part of the load. Why is this? Is the model completely incorrect? I tried linear regression, and the temperature component accounts for 20%, but it is still a low percentage. Why is this? 3) I am obviously making some mistakes in the use of `forecast.Arima` or the plot function parameters since when I plot the forecasts, I get a picture of the original time series which is continued (merged) with the forecasts for the whole time series period (from 2004 until 2015). I don't know how to explain this better, I tried to paste the picture, but it seems I cannot paste pictures here.


title: Two Way ANOVA repeated measurement

text: I am trying to compare the three replicate samples taken at different time point to see if the replicate one differ from replicate two and three and visa verse. My interest is to see if the three replicates represent the same community by have the same species abundance or not? From the reading I have done, I think that Two Way ANOVA repeated measurement would do the job for me but I am not quite sure that would be right? Please see the link of my data by following the link: https://gist.github.com/plxsas/9ec2a528f94bd2dff46d Can anybody advice me on that and how to achieve it in R, please?


title: Measuring Similarity of two Transient Signals

text: I am trying to measure the similarity between two signals, i.e. how identical they are. A complete overlap (two lines forming a single line) should give me a value of 1, the greater the distance, between the lines and the less similar the overall trend, the closer the value should go to 0. I have tried measuring standard deviation $\sigma$ and different variants of correlation $r$ using MATLAB functions `std()`, `xcorr()` and `corr2()`. They did not do the trick, as ,e.g. $\sigma_{\phi_{A}}$ was $>$ than $\sigma_{\phi_{B}}$ (see figure below). I was thinking about using some sort of _weighted average_ between $\sigma$ and $r$. Another possibility would be to measure the difference between areas under the curves (but then different intersecting graphs could have the same area). There are related questions to this, but they are either: (i) discussing the generic concept of _similarity_ or (ii) suggest using `xcorr()` (different variants of which did not work for me). ![Link to an image of two signals](http://oi43.tinypic.com/2s8h95k.jpg)


title: How to analyse a risk assessment questionnaire where each risk is rated for both probability and impact?

text: I would like your opinion on how to analyze a questionnaire on risk assessment in construction projects. For each question, concerning a specific risk for the project, there are two answers: An estimate for the probability of occurrence for the risk (on the scale from 1 to 5) and one for the estimated impact of the risk to the project (also in the range 1-5). A common way to combine these two values is by multiplying them. Then, according to a chosen (by the organization) "Probability/Impact" matrix, I could determine which combinations of probability and impact result in a risk’s being classified as "high risk", "moderate risk", or "low risk". I could then use these qualitative responses to analyze my data. However, I do not want to lose the information about a risk's probability because a risk classified as "high" should be treated differently, if this arises from high probability instead of high impact. The impact cannot be changed but the probability might be reduced by taking appropriate measures. What would the form of the data matrix be, if I want to keep a qualitative and a quantitative response for each question?. What methods do you suggest for the analysis?


title: Difference between binomial, negative binomial and Poisson regression

text: I am looking for some information about the difference between binomial, negative binomial and Poisson regression and for which situations are these regression best fitted. Are there any tests I can perform in SPSS that can tell me which of these regressions is the best for my situation? Also, how do I run a Poisson or negative binomial in SPSS, since there are no options such as I can see in the regression part? If you have any useful links I would appreciate it very much.


title: Correlating time series for 20 regions (SPSS)

text: I have a question to which I can't find an answer although I spent really awfully lot of time searching. I have time series data for about 20 regions of a country. Each time series covers 20 years. I want to measure to what extent the evolution of the variable x followed one pattern in the regions. I'm sure it did to some extent but I want to argue in my paper that they follow a similar pattern, so I need an exact measure. I was thinking about just correlating the data for the regions, but the cross- table I receive is not enough for me to conclude how high the similarities were. Any ideas? I saw the option autocorrelation and crosscorrelation but I guess that's not the thing I'm looking for. It seems that the forecasting option is also not the best way to do that. I want to mention the statistical analysis only at the very beginning of my paper - then I go into qualitative case studies to prove my point, so I don't need to propose a huge model for the the statistical data. I will appreciate your help.   * EDIT - My hypothesis is theory-driven. I assume x and y are causaly linked. Theoretically it makes sense. I've investigated the country qualitatively (historical analysis) and it seems to make sense. I've also checked it with country-level statistical data and there's a huge correlation between x and y, but I only have data for 22 years so I'm afraid the correlation is just a coincidence. I'm asking myself for example whether the correlation is not caused by changes in some, few, regions. That would invalidate by causal explanation. That's why I want to see if the data series are correlated - I've seen charts and all the regions seem to follow the pattern visible in the country-level data, but I want to be sure. I would love to check correlations of x and y on the regional level too! No idea how to do that (but that's less important, simply because I have some doubts whether the regional data I found on y is reliable and have no other data source).


title: R: product of dummy variables when using factor in a linear regression

text: I have data in a dataframe with the following columns: date, time, symbol, price I am attempting to run the following model in R                price ~ factor(time) + factor(symbol)*factor(time) + 0      where the first coefficient comes from a dummy variable of the time column in the dataframe and the second coefficient comes from the product of the dummy variable of the time column and the symbol column. I used lm() to attempt to do the model.                fit<-lm(price~factor(time) + factor(symbol)*factor(time) + 0, data=mydata)       However, that didn't work as it gets me some really weird coefficients that didn't say anything. I guessing as that was matching effects in some way and not just the straight forward product. I tried making a new variable by multiplying the product of the factors and then putting it in the regression but that didn't work as r told me that * is not in the possible operations for factor(). So, how can you multiply two factors in a linear regression? Your help is much appreciated, thank you in advance. edit: when I say weird I guess I should more say not what I want. I get coefficients for the firm, then coefficients for time then coefficients for firm and time. However, I only want the coefficients for firm and time to impact the model and ultimately I am really interested in what this does to the intraday effects which are approximated using the time dummy variables. If I am getting the interaction effects of the firm, time and then firm and time then this will impact the first coefficient as opposed to only firm and time affecting the first coefficient?


title: How is q value defined?

text: 1. From Applied Statistical Genetics with R by Foulkes (pFDR = positive false discovery rate): > The q-value is defined similarly as the greatest lower bound of the $pFDR$ > that can result from rejecting a test statistic $T$ and a sample $x$ with > $T(x)=t$ based on a rejection region $Γ$ in the set of nested rejection > regions. Formally, the q-value is written $$q(t) =\inf_{Γ :t∈Γ } [pFDR(Γ > )]$$ Since $pFDR$ is involved here, there must be several null hypotheses, and different null may have different testing rule. But here does it assume that there is one common testing rule $(T, Γ)$ for all the nulls?   2. From Wikipedia > The q-value is defined to be the FDR analogue of the p-value. **The q-value > of an individual hypothesis test** is the _maximum_ FDR at which the test > may be called significant. One approach is to directly estimate q-values > rather than fixing a level at which to control the FDR.     * Why are we talking about the q-value of **an individual** hypothesis test? Isn't the q-value for multiple nulls and their test rules?     * Should " _maximum_ FDR" be "minimum FDR" instead?     * What does "One approach is to directly estimate q-values rather than fixing a level at which to control the FDR" try to mean? What is "one approach" trying to solve?   3. Am I right that pFDR depends on and changes with the sample's distribution, so does q-value? Thanks and regards!


title: Why is the number of FPs less than reported positives in this case?

text: From page 17 of some slides, In multiple comparison, when "not all null hypotheses are true", why is the number of false positives $V$ less than the the number of reported positives $R$, i.e. $V < R$? There can still be true positive(s), right?


title: Use of different outputs in a meta-analysis?

text: I'm interested in performing a meta-analysis of instrumental variables analyses from a handful of different studies, asking for the same set of model results from all study investigators. However, at the moment I’m not familiar with Stata and the output `ivregress` generates. Is it possible to use separate outputs to perform a meta-analysis?


title: Good textbook for a one semester "probability and statistics" course for computer enginering 1st year

text: I am looking for a good textbook for a one semester "probability and statistics" course for computer engineering 1st year university students (e.g. no calculus). Can you suggest? Calculus should not be a prerequisite: Our program has probability and statistics on the first semester of first year. That means students will be taking this course alongside calculus. Hence they will not know calculus. This rules out a lot of excellent probability and statistics textbooks. No high level programming: Even though our students are computer engineering students, they will not know programming. If possible, the book should be related to the use cases of computer engineers. I am aware of similar "good book" questions in this site, but I think my situation is specific enough to ask another question.


title: Propagation of errors with median absolute deviation from the median?

text: Is there a theoretically-sound way to perform propagation of errors with robust statistics? I am trying to characterize the errors inherent in a measurement and propagate the uncertainty through calculations involving that estimate. In the past, I have assumed that my errors are normally distributed and then used uncertainty propagation. However, this is clearly a bad assumption (my data is symmetric but guaranteed to have significant outliers), and I am intrigued by using a measure like Median Absolute Deviation. Is it valid to treat MAD like a standard deviation? Is there a different measure I should be using for this application?


title: ANOVA on ranks vs. Mann-Whitney U test

text: What exactly is the difference between performing a t-test on data that has been ranked and performing a Mann-Whitney U test on the unranked data? I keep reading things like: For example, the Wikipedia page for ANOVA on Ranks states: > This rank-based procedure has been recommended as being robust to non-normal > errors, resistant to outliers, and highly efficient for many distributions. > It may result in a known statistic (e.g., in the two independent samples > layout ranking results in the Wilcoxon rank-sum / Mann–Whitney U test). The reason I am confused is that the Mann-Whitney U-test for differences in location assumes that samples come from identical distributions. And yet performing ANOVA on ranks has no such assumptions. Thanks!


title: Importance of intercept term in regression equation

text: What is the actual importance of an intercept term in a regression equation? If regressors include variables like age, height, weight and gender(1= male, 0= female), should the equation have an intercept term?


title: Is the Spearman rank correlation matrix positive definite or not?

text: I've always believed that the Spearman rank correlation matrix doesn't need to be positive semidefinite because the correlations are estimated pairwise so there's always a chance that it may not be. Recently, however, I have been looking deeper into this issue and I am confused, particularly after reading this: https://www.usenix.org/legacy/event/sysml07/tech/full_papers/sabato/sabato.pdf On page #4, Theorem 2 it says: _A Spearman rank correlation matrix is PSD_ **Proof:** _A Spearman correlation is a Pearson correlation applied to ranks. Therefore the Spearman rank correlation matrix is PSD_ So... can the Spearman correlation matrix be not positive semidefinite? Is there a counter-example people know? Or is it a misunderstanding to think that the Spearman correlation matrix could be non-positive definite?


title: Determining the significance of a rule

text: Let's say I establish the following rules from a data set: **Rule A:** If you exercise daily, you have a 70% chance of having a BMI under 28, based on 100 cases from the data. **Rule B:** If you eat fast food less than once per week, you have a 90% chance of having a BMI under 28, based on 10 cases from the data. The challenge for me is that Rule B seems to show a stronger correlation between a factor and a BMI under 28, but it is observed in much fewer cases than in Rule A, and therefore has a higher variance. How can I mathematically determine which rule is more significant?


title: Pooled time series regression in R

text: I've got data on mail volume sent by household for seven age groups, with 12 years of data for each age group. I originally ran a simple regression on each age group individually and realized I needed to dig deeper. My aim now is to pool the data (giving me 84 observations) and try to identify some period effects (or year effects, whichever you prefer). My pooled data are currently organized like this (PPHPY stands for Pieces per Household Per Year):               Age Group    Year   PPHPY     1            2001   127.62     1            2002   144.47     1            2003   111.70     1            2004   95.96     1            2005   96.46     1            2006   139.91     1            2007   85.52     1            2008   75.43     1            2009   109.34     1            2010   53.16     1            2011   64.09     1            2012   50.94             2            2001   176.48     2            2002   172.86     2            2003   137.79     .              .      .     .              .      .     .              .      .     7            2012   163.39      I first regressed PPHPY on year and year dummies (leaving the intercept as 0 to avoid perfect multicollinearity). This gave me period effects for the aggregated data (ie something like a period effect across all age groups, I think). This looked like the following:               > ## Generate YearDummy using factor()     >     > YearDummy <- factor(YearVar)     >     > ## Check to see that YearDummy is indeed a factor variable     >     > is.factor(YearDummy)     [1] TRUE     >     > ## (...+0) ensures intercept is left out and thus YearDummy1 remains in.         ## One or the other must be subtracted out to avoid perfect mutlicollinearity     >     > LSDVYear <- lm(PPHPY ~ YearVar + YearDummy + 0, data=maildatapooled)     > summary(LSDVYear)     Call:     lm(formula = PPHPY ~ YearVar + YearDummy + 0, data = maildatapooled)     Residuals:     Min 1Q Median 3Q Max     -99.658 -39.038 8.814 43.670 82.300     Coefficients: (1 not defined because of singularities)     Estimate Std. Error t value Pr(>|t|)     YearVar 5.743e-02 9.851e-03 5.830 1.45e-07 ***     YearDummy2001 1.099e+02 2.795e+01 3.930 0.000193 ***     YearDummy2002 1.209e+02 2.796e+01 4.324 4.85e-05 ***     YearDummy2003 7.791e+01 2.797e+01 2.786 0.006819 **     YearDummy2004 8.053e+01 2.797e+01 2.879 0.005251 **     YearDummy2005 6.887e+01 2.798e+01 2.461 0.016236 *     YearDummy2006 6.572e+01 2.799e+01 2.348 0.021618 *     YearDummy2007 5.975e+01 2.799e+01 2.134 0.036210 *     YearDummy2008 5.836e+01 2.800e+01 2.084 0.040696 *     YearDummy2009 4.119e+01 2.801e+01 1.471 0.145745     YearDummy2010 3.056e+01 2.801e+01 1.091 0.278990     YearDummy2011 1.472e+01 2.802e+01 0.525 0.600951     YearDummy2012 NA NA NA NA     ---     Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1     Residual standard error: 52.44 on 72 degrees of freedom     Multiple R-squared: 0.9316, Adjusted R-squared: 0.9202     F-statistic: 81.71 on 12 and 72 DF, p-value: < 2.2e-16      What I want, however, is to tease out period effects for each age group individually. This is what I'm not sure how to set up. I was hoping someone might help me devise some code in R that would kick out those period effects for _each_ of the seven age groups using the pooled data, as well as help me understand the problem conceptually. EDIT: I forgot to mention that I see I must include an interaction term involving the time dummies to allow the coefficients to vary across age groups. I'm just having difficulty constructing the proper interaction term and resulting regression equation. EDIT 2: I came up with two models and ran them. I felt like the question had evolved at this point and might merit a new post, which is can be found here.


title: Type 1 error Bayesian Stats

text: Is there anybody who can tell me or explain to me about finding the Type 1 error for a Bayesian clinical trial. Of a 100 patients I need 59 successful patients to exceed my posterior threshold. I have been able to work this through R but now am stuck trying to find the probability Type 1 error for these 59 successes. I know the value should be ~ 0.044, but has to how to find this value I am unsure. Any help please.


title: Tools for filling in missing values in a data set

text: I created an algorithm to fill missing values in a list of numbers. I have a list for example [2,4,5,3,6,8,4,NaN,5,NaN,4,6,4,2,4,6,NaN.....]. Think of it much much bigger. I am looking for tools where they can do the same thing. Fill those missing values. But not in random. In my case, I use a neural net to find the most likely number to put there. I have 90% correct with my algorithm and I want to compare it with other tools in the same data. Can you suggest any tool to try it?


title: How to check adequacy of a time series model using residual autocorrelations?

text: This was in an past exam question I came across. A first-order autoregressive model has been fitted to a time series of 50 observations giving $\hat\mu = 15$ and $\hat\alpha =0.6$. The first 12 residual autocorrelations were given and the estimated standard deviation of $r_k(\hat Z_t)=0.15 $ $ k=1,2...12.$ The question was to check the adequacy of the fitted model. The only adequacy test we were given was the Ljung-Box test which I think uses sample autocorrelations, not residual autocorrelations.


title: Algebra for data confidence

text: Very often, we use data which are derived from some measurements. These measurements usually have a confidence measure associated which tells how reliable or confident we are about the measure. For example, we often see some confidence intervals associated with various polls. I was wondering if there is a theory or algebra about combining multiple measurements and the resulting confidence measure of the aggregate. For example, if I measure X with a confidence interval +-x%, Y with confidence interval +-y%, what can I say about the confidence interval of X+Y? Similarly, what can we say about the confidence intervals for other operators? Is there an algebra for this?


title: In the glm function for logistic regression, where is the likelihood function stored? Is it in family?

text: I am currently trying to run a logistic regression on my own, using the functions `optim`, `nlm`, etc. However, I am somehow getting different coefficient values for the custom model as compared to when I run glm. I suspect that my likelihood function might be different, but cannot find where the logit likelihood function is stored in `glm`. Does anyone know? Thanks!


title: Nested random factors advice

text: I have four (or 6 depending on site) transects nested within 10 sites, with data repeated over 4 years, with 5 monthly visits. My response is count data and I am using a Poisson (and perhaps a Neg Bin) GLMM. Is it advisable to treat transects as a random effect? Is it perhaps better to look at the site x transect random effect interaction, or even to treat the number of transects within a site as an offset? I have tried running a model with the site x transect random effect interaction, but I receive a NaN message for the correlation between the slope and intercept, when I introduce covariates into the model (which are related to transects). Why might this be? Any advice is greatly appreciated. Thank you.


title: How to select no of rows from Data frame

text: I am reading CSV file and creating a data set Yourdata<-read.csv(‘C:/Work/mydat.csv’) This have five columns and nearly 100K rows and I want to select only few rows. Like, if I want to select only first 1000 rows then how can extract it in vector to do further processing on it.


title: Random effect with zero SD in LMM

text: In my mixed-effects model there are one fixed effect and two random effects (subject and time of measurement).               fit <- lmer(DV ~ group + (1|subject) + (1|time))      Now for some DVs variance and standard deviation for random effect time are 0. On pages 10/11 of his book (lme4: Mixed-effects modeling with R) Douglas Bates wrote: > An estimate of 0 for σ does not mean that there is no variation between the > groups. (...) The estimate, σ1 = 0, simply indicates that the level of > “between-group” variability is not sufficient to warrant incorporating > random effects in the model. The important point to take away from this > example is that we must allow for the estimates of variance components to be > zero. We describe such a model as being degenerate, in the sense that it > corresponds to a linear model in which we have removed the random effects > associated with Batch. Degenerate models can and do occur in practice. Even > when the final fitted model is not degenerate, we must allow for such models > when determining the parameter estimates through numerical optimization. So I do not entirely understand what that means for me. On the hand it sounds like it's not justified to include time as a random effect ("..not sufficient to warrant incorporating random effects in the model"), on the other hand it sounds as if it's OK for variance/SD to be zero (as long as interpreted correctly). Can you help me with the correct interpretation of the paragraph?


title: What is contingent in a contingency table?

text: The Merriam-Webster dictionary defines a _contingent_ event or situation as               1 : likely but not certain to happen : possible     2 : not logically necessary; especially : empirical     3 a : happening by chance or unforeseen causes       b : subject to chance or unseen effects : unpredictable       c : intended for use in circumstances not completely foreseen     4 : dependent on or conditioned by something else     5 : not necessitated : determined by free choice      Regarding current statistical terminology, since contingency tables are used to represent a large variety of data in many different situations, why do we call them "contingency" tables? In which of the five senses above is the word "contingency" used in this terminology?


title: Number of kernel evaluations in SVM training

text: What is the typical number of kernel evaluations (between two training vectors) performed during a (kernelized) Support Vector Machine (SVM) training? I am asking this question because I need to determine how much I need to optimize the kernel in order to have some hope to do the training in a reasonable amount of time. The current kernel calculation time is prohibitive (1+ hour), but there is room for improvement (by many orders of magnitude). Note: There are about 60k vectors in one class and 10k vectors in the other one.


title: Number of interactions in ANOVA with 4 independent variables

text: I'm reading an article that states in ANOVA with 4 independent variables, there will be 4 main effects and 11 interactions. I was wondering how the number, 11 interactions, is obtained?


title: Balance between randomization, interspersion and even initial conditions of experimental units

text: We are setting up a manipulative experiment in forest where we will test the effect of management on a number of health indicators.   We have 1 factor, 3 levels, 3 plots each replicated 3 times; experimental design is completely randomized. Our concerns and our two questions are about the best way to assign the factor levels to the plots. Plots are scattered in the forest but may be adjacent with one or more side shared between each one. Reducing plot size is not an option nor it is increasing replicates number. To assure independency of plots we are going to put a few constraints to the random assignment of levels to plots in order to assure that different factor levels be assigned to adjacent plots. **First question**   Do you think this compromise between randomization and interspersion would bias further analyses? We would like to guarantee that initial conditions of plots, as far as our response variables are concerned, are as much similar as possible. We have enough data on initial conditions and may test for equality of a response variable, for each combination of level assignment configuration to plots, before starting the experiment. Then we could pick up the levels configuration on plots that exceeds a significancy threshold (or simply maximizes the _p_ value), or minimizes the variance of the response variable explained by the factor level. **Second question**   Where is the correct balance between randomization and assuring even initial condition between the experimental units?


title: Does an unbalanced sample matter when doing logistic regression?

text: Okay, so I think I have a decent enough sample, taking into account the 20:1 rule of thumb: a fairly large sample (N=374) for a total of 7 candidate predictor variables. My problem is the following: whatever set of predictor variables I use, the classifications never get better than a specificity of 100% and a sensitivity of 0%. However unsatisfactory, this could actually be the best possible result, given the set of candidate predictor variables (from which I can't deviate). But, I couldn't help but think I could do better, so I noticed that the categories of the dependent variable were quite unevenly balanced, almost 4:1. Could a more balanced subsample improve classifications?


title: What if a binary variable is very unbalanced?

text: For example, if I have gender in the regression formula. My data has 90 Male and 10 Female. I know it is not good to have two groups so unbalanced. The estimate of gender will have relative big standard error. Is there any other downside for including it as explanatory variable? What if it is a dependent variable in logistic regression?


title: Logistic regression and ratio cases to non-cases in sample

text: I have a dataset with approx. 13,000 cases in total, and approx. 1,000 cases have the outcome I am interested in predicting (sentenced to imprisonment). I have run a logstic regression on the whole sample and the model is terrible at predicting 'imprisonment' although quite good at predicting 'no imprisonment'. I tried drawing a random sample of 'non-imprisonment' cases and joining them with the 'imprisonment' cases to achieve a more balanced sample. When I have a matched number of cases (i.e. 50% imprisonment 50% no imprisonment) I achieve quite good prediction of both outcomes. However, when I double the number of 'non-imprisonment' cases my prediction overall increases but the ability to prdict imprisonment in particular decreases. It is valid to use a sample to best predict imprisonment or am I somehow influencing the results? Does the fact that my initial model was poor at predicting imprisonment just mean I do not have a good model? (diagnostics generally indicate it is not great but OK). Thanks


title: How should decision tree splits be implemented when predicting continuous variables?

text: I'm actually writing an implementation of Random Forests but I believe the question is specific to decision trees (independent of RFs). So the context is that I'm creating a node in a decision tree and both the prediction and target variables are continuous. The node has a split threshold to partition data into two sets, and I create a new prediction for each subset based on the average target value in each set. Is this the correct approach? The reason I ask is that when predicting _binary_ variables I believe the typical (correct?) approach is to divide the data into 0 and 1 subsets without taking an average over the data rows in each subset. Subsequent splits will divide into finer grained subsets and taking an average at each split results subsequent splits (lower down the decision tree) operating on what are now continuous variables rather than binary variables (because we are operating on the residual error values instead of the original targets). Side question: Is the distinction between the two approaches (binary vs continuous) significant - or will they actually give identical results for a complete decision tree?


title: How to implement bootstrap efficiently in a matrix-based language?

text: I am working on a Monte Carlo study of bootstrapping in an AR(1) model for a homework assignment (I'm using Matlab). The goal is to say something about the empirical rejection probabilities of the bootstrap in this specific context. However, I've been running into some computational problems. Some context first. Suppose one has a regression $y_t = \gamma + \rho y_{t-1} + \varepsilon_t$, and wants to test the two-tailed hypothesis $\rho = \rho_0$ by the bootstrap. To do this, one first estimates the regression by, say, OLS, and then uses the obtained residuals $\hat \varepsilon_t$ and coefficients $\hat \gamma, \; \hat \rho$ for constructing the bootstrap samples. However, in this case, the bootstrap sample has to be constructed recursively. What is even worse for the computations, I have to repeat these calculations many times using different simulated data sets _and_ different values of $\rho$. (The reason is that this allows one to get some sense of the empirical level of the test at different values of $\rho$.) So far I have been doing it in the naive way of just estimating many regressions, generating the bootstrap samples, and so on. I was able to eliminate quite a few loops by using the `filter` function when generating bootstrap samples but still computational times are quite huge. I compared it with the built-in bootstrap in gretl, and it is many times slower. (I'm aware that gretl is written in C, but I sense that there is more to this than that.) Hence, my question would be: what general advice could you give to implement bootstrapping efficiently in a matrix based language? Are there some smart ways to avoid estimation of many OLS regressions and generation of the bootstrap samples? Although I'm using Matlab for the computations, but if you have some general advice on say, doing this with R, all the help would be very much appreciated. Also, since this is a homework assignment I have to do, using existing libraries is not really an option. :)


title: Test for independent but not identically distributed time-series?

text: So far, runs test or bds test are doing OK with i.i.d data but not so for independent but **not** identically distributed data set based on my practical experiences with them. Are there any other tests that have good properties at test these non-i.i.d data stream?


title: Solving logarithms using variables

text: LogQ= a+bLogP+cLogI+dLogPm b= -2.174 c= 0.461 d= 1.909 Determine price elasticity of demand, income elasticity, and cross price elasticity


title: Fit a moving average (MA) time series model to the data (R:stats::ar equivalent)

text: I am looking for some tools for **automatic fitting of moving average (MA) time series model** to my data in `R`. I know `R:stats::ar` method to fit an autoregressive time series model to the data (by default selecting the complexity by AIC) but I keep failing in finding its equivalent for MA model. There is `armaFit` method but it requires pre-defining an order of the MA model. I would like to have something that would do it automatically.


title: Significance of mediation

text: I have some trouble understanding the mediator analysis I just conducted. I used the Sobel macro of Preacher & Hayes (2004) for SPSS. The results look like this: DIRECT AND TOTAL EFFECTS                          Coeff   s.e.   t     Sig(two)       b(YX) ,1626 ,0526 3,0908 ,0035 b(MX) ,2954 ,0897 3,2948 ,0020 b(YM.X) ,2116 ,0959 2,2064 ,0327 b(YX.M) ,1263 ,0538 2,3462 ,0236 TOTAL EFFECT, b(YX), STANDARDIZED TO THE METRIC OF THE DIRECT EFFECT, b(YX.M)                      Coeff      b(YX) ,2310 INDIRECT EFFECT AND SIGNIFICANCE USING NORMAL DISTRIBUTION                       Value      s.e.  LL95CI    UL95CI   Z    Sig(two)       Effect ,0625 ,0352 -,0064 ,1315 1,7776 ,0755 So what I understand is that all paths got significant - which I thought means that there is mediation. But I also see that the Z-test and the bootstrap analysis tell me, that the effect is not significant. I don't really understand how all paths could become significant but then there is still no indirect effect? Would be right to assume that there might be an effect if the sample size was larger? Is there any other explanation?


title: Confused by Derivation of Regression Function

text: I just got a copy of _The Elements of Statistical Learning_ by Hastie, Tibshirani, and Friedman. In chapter 2 (Overview of Supervised Learning) section 4 (Statistical Decision Theory), he gives a derivation of the regression function. > Let $X \in \mathbb{R}^p$ denote a real valued random input vector, and > $Y\in\mathbb{R}$ a real valued random output variable, with joint > distribution $Pr(X,Y)$. We seek a function $f(X)$ for predicting $Y$ given > values of the input $X$. This theory requires a loss function $L(Y,f(X))$ > for penalizing errors in prediction, and by far the most common and > convenient is squared error loss: $L(Y,f(X))=(Y −f(X))^2$. This leads us to > a criterion for choosing $f$, > > $$\begin{align*} EPE(f) &= E(Y-f(X))^2 \\\ &= \int [y - f(x)]^2Pr(dx, > dy)\end{align*}$$ the expected (squared) prediction error. I completely understand the set up and motivation. My first confusion is: does he mean $E[(Y - f(x))]^2$ or $E[(Y - f(x))^2]$? Second, I have never seen the notation $Pr(dx,dy)$. Can someone who has explain its meaning to me? Is it just that $Pr(dx) = Pr(x)dx$? Alas my confusion does not end there, > By conditioning on $X$, we can write $EPE$ as $$\begin{align*}EPE(f) = > E_XE_{Y|X}([Y-f(X)]^2|X)\end{align*}$$ I am missing the connection between these two steps, and I am not familiar with the technical definition of "conditioning". Let me know if I can clarify anything! I think most of my confusion has arisen from unfamiliar notation; I am confident that, if someone can break this derivation down into plain English, I'll get it. Thanks stats.SE!


title: Cross validation and neural network performance

text: Is cross validation used for model selection only, or can I use cross validation to test the performance accuracy of a neural network?


title: Does thetaf() in the forecast package in R detect seasonality?

text: I had a question on the theta method in the `forecast` package in `R`. I tried using the `AirPassenger` example with `thetaf()`, but it just provided a trend line without seasonality in the forecast. Does `thetaf()` detect and forecast seasonality? This is the `R` code that I used:               library(forecast)     x.fit4 <- thetaf(AirPassengers, h=24)     plot(x.fit4)


title: Binary Logistic Regression – Extremely Low Log-Likelihood Value

text: I am a novice statistician using R. I am running a BLR to determine the probability of good or bad sites to build a structure. Each location has a composite "score" that was generated from index values created in a previous model. When I run the BLR, my intercept and coefficient seem out of whack. The intercept registers a value of -91.197, and the coefficient is 18.188. Likewise the log-likelihood does not seem right, as it is -2.1503e-07. Therefore when calculating the predicted probability from the Log Units and Raw Units, I am getting extremely low values like 1.024419e-166. Can anyone shed light on what the issue might be? Do I have a bad data set? Am I using the wrong statistical method? Any help is greatly appreciated.


title: Accuracy of aggregate vs. disaggregate forecasting

text: I've found a few interesting articles online on this topic, but none which appear to be too cut and dry. My question is coming up with an accurate predictive forecast based on forecasting individual component parts, then adding then up (or whatever function they comprise) ... vs. forecasting at an aggregate level (which will miss some trends but may have less variance). The example I'm working on: I'm trying to forecast phone and email contacts. The source of these contacts is customers ordering software products, and then have technical or payment issues, or what-have-you. I could, at the most bottom-up level, have a forecast for the contacts for each disparate software group (Software A has much different contact patterns than Software B). I could also split up forecasts between orders placed and contact rate (we may see contacts flat, but orders are increasing, and contact rate is decreasing). These all really point to the same question – I have all the granular data – is forecasting more accurate at the aggregate level (what are my total contacts this month?) or would forecasting the disaggregated parts, then adding them be better?


title: Required number of simulations for Monte Carlo analysis

text: My question is about the required number of simulations for Monte Carlo analysis method. As far as I see the required number of simulations for any allowed percentage error $E$ (e.g., 5) is $$ n = \left\\{\frac{100 \cdot z_c \cdot \text{std}(x)}{E \cdot \text{mean}(x)} \right\\}^2 , $$ where $\text{std}(x)$ is the standard deviation of the resulting sampling, and $z_c$ is the confidence level coefficient (e.g., for 95% it is 1.96). So in this way it is possible to check that the resulting mean and standard deviation of $n$ simulations represent actual mean and standard deviation with 95% confidence level. In my case I run the simualtion 7500 times, and compute moving means and standard deviations for each set of 100 sampling out of the 7500 simulations. The required number of simulation I obtain is always less than 100, but % error of mean and std compare to mean and std of entire results is not always less than 5%. In most cases the % error of mean is less than 5% but the error of std goes up to 30%. What is the best way to determine number of required simulation without know actual mean and std (in my case subjected outcome of simulation is normally distributed)? Thanks in advance for any help. * * * In order to have an idea about what may distribution of simulation results look like when iteration is run infinite number of times: Instead of using resulted mean and variance after n number of simulations, I've decided to find a fit function of resulted distribution, but here n has to fullfill allowed % error. I think on that way I can find more correct results on cumulative distrubution function that is related with e.g. 97.5%. Because when I compare the results of 400 and 7000 simulation, fit functions of distribution for both sampling looks like each other only curve of 2nd one is smoother. Also, hence model in MATLAB/Simulink is nonlinear, although generated input parameters are normal distributed, resulted histogram of simulations are not normal for that reason I used "generalized extreme value distribution", which is named as 'gev' in MATLAB. But still, I am quite not sure about this methodolgy, thanks for any command in advance


title: Negative Offset in Rate (Poisson or Negative Binomial) models

text: I have a dataset that contains:   * the counts of successes, $Y_i$   * the length of observation, $\text{length}_i$   * few predictors, $X_1, X_2, \text{etc}$ Since the counts are observed at different lengths I need a rate model. For the Poisson regression we achieve this by writing: $\log(\mu_i /\text{length}_i) = \beta_0 + \beta_1 X_1 +\beta_2X_2$ which is then rewritten as: $\log(\mu_i) = \beta_0 + \beta_1 X_1 +\beta_2X_2 + \log(\text{length}_i)$ The $\text{length}_i$ is the offset and it has a fixed coefficient of 1. In R we specify it using the `offset` method. Something like this: `model <- glm(Y ~ offset(log(LENGTH)) + x1 + x2, data = data, family = "poisson")` However, I get better prediction results when I specify the offset as having a coefficient -0.5. That is: `model <- glm(Y ~ offset(-0.5 * log(LENGTH)) + x1 + x2, data = data, family = "poisson")` This would mean that the mean is related to the predictor variables like: $\log(\mu_i * \sqrt{\text{length}_i}) = \beta_0 + \beta_1 X_1 +\beta_2X_2$ I have a few questions:   1. Why am I getting better results using negative coefficient?   2. What are the possible theoretical explanation for a negative offset?   3. Should I try some other model?


title: t-test or not, independent or not

text: I have searched, but I'm unable to understand: please help a beginner. 1) Scenario - 10 people (random sample) are asked 6 questions. The 6 questions are divided into 2 categories (i.e. 3 long words category questions and 3 short words category questions) which are the independent variable to test whether there is a difference between the long words category performance and the short words category questions (based on the word length effect on memory span for working memory). So the 10 subjects have individual scores based on long word performance and short word performance (all derived from being asked one questionnaire). I further go on and calculate the means for all the subjects: for example, the long word mean is 1.3 and the short word mean is 2.3. The standard deviations of the groups will also be calculated. 2) Question 1. To test if the two means from the same group with results based on different question categories of a single questionnaire administered in full to an individual, can I use t-test to establish if the means differ statistically (all material I have found talks about two different samples {normally as independent samples}, or a single sample but having undergone two separate tests {one before intervention then another after intervention, but the sample being referred to as dependent}? So I am not too sure about applying a t-test when the means are derived from one questionnaire and administered once to one group. 3) Question 2 - if I am to use a t test, should it be of independent or dependent?


title: Does this model make any sense?

text: I have fitted data of user activity per day to an arima(1,0,1)x(1,0,0). It's a nice fit. But does it make any sense?


title: Confidence band in correlogram

text: Could some one explain me, what is the confidence band in correlogram and if an autocorrelation coefficient is out of this line, what does it mean? For example in my Diagramm, lag 14 is significantly high, what does it mean? ![enter image description here](http://i.stack.imgur.com/pTqbG.png)


title: Random effect with 0 variance in GLMM

text: In this study, subjects are measured continuously over the day via electrocardiography (ECG). During the day, certain trigger events occur randomly. Once all the data is collected, the trigger events are examined and everytime a trigger event occured, the ECG is used to check for the presence of a certain medical condition. This means, that the number of observations per subject are different. In total there are $n=2240$ observations and $m=42$ subjects. The number of observations per subject range from 13 to 216. This study has already been analyzed and I'm trying to understand the results. To analyze the data a logistic regression with random intercept was choosen. I will use `glmer` from `lme4` in `R` to reproduce the results. Here $Y$ is a binary 0/1 variable and $X1, \ldots, X3$ are covariates of interest and `id` is the subject indicator. If I fit a model with just one variable I get this:               Generalized linear mixed model fit by maximum likelihood ['glmerMod']      Family: binomial ( logit )     Formula: Y ~ X1 + (1 | id)         Data: data                 AIC       BIC    logLik  deviance      1548.1535 1565.2962 -771.0768 1542.1535           Random effects:      Groups Name        Variance Std.Dev.      id     (Intercept) 0.02596  0.1611       Number of obs: 2240, groups: id, 42          Fixed effects:                   Estimate Std. Error z value Pr(>|z|)         (Intercept) -2.1235717  0.0756041 -28.088  < 2e-16 ***     X1           0.0007553  0.0001924   3.927 8.62e-05 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Correlation of Fixed Effects:        (Intr)     X1 -0.166      What we see is that the variance of the random subject effect is 0.026. It is small, but not equal to zero. If I add two more covariates, I get:               Generalized linear mixed model fit by maximum likelihood ['glmerMod']      Family: binomial ( logit )     Formula: Y ~ X1 + X2 + X3 + (1 | id)         Data: data                 AIC       BIC    logLik  deviance      1514.2474 1542.8186 -752.1237 1504.2474           Random effects:      Groups Name        Variance  Std.Dev.       id     (Intercept) 1.021e-11 3.195e-06     Number of obs: 2240, groups: id, 42          Fixed effects:                   Estimate Std. Error z value Pr(>|z|)         (Intercept) -2.1770583  0.0727357 -29.931  < 2e-16 ***     X1           0.0009218  0.0001914   4.816 1.46e-06 ***     X2          -0.3782134  0.0625888  -6.043 1.51e-09 ***     X3           0.0041024  0.0045297   0.906    0.365         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Correlation of Fixed Effects:        (Intr) X1     X2         X1 -0.229                   X2  0.281 -0.131            X3 -0.044  0.013  0.054      Now the variance it technically 0. How would I explain this? My interpretation is that due to the added covariates, all the variation is accounted for and there is no more subject specific variation left. Is this reasonable? Additional info about the covariates: $X1$ and $X2$ are trigger specific measurements (length and intensity).Their value usually changes for different trigger events. Hence they are different for each measurement even for the same subject. $X3$ is the age of the subject.


title: Non-normality of residuals in linear regression of very large sample in SPSS

text: I have a dataset of ~17,000 cases in SPSS 21 with which I am trying to run multiple linear regression. I have plotted the Studentised residuals against the unstandardised predicted values and also against each predictor that was included in the model and these plots indicate a fair degree of skewness and a few possible outliers. The outcome variable itself is highly skewed so I also fitted a model based on a Log10 transformation of this data which produced less non-normal (though now slightly skewed in the opposite direction) residuals. The data have been extensively cleaned and I'm quite confident that any outliers are not errors. I'm unsure as to how to proceed. I read somewhere that the normality of residuals is actually of minor importance so long as the other assumptions (validity and linearity) are met. Especially consdering the large sample size, do I need to worry too much about the non-normality of the residuals? And if I were to use the model based on the transformed data, how would I properly interpret the output? Thank you


title: Find specific 'y' value from GLM-fitted segmented multivariate regression

text: I loaded the data to R and fitted it using GLM function. > fit.glm = glm(y~ aX+bZ+cW) Then, I found the cuttinf-point using "segmented" tool of R. > o<-segmented(fit.glm,seg.Z=~X,psi=10) Now I have the cut-point and two different slope of X. > Call: segmented.glm(obj = fit.glm1, seg.Z = ~PTH, psi = 10) > > Meaningful coefficients of the linear terms: (Intercept) X Z W > > 19.43840 2.29574 0.08701 8.75784 > > Estimated Break-Point(s) psi1.PTH : 8.5 > > Degrees of Freedom: 324 Total (i.e. Null); 314 Residual Null Deviance: 17320 > Residual Deviance: 7645 AIC: 1971 However, I'm trying to get estimated y value of two linear regression models using the result of 'segmented'. Is it possible using R, or I have to substitute Z and W with mean values of Z & W, and calculate the y value myself?


title: Comparing error distributions to determine which model yields better results

text: I am a physicist who is a few years out of doing his last course in statistics, so I am hoping to get some advice when comparing some data I recently generated. I asked this question on the Mathematics Stack Exchange website, but received no answers, so I thought I would try here as well. The context is as follows. I have two slightly different theoretical models which I used to generate two sets of data (the numbers are of radiative strengths in multi-electron atoms for various transitions, but not exactly relevant for the mathematics). I also have a corresponding set of "accepted" or "literature" values for these transitions, and I would like to see which of my two models produces results that are closer overall/on average to the "accepted" model. My current metric of comparison is relative error for each transition (i.e. $\delta_{rel}=\frac{|T_{i,j}−A_j|}{A_j}$ where $T_{i,j}$ is my theoretical result for transition $j$ and theoretical model $i$ and $A_{j}$ is the accepted result for transition $j$). Seemed like the most obvious place to start. I am now stuck as to how to get some sort of meaning from these relative errors. Currently, I have just taken the arithmetic mean of the relative errors, but my intuition tells me this isn't a rigorous or correct way about this. From what I can remember, each relative error can be considered to be an independent random variable. If I then take the square of the difference in the numerator of my relative difference, then I believe I should get a $\chi^2$ distribution. Thus, I would have two samples from two theoretical $\chi^2$ distributions, one for each of my models. How should I go about comparing these two distributions? My hypothesis testing knowledge is a bit rusty, and I can't remember ever comparing two $\chi^2$ distributions. Any advice or indication of whether I am barking up the wrong tree would be appreciated.


title: Dependent proportion, and independent 5-point likert scale, what regression to use?

text: The dependent variable is a proportion as a measure of participation in an experiment (questions answered/questions asked). The independent variable is a 5 point likert scale that depicts the interest in the experiment. What sort of regression analysis is applicable here?


title: What are disadvantages of state-space models and Kalman Filter for time-series modelling?

text: Given all good properties of state-space models and KF, I wonder - what are **disadvantages** of state-space modelling and using Kalman Filter (or EKF, UKF or particle filter) for estimation? Over let's say conventional methodologies like ARIMA, VAR or ad-hoc/heuristic methods. Are they hard to calibrate? Are they complicated and hard to see how a change in a model's structure will affect predictions? Or, put another way - what are advantages of conventional ARIMA, VAR over state-space models? I can think only of _advantages_ of a state-space model:   1. It easily handles structural breaks, shifts, time-varying parameters of some static model - just make those parameters dynamic states of a state-space model and model will automatically adjust to any shifts in parameters;   2. It handles missing data very naturally, just do transition step of KF and don't do update step;   3. It allows to change on-a-fly parameters of a state-space model itself (covariances of noises and transition/observation matrices) so if your current observation came from a little different source than others - you can easily incorporate it into estimation without doing anything special;   4. Using above properties it allows easily handle irregular-spaced data: either change a model each time according to interval between observations or use regular interval and treat intervals without observations as missing data;   5. It allows to use data from different sources simultaneously in the same model to estimate one underlying quantity;   6. It allows to construct a model from several _interpretable_ unobservable dynamic components and estimate them;   7. Any ARIMA model can be represented in a state-space form, but only simple state-space models can be represented exactly in ARIMA form.


title: Log-normal model learning optimization

text: I am trying to model a variable (V) as a log normal distribution. $$ V = \mu \cdot\eta\,, \ \text{where } \eta \text{ is log normal} $$ Then $$ \ln V = ln(\mu) +\xi \\\\\\\\\text{where} \ \ \ \xi \ \text{ is }\ N(0,\sigma^2) $$ V is a function of the binary inputs $x_i, i\in(1,...,n)$ . I want to make the model depend on the inputs by assuming that $\mu=\Pi{f_1(x_i)}$ and $\sigma=\Pi{f_2(x_i)}$ where $f_1(x_i)=c_i$ if $x_i=1$ or 1 if $x_i=0$ and $f_2(x_i)=d_i$ if $x_i=1$ or 1 when $x_i=0$. My aim to obtain the best estimates for the $c_i,d_i $ values given an set of $(V_i,\\{x_j\\}_{j\in(1,n)})$ measurements. the log-likelihood in this case will become: $$ -\log L = \frac{1}{2}\left(\frac{\ln V - \Sigma \ln c_i}{\Pi{d_i}}\right)^2 + \log\sqrt{2\pi}\Pi{d_i} $$ What will be the best way to tackle this problem, to solve and obtain the best values for the $c_i$ and $d_i$ parameters? If I would assume that $\sigma$ is not a function of the inputs but rather a constant, I know I can solve the problem by a linear elastic-net regression, but when I take σ also as function of the inputs $x_i$ I am not sure which regression or method should I use to solve the problem? (how do I know the functional is convex?)


title: Neural network modeling sample size

text: I am preparing to conduct a study and was interested in using a neural network approach. I was wondering if there was a way to roughly work out the sample size that would be needed?


title: What are some data mining techniques for analyzing cause of disease

text: I have a dataset of 300 observations, of which 200 are normal, and the rest have the disease. I have the cognitive assessment scores of these 300 participants, and the assessment is divided into different sections: delusions, depression, anxiety, etc. I'm wondering what technique(s) would be good to look at the variation within the data. And whether a certain section's score would be a good predictor of the participant having the disease. I know of logistic regression, and with that I can construct a classification tree and an ROC curve. Are there any other techniques that I can use here? I'm open to different approaches (techniques for predicting, or just any other ways to look at the data in a meaningful way).


title: choosing $β_0$ and $β_1$ to minimize the residual sum of squares

text: I'm reading a book called An Introduction to Statistical Learning: with Applications in R, and I have a question in regards to the material inside. I understand that we can find the residual sum of squares by finding the difference between the $i^\text{th}$ observed response value and the $i^\text{th}$ response value that is predicted by our model. In other words the sum of $(y_i - \hat{y}_i)^2$ for every sample data point where $\hat{y}_i$ is the expected response from our model using constants $β_0$ and $β_1$. What I cannot understand is how to use calculus to find the minimizers $β_0$(this makes sense actually) and $β_1$ for the model. I was told that to come up with this equation that finds the minimizers we needed to find the minima of the RSS function using partial derivatives where both $β_0$ and $β_1$. I however have no clue how to do this, I am hoping that someone more experience could write it out for me. In the meantime I am going over partial derivatives to better grasp this. Here is a screenshot of the formula given in the book to find the minimizers: ![screenshot of the formula given in the book to find the minimizers](http://imgur.com/w388Dvq.jpg)


title: Log-likelihood (and AIC) of robust nlrob model differs from standard nls model

text: Comparing models generated by nlrob to ones generated by nls, I've noticed that even though the models might be nearly identical, the log-likelihood of the models is sometimes significantly different, e.g. (reproducible code below):               coefficients(fit.nls)             a         b         c       1.990388 -3.049477  1.019370      coefficients(fit.nlrob)             a         b         c       1.990650 -3.049251  1.019687       ![Fits plot](http://i.stack.imgur.com/KgnNJ.png) However:                logLik(fit.nls)     'log Lik.' 3.341953 (df=4)      logLik(fit.nlrob)     'log Lik.' 5.114307 (df=4)      Looking into the Loglik.nls code, this seems to arise from using the following line for obtaining the model's residuals:               res <- object$m$resid()      Using the same code for calculating log-likelihood but with the residuals() method instead, the obtained log-likelihood is much closer to the one generated for the nls model:               logLik.nlrob.method     [1] 3.340076      the nlrob documentation does relate to the different residuals, however I'm not clear if this has to do with the above: > residuals(.), by default type = "response", returns the residuals e_i, > defined above as e[i] = Y[i] - f(x[i], theta^). These differ from the > standardized or weighted residuals which, e.g., are assumed to be normally > distributed, and a version of which is returned in working.residuals > component. I'm not familiar with robust statistics, so I wonder:   1. Does the different log likelihood indeed reflects some deeper difference between the models? or is it just a technical difference arising from non-standard use of the resid() function?   2. Let's say I'd like to choose between nls and nlrob models based on the AIC criterion. for the sake of the question, I'll assume I don't know enough about the underlying model to say anything about the likelihood of outliers. Should I use the standard AIC/logLik functions for that? or should I use my own version, based on the residuals() function instead of the object\$m\$resid()? code (intended for interactive session):               library(robustbase)     library(ggplot2)          #generate data     set.seed(1357)     f <- function(x) 2*x*x-3*x+1     x=rnorm(50,0,2)     df <- data.frame(x=x,y=f(x)+rnorm(length(x),0,0.3))          #fit     fit.nls <- nls(y~a*I(x*x)+b*x+c,start=c(a=2,b=-3,c=1),data=df)     fit.nlrob <- nlrob(y~a*I(x*x)+b*x+c,start=c(a=2,b=-3,c=1),data=df)          coefficients(fit.nls)     coefficients(fit.nlrob)          #plot     s <- data.frame(x=seq(0,2,by=0.1))     ggplot(data=s,aes(x=x)) + geom_line(aes(y=f(x),color="f(x)")) + geom_line(aes(y=predict(fit.nls,s),color="nls")) + geom_line(aes(y=predict(fit.nlrob,s),color="nlrob"))          #log likelihood     logLik(fit.nls)     logLik(fit.nlrob)          #Let's compare the two versions of the residuals     res.obj <- fit.nlrob$m$resid()     res.method <- residuals(fit.nlrob)          N <- length(res.obj)          #all weights are 1 in both models, so the weight terms are cancelled     logLik.nlrob.obj <- -N * (log(2 * pi) + 1 - log(N) + log(sum(res.obj^2)))/2          logLik.nlrob.method <- -N * (log(2 * pi) + 1 - log(N) + log(sum(res.method^2)))/2                 #this one is identical to the original logLik.nls function     logLik.nlrob.obj     #and this one is very close to the fit.nls log likelihood     logLik.nlrob.method