title: Does the activation function of output layer differ during training and already trained network?

text: I'm creating and OCR app, and so far it seams to work. It's quite similar to example from Coursera - Machine learning course. Output layer of network has as many neurons as classes needed to recognize, 10 digits, 26 letter and alike. When the network is trained, it should basically output [1,0,0...0] for letter A, [0,1,0,...0] for letter B and so on. In the description of hidden layers (during Coursera lessons), it is said that they use sigmoid function as activation function, and I'm not sure if output layer also uses sigmoid function. So at the activation point, what differs when I train network, and when it is already trained? Is it the same activation function except I should be turning largest output value into 1 and rest to 0, or there is something I'm missing here? Explanation or links would be equally helpful. :)


title: Why do Lars and Glmnet give different solutions for the Lasso problem?

text: I want to better understand the R packages `Lars` and `Glmnet`, which are used to solve the Lasso problem: $$min_{(\beta_0 \beta) \in R^{p+1}} \left[\frac{1}{2N}\sum_{i=1}^{N}(y_i-\beta_0-x_i^T\beta)^2 + \lambda||\beta ||_{l_{1}} \right]$$ (for $p$ Variables and $N$ samples, see www.stanford.edu/~hastie/Papers/glmnet.pdf on page 3) Therefore, I applied them both on the same toy dataset. Unfortunately, the two methods do not give the same solutions for the same data input. Does anybody have an idea where the difference comes from? I obtained the results as follows: After generating some data (8 samples, 12 features, Toeplitz design, everything centered), I computed the whole Lasso path using Lars. Then, I ran Glmnet using the sequence of lambdas computed by Lars (multiplied by 0.5) and hoped to obtain the same solution, but I did not. [edit:Please ignore the following paragraph, the diffierent intercepts were caused by a bug in the code, which I fixed now, sorry for that. The rest of the question stays relevant, as fixing the bug did not solve the problem adressed.]   When I saw the (non-identical) outputs, I recognized that glmnet computed an intercept for each lambda, while lars did not. So I subtracted by hand the intercept of some glmnet solution from Y and recentered Y. Then I computed the Lasso path again using Lars, and then with glmnet using the new lambdas from Lars. One can see that the solutions are similar. But how can I explain the differences? Please find my code below. There is a related question here: GLMNET or LARS for computing LASSO solutions? , but it does not contain the answer to my question.               ### Load packages     library(lars)     library(glmnet)     library(MASS)          ###Set parameters      nbFeatures<-12     nbSamples<-8     nbRelevantIndices<-3     snr<-1     nbLambdas<-10     nlambda<-nbLambdas               ### Create data, not really important      sigma<-matrix(0,nbFeatures,nbFeatures)     for (i in (1:nbFeatures))       for (j in (1:nbFeatures))         sigma[i,j]<-0.99^(abs(i-j))     X<-mvrnorm(n = nbSamples, rep(0, nbFeatures), sigma, tol = 1e-6, empirical = FALSE)     relevantIndices <- sample(1:nbFeatures, nbRelevantIndices, replace=F)     X<-scale(X) #bug fixed, replaced 'scale(X)' by 'X<-scale(X)'     beta <- rep(0, times=nbFeatures)     beta[relevantIndices] <- runif(nbRelevantIndices,0,1)     epsilon<-matrix(rnorm(nbSamples),nbSamples,1)     simulatedSnr<-(norm(X %*% beta,type="F"))/(norm(epsilon,type="F"))     epsilon<-epsilon*(simulatedSnr/snr)     Y <- X %*% beta + epsilon;     Y<-scale(Y) #bug fixed, replaced 'scale(Y)' by 'Y<-scale(Y)'          ### launch lars      la <- lars(X,Y,intercept=TRUE, max.steps=1000, use.Gram=FALSE)     coLars <- as.matrix(coef(la,,mode="lambda"))     print(round(coLars,4))                    ### Launch glmnet with lambda=1/2*lambda_lars      glm2 <- glmnet(X,Y,family="gaussian",lambda=0.5*la$lambda,thresh=1e-16)     coGlm2 <- as.matrix(t(coef(glm2,mode="lambda")))     print(round(coGlm2,4))          ##########################     #edit: The rest of the program is not relevant any more, as the issue of the intercept was due to a bug.         ##########################          ### Remove intercept computed by glmnet for some lambda (adjusted by hand, depending on output)      Y<-Y--0.0863     Y<-Y/norm(Y,"F")          ### Launch Lars      la <- lars(X,Y,intercept=TRUE, max.steps=1000, use.Gram=FALSE)     coLars <- as.matrix(coef(la,,mode="lambda"))     print(round(coLars,4))          ### Launch Glmnet with lambda=1/2*lambda_lars      glm2 <- glmnet(X,Y,family="gaussian",lambda=0.5*la$lambda,thresh=1e-16)     coGlm2 <- as.matrix(t(coef(glm2,mode="lambda")))     print(round(coGlm2,4))


title: Find the test statistic of $H_0:m_1=\lambda m_2$

text: Three independent samples $(X_1,X_2,...,X_n)(Y_1,Y_2,...,Y_n)$ and $(Z_1,Z_2,...,Z_n)$ are drawn from normal population, each having the same unknown variance $\sigma^2$ but, $E(X_i)=m_1,E(Y_i)=m_2,E(Z_i)=m_1-m_2,\forall i$,show that the test of $H_0:m_1=\lambda m_2$ can be carried out by the test statistic $$T=\frac{[(2-\lambda)\bar X+(1-2\lambda)\bar Y+(1+\lambda)\bar Z]\sqrt n}{s\sqrt{6(\lambda ^2-\lambda +1)}}$$ where $s^2$ is the pooled sample variance , $\lambda$ is known constant and $T$~$t_{3n-2}$ , under $H_0$. First I think $m_1=\lambda m_2$ as $m_1-\lambda m_2=0$ and I can define $T=E(X_i)-\lambda E(Y_i)$ .But since there also given $E(Z_i)=m_1-m_2$ It should be used.Please help me.


title: Simulation of binomial random variables

text: I am getting started in R. I need help with a command to: "Simulate 10000 binom(20,0.3)-RV's". I would do it like:               binsim <- rbinom(10000, 20, 0.3)      Question cont.: "Let X be a binomial(20,0.3)-RV. Use simulated numbers to estimate: P(X<=5) and P(X=5)" Okay for the first I would use: `pbinom(???)` and for the second `dbinom(?????)` How do I connect the simulation results from `binsim and the X-RV's to get the required probabilities? Could I do it with `length(binsim(binsim<=5)/length(binsim)`?


title: Similarities and dissimilarities between optimal scaling procedures in CATPCA and MCA

text: I have seen that multiple correspondence analysis (MCA) can also be seen as a generalization of principal component analysis when the variables to be analyzed are categorical instead of quantitative. In SPSS I have found both CATPCA (Categorical Principal Components Analysis) and MCA under the optimal scaling procedure. I am quite interested in the alternating least square method and optimal scaling method that these procedures use. How do the optimal scaling of these procedures differ? Can someone give a mathematical illustration of it? It would be great if someone comes with some illustrative equations about optimal scaling, CATPCA and MCA. Can we get Cronbach's alpha for these two methods? Is the inertia more accurate measure than eigenvalues for MCA? Could you tell me the differences, possibly mathematically?


title: Correlation analysis: Spearman's $\rho$ is significant, but Pearson's $r$ is not - with likert scale

text: I am doing a correlation analysis testing summated likert scales. However, one of the correlation analysis shows an unsignificant pearsons r (p = 0.068), BUT spearman's rho is 0.044. What should I conclude from that? I know that likert scales are ordinal, but to run the correlation analysis I am assuming interval scales. I run regression analysis on the significant values, and the above mentioned variable has an ANOVA with p-value of 0.136 meaning no significant relationship between the DV and IV.


title: Not sure how to carry out this two-way repeated measures ANOVA

text: I have a scenario where both factors are within-subject. One is time (5 levels) and one is condition (2 levels). The time point of 0 was only measured once in each subject, and is a 'resting value' that is to be used as the time point of 0 for both conditions. I cannot find a way of dealing with this accept for entering in the value twice into the data set as t=0 for both conditions. I feel like this is not a valid approach but I don't know what else there is to do. I could run two separate one way ANOVAs but then I couldn't see an interaction between the two conditions. I use SPSS. UPDATE: Thanks for the responses. The nature of the methodology/study design includes assuming rest is the same for both conditions (to reduce tissue sampling required). Therefore, entering the rest measures twice in the two-way ANOVA is valid and is the nature of the methodology in this case, albeit a limitation.


title: What is the distribution of the 'achieved' $R^2$?

text: I am interested in the distribution of/performing inference on the 'achieved' $R^2$ coefficient in multiple linear regression. Suppose that $y\sim x\beta + \epsilon$ with $\epsilon \sim \mathcal{N}\left(0,\sigma^2\right)$, where $x$ is a $p$-vector. You observe $x_1,x_2,\ldots,x_n$ and corresponding $y_1,y_2,\ldots,y_n$, with independent errors. You perform linear regression to get $\hat{\beta}$, and then compute the _sample_ $R^2$ in the usual way: $R^2 = 1 - \frac{SS_\mbox{err}}{SS_{\mbox{tot}}}$. Nothing new here. There are known methods for hypothesis testing the _population_ analogue of $R^2$, as well as computing confidence intervals on it. (I am thinking of Lee's paper, and followups by Algina, O'Brien, _inter alios_ ). The population analogue is defined in terms of the (unknown) true regression coefficients $\beta$ and the true variance $\sigma^2$. I am interested in the _achieved_ $R^2$ which is the amount of variance explained by $\hat{\beta}$. Conditional on $x$ _and on_ $\hat{\beta}$, I would define it as $$ R^2_{\mbox{ach}} = 1 - \frac{\mbox{E}\left[\left(y - \hat{\beta}x\right)^2\right]}{\mbox{E}\left[\left(y-\bar{y}\right)^2\right]}. $$ Clearly this is less than the _population_ $R^2$, because the true $\beta$ minimizes $\mbox{E}\left[\left(y - \beta x\right)^2\right]$. I would guess it is less than the _sample_ $R^2$ because that (usually) overestimates the population $R^2$. I am having problems thinking about this quantity because it is both unobserved _and_ random ($\hat{\beta}$ is a random variable.) I am guesing that, up to some transforms, I can represent this as a non-central chi-square, abusing $\hat{\beta}\sim\mathcal{N}\left(\beta,\sigma^2\left(X^{\top}X\right)^{-1}\right)$, and noting that $\mbox{Var}\left[y - \hat{\beta}x\right] = \sigma^2 + \mbox{Var}\left[(\hat{\beta} - \beta)x\right]$, but I'm being dense and lazy. (Also, treating $x$ as fixed will later be generalized to $x$ being random.) Is this is a well-known problem? For predicting future performance of a linear model, it would seem to be a more important quantity than the population $R^2$, for example.


title: How can invariances be taken advantage of when selecting a kernel?

text: Some answers to other questions on kernel selection suggest using invariances known about the problem domain, and choosing a kernel that preserves such invariances. Can you give an explicit example of it's done?


title: ordered probit regression with only categorical variables

text: I would like to know if it makes sense to run an ordered probit regression (dependent variable is ordinal with three outcomes) with only categorical explanatory variables (some are dichotomous e.g. gender; others are categories that have been entered as dummy variables e.g age). Thanks


title: Suggest a method to assign weights

text: I am trying to prioritize customer tickets (issues/change requests). I have identified a few parameters to prioritize the tickets in the queue. The parameters are:   1. Criticality of the issue (minor - 100 points, major - 200 points and critical - 500 points).    2. Type of customer - Silver - 100, Gold - 200, Platinum - 500 points (based on the fees paid for support)    3. Turnaround time/efforts needed to fix the issue (to reduce wait time) etc.  The prioritization score would be a sum of points obtained for all the parameters   1. The numbers mentioned above for minor, major and critical (similarly for silver, gold and platinum) are just assumed numbers, not derived based on any model. How can I arrive at these scores mathematically.   2. Also, please suggest a method to assign weights to each of the parameters to decide the contribution of each parameter to the final score.


title: When to average samples in order to avoid pseudoreplication?

text: My question is when to average samples or when it is allowed to apply a mixed effects model. I have the following constellation: there is a cell-line experiment conducted in a 96-wellplate, meaning that we have cell material from the same source but incubate samples in 96 individual wells on one large plate (the 96-wellplate). This plate is repeated every day for seven days, so we have for each day 96 wells of each plate thus n=672 observations. On every plate, there are 12 treatments (including 4 controls), which are repeated 8 times. Since for every day the preparation of treatments were mixed again. I ask now whether the 8 observations per treatment per day need to be averaged, so I end up with only 7 observations per treatment or if I can use mixed effect modelling and use day as a random factor. In literature I found all possible combinations of averaging only, averaging and applying mixed effect models and also only mixed effect models... the data is "mainly" normal distributed (one treatment is a bit variable...) hence no need for glmms .


title: Contradicting interaction significance with graphical visualisation in glmer

text: I'm investigating environmental effects (wind) on acoustic receiver detection probability for two types of transmitters using a binomial glmer. While my model analysis indicates that there's a significant effect between wind speed and transmitter type, graphical visualisation does not confirm this. If I'm correct, an interaction should demonstrate different regression slopes.               m1 <- glmer(cbind(df$Valid.detections, df$False.Detections) ~ Transmitter.depth +                     Receiver.depth + Water.temperature + Wind.speed + Transmitter +                     Distance + Habitat + Replicate + (1 | Day) + (Distance | SUR.ID) + Transmitter:Distance + Transmitter.depth:Habitat +                     Receiver.depth:Habitat + Wind.speed:Transmitter, data=df, family=binomial(link=logit))      The model summary is as follows:               Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']      Family: binomial  ( logit )     Formula: cbind(df$Valid.detections, df$False.Detections) ~ Transmitter.depth +           Receiver.depth + Water.temperature + Wind.speed + Transmitter +           Distance + Habitat + Replicate + (1 | Day) + (Distance |           SUR.ID) + Transmitter:Distance + Transmitter.depth:Habitat +      Receiver.depth:Habitat + Wind.speed:Transmitter        Data: df               AIC      BIC   logLik deviance df.resid        3941.9   4043.8  -1953.9   3907.9     2943           Scaled residuals:          Min      1Q  Median      3Q     Max      -9.4911  0.0000  0.0000  0.5666  1.9143           Random effects:      Groups Name        Variance Std.Dev. Corr      SUR.ID (Intercept)  0.33414 0.5781                    Distance     0.09469 0.3077   1.00      Day    (Intercept) 15.96629 3.9958            Number of obs: 2960, groups:  SUR.ID, 20 Day, 6          Fixed effects:                                     Estimate Std. Error z value Pr(>|z|)         (Intercept)                      3.20222    2.84984   1.124  0.26116         Transmitter.depth               -0.35015    0.11794  -2.969  0.00299 **      Receiver.depth                  -0.57331    0.51919  -1.104  0.26949         Water.temperature               -0.26595    0.11861  -2.242  0.02495 *       Wind.speed                       1.31735    1.50457   0.876  0.38127         TransmitterPT-04                -0.68854    0.08016  -8.590  < 2e-16 ***     Distance                        -0.39547    0.09228  -4.286 1.82e-05 ***     HabitatFinger                   -0.23746    3.57783  -0.066  0.94708         Replicate2                      -0.21559    0.08009  -2.692  0.00710 **      TransmitterPT-04:Distance       -0.27874    0.08426  -3.308  0.00094 ***     Transmitter.depth:HabitatFinger  0.73965    0.28612   2.585  0.00973 **      Receiver.depth:HabitatFinger     3.02083    0.74546   4.052 5.07e-05 ***     Wind.speed:TransmitterPT-04     -0.15540    0.06572  -2.364  0.01806 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Correlation of Fixed Effects:                 (Intr) Trnsm. Rcvr.d Wtr.tm Wnd.sp TrPT-04 Distnc HbttFn Rplct2 TPT-04: Tr.:HF Rc.:HF     Trnsmttr.dp -0.024                                                                                Recevr.dpth -0.120 -0.267                                                                             Watr.tmprtr  0.019 -0.159  0.007                                                                      Wind.speed   0.130  0.073 -0.974  0.040                                                               TrnsmtPT-04 -0.015  0.027  0.020 -0.018 -0.024                                                        Distance     0.022 -0.080  0.151 -0.052 -0.141 -0.164                                                 HabitatFngr -0.813  0.010  0.241 -0.025 -0.253  0.009   0.029                                         Replicate2  -0.067  0.033  0.377 -0.293 -0.394  0.010   0.085  0.103                                  TrnsPT-04:D -0.006  0.043 -0.007 -0.050 -0.003  0.516  -0.373  0.004  0.006                           Trnsmtt.:HF  0.017 -0.352  0.021  0.055  0.049  0.026  -0.142  0.031 -0.088  0.025                    Rcvr.dpt:HF  0.103  0.189 -0.830  0.051  0.817 -0.036  -0.143 -0.224 -0.385 -0.003  -0.229            Wnd.:TPT-04 -0.002  0.026 -0.015  0.003 -0.009  0.176  -0.114 -0.002  0.016  0.306  -0.008  0.014      ![enter image description here](http://i.stack.imgur.com/DFxyQ.png) A side question: I noticed a strong negative correlation between the intercept and a dichotome categorical predictor. I wonder if this causes any problems for my data analysis. All the covariates are centered and scaled for numerical stability during modelling.


title: Threshold in precision/recall curve

text: While I was reading Torgo's Data Mining with R, I found that the description of precision/recall curve was different compared with other approaches. Usually, these curves are based on a threshold that determines which probability value is good enough to decide when an event has occurred, so we can classify future events depending on that value. However, Torgo's description is as follows: > Precision/recall (PR) curves are visual representations of the performance > of a model in terms of the precision and recall statistics. The curves are > ob- tained by proper interpolation of the values of the statistics at > different working points. These working points can be given by different > cut-off limits on a ranking of the class of interest provided by the model. > In our case this would correspond to different effort limits applied to the > outlier ranking produced by the models. Iterating over different limits > (i.e., inspect less or more reports), we get different values of precision > and recall. PR curves allow this type of analysis. The application the author has in mind is that of a fraud detection problem in which we have a classification task resulting in values `fraud`, `unknown` and `ok`. We would like to output probabilities, rank them, select the first $k$ reports and be able to inspect them. Is this an alternative measure of threshold in precision/recall curves? I think it is assuming that probabilities below 0.5 are to be classified as `ok`, 0.5 is equivalent to `unknown` and above 0.5 means `fraud`. Is that a correct assumption to make? Thanks a lot!


title: Granger Causality test with dummy variables

text: I intend to assess Granger Causality between three endogenous variables, where one of these variables is a dummy variable, indicating specific events during the continuous time frame. I was wondering whether normal Granger Causality tests can be used in the this case. I presume that I can, since the GC test is but an F-test of joint significance of other endogenous variables' lagged coefficients. So, if anything, the VAR model might not work with dummies. Anyhow, I was hoping I could get a second opinion on this: Can I safely use dummies in a VAR and subsequent Granger Causality test?


title: Compare differences or ratios of factors by time

text: I observe a process in which particles could be detected on a start point, and sometimes cross a border. So I count the number of particles during the days. Which non-parametric method would you suggest to clarify if the difference or ratio of start and cross border factors is different between two given days?                                 Days of observation     Particles         2013-07-22 2013-07-23 2013-07-24 2013-07-25 2013-07-26 2013-07-27 2013-07-28 2013-07-29 2013-07-30 2013-07-31     start             69         297        53         31         34         18         32         255        78         12     crossborder       14         231        45         24         26          6         47         193        66          4


title: Robust Residual standard error (in R)

text: I have a question regarding to the concept of robust standard errors. What I found about that topic is, that one can estimate the robust standard error for regression coefficients to eliminate problems with heteroscedasticity (when one wants to interpret a model). I want to know if there is a way not only to determine robust standard errors of coefficients but also of the standard error of the overall regression (residual standard error). When its possible, how can I calculate such a value in general? Because I'm using R its also interesting for me if there is a `R-function` for this problem (I only know the `sandwich-package` for the normal robust SE of the coefficients). Thanks.


title: Adjusted $R^2$ & F test are not shown in regression with robust standard errors in Stata

text: The adjusted $R^2$ is not shown when a regression with robust standard errors is calculated in Stata. This is surprising to me since the value of the $R^2$ is unaffected in regressions with robust standard errors. Is there any statistical reason for not quoting the adjusted $R^2$ when using robust standard errors in regression? Furthemore if I add more variables the F test disappears.(e.g. with 9 variables it shows but not with 13), Is this for the same reason? How can we report such results and deal with this issue?


title: Adjust estimate precision by sample size

text: I wish to find comparable SDs/Confidence Intervals for two parameters. One parameter has been measured by only one (new) study that reports the standard error. Here, _t_ value, sample size and SD are provided. The other (a correlation coefficient) has been measured by 10 (older) studies, in some as a mean across subjects, in some across all trials (ignoring subjects). The only values available for all studies are n of participants and the mean correlation coefficient; many don't report sufficient information to compute comparable SDs. I believe that all in all, _r_ is roughly comparable across the 10 studies in this set. I don't think I can directly compare the sample SD between subjects in set 1 with the between-study SD in set 2. How can I find comparable CIs/SDs for the two values while somehow controlling for the low sample size in the first study? I'm not trying to compare the two values using a test (they're physiologically incommensurable), I simply want to look at them on a similar scale. If I cannot come up with a proper quantification, I would simply be informal and focus on one SD for the 10-study set and two SDs for the one-study set, to reflect that one of them has been estimated "a lot more" precisely than the other.


title: The effective sample size from WinBUGS results

text: I am running a Bayesian regression model by WinGUBS via R2WinBUGS package in R. Everything looks fine except for one parameter:               nlssim$summary[37,][c(1,2,8,9)]         mean           sd         Rhat        n.eff      3.054326e-05 9.523965e-06 1.000000e+00 1.000000e+00       The number of effective size is 1! It seems to indicate that the autocorrelation between samples are extremely large, but the trace plot looks all fine. And I then try the effectiveSize function in code package and the result is 4590. I am curious how WinBUGS calculate the n.eff statistic, and if n.eff=1 is symptom of some of my mistakes? I run 3 chains in parallel, and each chain has 60000 iterations. n.burnin=200, n.thin=30 Thank you very much in advance!


title: Shall I use weekly or monthly data for forecast?

text: I seem not to find this in any textbooks. So I post these questions.   1. Is monthly data better than weekly data for forecasting?   2. Can there be seasonality in weekly data? Most software/methods don't seem to find seasonality in forecasting data.   3. Is there a way to aggregate weekly data to monthly data?   4. How do methods like ARIMA/Exponential Smoothing handle seasonality with weekly data?


title: What is the $\mu^2$ squared effect size?

text: I was recently looking at a paper in the journal _Psychological Science_ and came across this: _F_ (1, 71) = 4.5, _p_ = .037, $\mu^2$ = .06 _F_ (1, 71) = 0.08, _p_ = .78, $\mu^2$ = .001 I was wondering what the $\mu^2$ is in the above. Typically in APA the third thing should be either the MSE or it should be a standardized effect size (or you should have all 4). I'm guessing it's a standardized effect size of some sort but I'm not familiar with it and searching the net has turned up nothing. The actual effect, as near as I can tell from the graph, is about 12 for the first one. Is this an effect size I haven't heard of yet or a typo in the article? Farrelly, D., Slater, R., Elliott, H. R., Walden, H. R. and Wetherell, M. A. (2013) Competitors Who Choose to Be Red Have Higher Testosterone Levels. _Psychological Science_ , DOI:10.1177/0956797613482945 Here's a screen shot of the text (p.2) ![enter image description here](http://i.stack.imgur.com/h5BtM.png)


title: Weighted Wilcoxon ranksum test

text: I am using Stata for a survival analysis project involving inverse probability weighting (IPW). The question has arisen as to how to analyze weighted continuous data between two groups with a Wilcoxon rank sum test? Is this easily implemented in any other package?


title: AIC criterion: definition

text: I have two questions regarding the AIC criterion : AIC=$2k-2ln(l)$ Where does the number 2 comes from? As we usually minimize it why don't we consider only : $k-ln(l)$. (Maybe I am missing something.) Second question : Can we consider that AIC selection follow the same "philosophy" as lasso/ridge selection ? Indeed, if we consider that $0^{0}=0$ we would then have: AIC selection: max[$ ln(l)-\lambda\sum|\beta|^{0}$] (with $\lambda=1$) Lasso selection: max[$ ln(l)-\lambda\sum|\beta|^{1}$] Ridge selection: max[$ ln(l)-\lambda\sum|\beta|^{2}$]


title: Stuck on a proposition

text: Suppose four numbers $\\{a,b,c,d\\}$ where $a$ and $c$ random variables from a continuous distribution with support on $\mathbb{R}$. Does $b\neq d$ imply $|a−b|−|c−d|+a−c\neq 0$ almost surely? I really want to believe this proposition is true but that is typically not enough to make it so. # Edit: MickMack's counter example made me realize I need to refine the proposition to make it: Does $|b|\neq|d|$ imply $|a−b|−|c−d|+a−c\neq 0$?


title: Regression analysis if one of dependent variable is almost constant

text: While doing multi-variate regression, i am encountering cases like below image, where for one of dimension, apart from couple of data points all values are same. Cause of this, model is going haywire. ![enter image description here](http://i.stack.imgur.com/G6NjO.jpg) What's right way to handle these cases ? Also these outlier points are not explained by other dimensions also.                1. remove these points as outlier ?      2. ?


title: What is the difference between GLM and splines?

text: Suppose we want to predict $Y$ given the following $X$ observations:               x = c(abs(rnorm(2500, 0.1, 0.25)), abs(rnorm(2500, 0, 0.05)))     y = (x^0.35) + rnorm(length(x), 0, 0.25)     x = c(x, -x)     y = c(y, -y)      Clearly we have an exponential relationship between the predictor $X$ and $Y$. It's more obvious looking at a local regression plot:               regrplot = function(x, y, main.par="", regr.par=T, lowess.par=T, xlab.par=NULL, ylab.par=NULL)     {       plot(x, y, cex=0.5, col="red", main=main.par, xlab=xlab.par, ylab=ylab.par, pch=".")       if( regr.par )         abline(lm(y ~ x), col="blue", lwd=1)       if( lowess.par )         lines(lowess(x, y), lwd=1, col="darkgreen")     }      ![enter image description here](http://i.stack.imgur.com/wVPxy.png) This is a replica of an empirical distribution I have to model and predict. What are the difference in using splines vs. GLM? Which GLM regression better predicts this kind of distributions? From what I understand splines are better suited to study empirical distributions, because of the flexibility of the polynomials based on a guided number of knots. GLMs seems to be better to study theoretical distributions. Of course, the drawback with splines is overfitting, especially for prediction.


title: Are eigenvectors obtained in Kernel PCA orthogonal?

text: As Kernel PCA is the same as PCA in higher dimension space, shouldn't the eigenvectors obtained be orthogonal? Suppose, I have $n$ data points and let $a$ and $b$ be two eigenvectors of covariance matrix of mapped data, and $\alpha \in \mathbb{R}^n$ and $\beta \in \mathbb{R}^n$ be the corresponding eigenvectors of Kernel PCA problem (obtained by doing eigenanalysis of kernel matrix $K$). Then $a$ and $b$ can be written as linear combination of input data, i.e. $$a = \sum_{i=1}^n \alpha_i \phi(x_i)$$ and $$b = \sum_{j=1}^n \beta_j \phi(x_j).$$ Their inner product is: $$\langle a,b \rangle \,=\, \left\langle \sum_{i=1}^n \alpha_i \phi(x_i) ,\sum_{j=1}^n \beta_j \phi(x_j)\right\rangle \, = \, \sum_{i=1}^n\sum_{j=1}^n \alpha_i \beta_j K_{ij}.$$ Why is this value zero?


title: Comparing estimators of location of the Cauchy distribution

text: I'm comparing the following 4 estimators of location of the Cauchy distribution: Let $x_{1},..x_{n}$ be observations and $l$ be the log likelihood function. $x=median(x_{1},..x_{n})$, $y=x+\frac{l'(x)}{n/2}$, $z=mean(x_{1},..x_{n})$ and $s=MLE$ I did simulations in R of the mean squared error and probability coverages. I found that the mean was a very bad estimator. The best estimator was $y$ in terms of giving the smallest MSE and the greatest probability coverage. $x$ and $z$ were pretty good estimators as well and in fact when I increased $n$ large enough $x$ and $s$ gave the same results. I know the reason $x$ is a bad estimator. What are the mathematical reasons for my other findings? Thanks


title: How to combine the information of two probabilities?

text: I have a segmentation problem that using two different models I can assign a probability to each pixel of the image to be a background part. Two models produce slightly different results. Some pixel may have zero probability using one model, but have some value using the other model. But the idea is try to combine the outputs of two models. The question is how to combine these two probabilities for each pixel? For example, using mutual information? I think the two probabilities I got are two marginal probabilities, then how to get the joint probability in order to calculate the mutual information? Or you have a better idea to combine the two probabilities please tell me. Thanks very much. A.


title: R rpart cross validation and 1 SE rule, why is the column in cptable called "xstd"?

text: The rpart() function in R returns cptable that includes columns xerror and xstd. Here is an arbitrary example.                           CP nsplit rel error    xerror       xstd     1  0.161992664      0 1.0000000 1.0002790 0.01853630     2  0.043985638      1 0.8380073 0.8385070 0.01749290     3  0.030278222      2 0.7940217 0.7963870 0.01709283     4  0.013881619      3 0.7637435 0.7695997 0.01653832     5  0.010181164      4 0.7498619 0.7560406 0.01606136     6  0.008004043      5 0.7396807 0.7466449 0.01600352     7  0.007026176      6 0.7316767 0.7356289 0.01549501     8  0.006614587      8 0.7176243 0.7388091 0.01559568     9  0.005312278     10 0.7043951 0.7254237 0.01522645     10 0.004883811     11 0.6990828 0.7248227 0.01526605      Some argue that the tree should be pruned based on the minimum cross-validated error (xerror) and thus would prune at row 10, where the minimum xerror occurred. Other argue that "1SE rule" advises to look for the minimum but then go up 1SE because that tree is less complex. Using column xstd, that would suggest using 0.7248227 + 1*0.01526605 = 0.7400887 and thus pruning should occur at row 7. See also this post: How to choose the number of splits in rpart()? **My simple question: why is the column labeled "xstd" (presumably meaning cross-validated standard deviation), and yet people refer to this as the 1SE rule and not the 1SD rule.**


title: Classifying a set of photos to places

text: I want to cluster photos and map them to places. As input I have   * Photos with locations (lat, long)   * Places - some as (imprecise) bounding boxes, some just as points, maybe others as bounding geometries However, all are not exact - photos relating to a building will be taken near, but outside the building. Also, bounding boxes are overlapping (e.g., a photo of Golden Gate Bridge may be outside the bounding box of Golden Gate Bridge and inside California, but should be labeled as Golden Gate Bridge). There are also tags which can serve as additional input - similarity in tags which are not 'generic' but 'local' could serve as another distance metric between photos, I think. But many photos are not tagged well, so this doesn't help much. I'm now using K-means / K-means-plus-plus and getting photo clusters, and now plan to somehow match their bounding box to place bounding boxes and assign them this way (probably by finding the bounding box which has the biggest intersection with photo cluster bounding box). Another option seems to be KNN, but I'm not sure how to best adjust the algorithm for the fact that the places are mostly as bounding boxes, and overlapping. Maybe generate synthetic points for each of the place bounding boxes, then apply KNN? Questions:   * Am I missing some proven algorithms to do this?   * What is the best way to map clusters of points to bounding boxes?   * How to best adapt KNN for the fact that the 'training set' are actually bounding boxes, and that they overlap?


title: Approximation by uniform distribution

text: Imagine I have a random variable $X$ drawn from the sample space $S$ such that $Prob[X = x] = p_x$. Let $A \subseteq S$ be a subset of values from the sample space. My question is if we have a lower bound for the maximum value of $\sum_{x \in A} \frac{p_x}{|A|}$, maximizing over all subsets of $S$. This would mean for me that we could approximate (in some sense) the statistical distribution of $X$ by a uniform distribution over the elements of $A$.


title: Value of r, i.e. correlation coefficient, of data points on a horizontal straight line

text: 1) If we use the formula $$r = \displaystyle \frac{\displaystyle\sum_{i=1}^n \left( x_i - \bar{x} \right) \left( y_i - \bar{y} \right)}{\sqrt{\displaystyle \sum_{i=1}^n \left( x_i - \bar{x}\right)^2 \cdot \sum_{i=1}^n \left( y_i - \bar{y}\right)^2}}$$ for calculating $r$ we get $r = \frac{0}{0}$. 2) If we see $r$ as a measurement of the strength and nature of the linear relationship between two variables, we get $r = 0$ because there is no linear relationship in this case. 3) If a perfect positive fit gives $r = 1$ and a perfect negative fit gives $r = -1$, then it seems like we should get $r = 1$ or $-1$. (I think this is wrong though.) Which is correct? ![enter image description here](http://i.stack.imgur.com/k81lY.png)


title: The size of the sample for split validation

text: At this moment I have a dataset with 4000 samples (50% positive and 50% negative). Normally I would do cross validation for this approach, however besides normal data mining techniques I am also trying an alternate ILP approach. Since I can't implement cross validation using the ILP system I am using, to maintain dataset coherence between the two types of techniques I decided to instead do a split validation. I keep hearing about split validation should only be used in large datasets, but how large should they be? Would implement it on this one produce unacceptable "bad" results? Should I really avoid this approach?


title: Covariance vs Bandwidth of Kernel Density Estimate

text: I've been working with the scipy gaussian kde implementation (here), but I don't quite understand the difference between the bandwidth factor and the covariance matrix. I'm using it for a single dimension here, so the covariance matrix only contains a single value. The bandwidth factor is estimated using either Scott's or Silverman's method, and then apparently it is multiplied with the covariance matrix; but why? What does this covariance matrix actually represent? I understand the bandwidth is the 'width' of the kernel being used, but I don't see the role of the covariance matrix.


title: Constrain regression to correctly rank-order certain cases

text: I am working on a logistic regression model that attempts to predict the probability that individuals will do a certain activity, where each individual can be described by a number of categorical variables. In other words, $$ \log\left(\frac{y}{1-y}\right) = X\beta + \epsilon$$ where $y$ is a binary variable and $X$ is a matrix of categorical variables. In addition, I have a number of "typical" individual profiles that should always be ranked in a specified order. These typical profiles form a constraint, such that the model coefficients should be chosen so as to maximize the likelihood function, but with the constraint that when the model is used to predict the typical profiles, their probabilities are correctly ranked. To take a fanciful example, the independent variables might be hair color and eye color, of which there are many combinations. However, we would like the model to always rank people with (hair=red, eye=green) higher than (hair=yellow, eye=purple), which in turn is ranked higher than (hair=blue, eye=gray), regardless of what the data say. I'm unclear on whether a solution to the problem is actually feasible, that is, whether coefficients can be chosen such that the ranking is guaranteed. If such a guarantee is infeasible, a second-best option would be to have the ranking be "mostly" correct. More important, it's unclear to me how to actually implement this, as I don't know how to characterize it in the typical constrained regression framework, where we are setting linear restrictions on model coefficients (maybe this is not the right way to approach it).


title: In regression, if the residuals are correlated with X then what assumption has most likely been violated?

text: 1. independence   2. homoscedasticity * * * This question is from an online quiz with one "correct" answer, but I'm not sure I agree. I understand that the answer is homoscedasticity, but isn't independence violated here as well? Since heteroskedasticity indicates that residuals are correlated with X, then there is a _dependency_ with x. Therefore, independence is also violated. Could someone please explain?


title: What is an intuitive explanation of why we want homoskedasticity in a regression?

text: I've read that homoskedasticity means that the standard deviation of the error terms are consistent and don't depend on the x-value. Question 1: Can someone explain intuitively why this is necessary? (An applied example would be great!) Question 2: I can never remember whether it's _hetero-_ or _homo-_ that is ideal. Can someone explain the logic of which one is ideal? Question 3: Heteroskedasticity means that x is correlated with the errors. Can someone explain why this is bad? ![Graph of *homo-* vs. *heteroskedasticity*](http://i.stack.imgur.com/RXQv9.png)


title: Estimating a VAR model with variable coefficients

text: I want to estimate a VAR model based on the Dufour and Engle paper "Time and the Price Impact of a Trade" (2000). There, the parameter $ b_{i} $ of the endogenous variable $ x_{i} $ is dependent on the strictly exogenous variable $ T_{t-i} $ : $$ b_{i}x_{t-i}=[ \gamma_{i} +\delta_{i} T_{t-i} ] x_{t-i} \ \ (1) $$ The authors state, that they estimated this system by OLS, but I am not quite sure how the estimator looks like in this case. I thought of multiplying out the equation and defining a new variable, say $ z_{t-i} $ $$ \gamma_{i} x_{t-i} +\delta_{i} [ T_{t-i} x_{t-i}] = \gamma_{i} x_{t-i} +\delta_{i} z_{t-i} $$ which would consist of an exogenous and endogenous part though, which I think makes this procedure methodically incorrect. Can anybody give an explanation how such a system as in (1) is estimated using OLS? I plan to use the vars- Package in R btw.


title: Generlized Cross Validation in Matlab

text: I want to use Generalized Cross Validation to compute the optimal ridge parameter for a ridge regression. The formula due to Golub et al (1979) is GCV(k)=(1/n Sum_i=1^n e_i(k)^2)/(1/n Sum_i=1^n(1-h_ii(k)))^2 where e_1(k)....e_n(k) are the residuals from applying the ridge regression for a given k. H(K)=X(X'X+kI)^{-1}X'. where X is the design matrix. and h_ii(k) are the digonal elelments of H(k). I have the ridge regression formula so I can easily compute the numerator. However I am not sure how to compute H? I found this set of tools which features a function gcv() but I can't understand how it it calculating the G function? http://www2.imm.dtu.dk/~pcha/Regutools/ Indeed it does not appear to be using x in it's calculations? Is this function going what I want above? If so can someone please explain to me how it works?


title: Handling Infrequently Occurring Categorical Variables

text: I'm dealing with some data where there are some infrequently occuring categorical variables related to a binary prediction target. For example marketing partners... some send 1000s of leads but many others send less than 10 a minuscule number in a data set with 20K+ examples. Typically, we deal with this by grouping all the infrequently occurring partners under the variable "Else" and using this as the reference category in logistic regression. I'm wondering what the pros and cons are of handling the issue like this and if there is a better way to deal with it in pre-processing phase.


title: Proof of Convergence in Probability

text: **Theorem** Suppose $X_n \stackrel{P}\to X\quad \text{and}\quad Y_n \stackrel{P}\to Y $. Then $X_n+Y_n \stackrel{P}\to X+Y $. **Proof** Let $\epsilon >0 $ be given. Using the triangle inequality, we can write $$ | X_n-X|+|Y_n-Y| \geq |\left( X_n+Y_n \right) -\left( X+Y \right) | \geq \epsilon $$ Since **$P$ is monotove relative to set containment** , we have $$ \begin{align} &P \left[| \left(X_n+Y_n \right)- \left(X+Y \right)| \geq \epsilon \right] \\\ & \\\ &\leq P \left[|X_n-X|+|Y_n-Y| \geq \epsilon \right] \\\ & \\\ &\leq P \left[|X_n-X| \geq \epsilon /2 \right]+P \left[|Y_n-Y| \geq \epsilon/2 \right] \end{align} $$ Then by the hypothesis of the theorem the last two terms go to $0$ which gives us the desired result. First of all sorry for the messy LaTeX notation, I did not know how to align those two inequalities in order for the second to appear right below the first, although I tried. My questions are what exactly does **"P is monotone relative to set containment"** mean and also how is the last inequality derived? Thank you.


title: Why do people like smooth data?

text: I am to use the Squared Exponential kernel (SE) for Gaussian Process Regression. The advantages of this kernel are: 1) simple: only 3 hyperparameters; 2) smooth: this kernel is Gaussian. Why do people like 'smoothness' so much? I know that the Gaussian kernel is infinitely differentiable, but is that so important? (Please let me know if there is other reasons why the SE kernel is so popular.) PS: I was told that most signals in real world (without noise) are **smooth** , so it is reasonable to use smooth kernels to model them. Could anyone please help me understand this concept?


title: How are the standard errors computed for the fitted values from a logistic regression?

text: When you predict a fitted value from a logistic regression model, how are standard errors computed? I mean for the _fitted values_ , not for the coefficients (which involves Fishers information matrix). I only found out how to get the numbers with `R` (e.g., here on r-help, or here on Stack Overflow), but I cannot find the formula.               pred <- predict(y.glm, newdata= something, se.fit=TRUE)      If you could provide online source (preferably on a university website), that would be fantastic.


title: Hettest Not Appropriate for Robust Cluster

text: I am using stata, and running a regression with robust standard errors. I want to test the Heteroscedasticity after running the regression with robust standard errors, and to compare this to before the regression was run using this facility? Is this possible in anyway?


title: Comparing levels of factors after a GLM in R

text: Here is a little background about my situation: my data refer to the number of prey successfully eaten by a predator. As the number of prey is limited (25 available) in each trial, I had a column "Sample" representing the number of available prey (so, 25 in each trial), and another called "Count" which was the number of success (how many prey were eaten). I based my analysis on the example from the R book on proportion data (page 578). The explanatory variables are Temperature (4 levels, which I treated as factor), and Sex of the predator (obviously, male or female). So I end up with this model:               model <- glm(y ~ Temperature+Sex+Temperature*Sex data=predator, family=quasibinomial)       After getting the Analysis of Deviance table, it turns out that Temperature and Sex (but not the interaction) have a significant effect on the consumption of prey. Now, my problem: I need to know which temperatures differ, i.e., I have to compare the 4 temperatures to each other. If I had a linear model, I would use the TukeyHSD function, but as I am using a GLM I can't. I have been looking through the package MASS and trying to set up a contrast matrix but for some reason it doesn't work. Any suggestions or references? Here's the summary I get from my model, if that helps making it clearer...               y <- cbind(data$Count, data$Sample-data$Count)     model <- glm(y ~ Temperature+Sex+Temperature*Sex data=predator, family=quasibinomial)      > summary(model)          Call:     glm(formula = y ~ Temperature + Sex + Temperature * Sex, family=quasibinomial, data=data)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.7926  -1.4308  -0.3098   0.9438   3.6831            Coefficients:                                            Estimate Std. Error t value Pr(>|t|)         (Intercept)                             -1.6094     0.2672  -6.024 3.86e-08 ***     Temperature8                             0.3438     0.3594   0.957   0.3414         Temperature11                           -1.0296     0.4803  -2.144   0.0348 *       Temperature15                           -1.2669     0.5174  -2.449   0.0163 *       SexMale                                  0.3822     0.3577   1.069   0.2882         Temperature8:SexMale                    -0.2152     0.4884  -0.441   0.6606         Temperature11:SexMale                    0.4136     0.6093   0.679   0.4990         Temperature15:SexMale                    0.4370     0.6503   0.672   0.5033         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for quasibinomial family taken to be 2.97372)             Null deviance: 384.54  on 95  degrees of freedom     Residual deviance: 289.45  on 88  degrees of freedom     AIC: NA        Number of Fisher Scoring iterations: 5


title: Is the "test statistic" a value or a random variable?

text: I am a student taking my first Statistics course now. I am confused by the term "test statistic". In the following (I saw this in some textbooks), $t$ seems to be a specific value calculated from a specific sample. $$ t=\frac{\overline{x} - \mu_0}{s / \sqrt{n}} $$ However, in the following (I saw this in some other textbooks), $T$ seems to be a random variable. $$ T=\frac{\overline{X} - \mu_0}{S / \sqrt{n}} $$ So, does the term "test statistic" mean a specific value or a random variable, or **both**? Thank you!


title: How does variance change as sample size increases

text: ## Situation: n possibilities each have their own probability of happening, and their own payout when they do. So expected payout $E_n$ is $\sum\limits_{i=1}^n \text{probability}_i*\text{payout}_i$ And the expected variance is $\sum\limits_{i=1}^n ( \text{payout}_i- \text{expected payout})^2* \text{probability}_i$ ## Question: If this activity is repeated 5 times, the expected payout is $5* \text{expected payout}$ But the variance is not 5 times the expected variance; it should get proportionally smaller as this activity goes on. But how do I calculate the variance through iterations?


title: Why do we use fixed effects and controls at the same time?

text: When we run a regression model (say OLS for simplicity) of y ~ x, we might have to use several control variables say z1 and z2. Now our model is y ~ x + z1 + z2, we may believe that z1 and z2 are not perfect controls then we introduce fixed effect dummies. Now our model is y ~ x + z1 +z2 + fixed effects; but if the fixed effects take into account all control information then surely they must be seriously correlated with z1 and z2. Otherwise, by definition z1 and z2 are junk. Is this not going to introduce multi-collinearity?


title: What's the appropriate way to calculate an effect size when there is an interaction?

text: I have a regression model with two significant main effects and a significant interaction. I'd like to calculate an effect size for one of the main effects. I'm using the package `compute.es` in R. I'm not sure whether I should use the empirical means and standard deviations or the t-value from the regression model to calculate the effect size.


title: Dealing with perfectly confounded microarray experiment

text: I need to compare microarray data, where all of the "cases" were hybridized in one batch and all of the "controls" in another, so I have no way of removing this batch effect. What would be the best strategy to detect differentially expressed genes (if at all possible)?


title: Estimating the slope of the straight portion of a sigmoid curve

text: I have been given this task and was stumped. A colleague asked me to estimate the $x_{upper}$ and $x_{lower}$ of the following chart: ![enter image description here](http://i.stack.imgur.com/FLDQU.png) The curve is actually a cumulative distribution, and x is some kind of measurements. He is interested to know what are the corresponding values on x when the cumulative function started to become straight and deviate from being straight. I understand that we can use differentiation to find the slope at a point, but I am not too sure how to determine when can we call the line straight. Any nudge towards some already existing approach/literature will be greatly appreciated. I know R as well if you happen to know any relevant packages or examples on this kind of investigations. Thanks a lot. * * * _UPDATE_ Thanks to Flounderer I was able to expand the work further, set up a framework, and tinker the parameters here and there. For learning purpose here are my current code and a graphic output.               library(ESPRESSO)          x <- skew.rnorm(800, 150, 5, 3)     x <- sort(x)     meanX <- mean(x)     sdX <- sd(x)     stdX <- (x-meanX)/sdX     y <- pnorm(stdX)          par(mfrow=c(2,2), mai=c(1,1,0.3,0.3))     hist(x, col="#03718750", border="white", main="")          nq <- diff(y)/diff(x)     plot.ts(nq, col="#6dc03480")          log.nq <- log(nq)     low <- lowess(log.nq)     cutoff <- .7     q <- quantile(low$y, cutoff)     plot.ts(log.nq, col="#6dc03480")     abline(h=q, col="#348d9e")          x.lower <- x[min(which(low$y > q))]         x.upper <- x[max(which(low$y > q))]     plot(x,y,pch=16,col="#03718750", axes=F)     axis(side=1)     axis(side=2)     abline(v=c(x.lower, x.upper),col="red")     text(x.lower, 1.0, round(x.lower,0))     text(x.upper, 1.0, round(x.upper,0))      ![enter image description here](http://i.stack.imgur.com/F38Uk.png)


title: Simplify calculation of mortgage repayments in R

text: I've just started learning R, but I feel that I'm using it only as a normal procedural language without taking advantage of some built-in functions and data structures. As an example I wrote a function that calculates repayments for mortgage with an option of specifying list of overpayments to reduce monthly payments. I would appreciate if you could suggest any improvements how to make it more native to R.               source("http://pastebin.com/raw.php?i=q7tyiEmM") ## for pmt()     mortgage <- function(price=1000000, deposit=0.10, rate=0.04, years=25, over=c()) {         loan <- price - price * deposit         months <- years * 12         capital <- c()         interest <- c()         payments <- c()         currentCapital <- loan         monthlyRate <- rate / 12         payment <- -pmt(monthlyRate, nper = months + 1, pv = loan)              cat("Initial capital: ", loan, "\n")         cat("Months: ", months, "\n")         cat("Payment: ", payment, "\n")              totalShouldPay <- payment * months         totalPaid <- 0              for(currentMonth in 1:months) {             interest <- currentCapital * monthlyRate             repaid <- payment - interest             currentCapital <- currentCapital - repaid                  overPayments <- subset(over, over$month == currentMonth)           if (length(overPayments) > 0 && nrow(overPayments) > 0) {            overPayments <- sum(overPayments$payment)                 currentCapital <- currentCapital - overPayments                 payment <- -pmt(monthlyRate, nper = months + 1 - currentMonth, pv = currentCapital)             }                  capital <- append(capital, currentCapital)                  totalPaid <- totalPaid + payment             payments <- append(payments, payment)         }         cat("Would have paid: ", round(totalShouldPay, 2), "\n")         cat("Will pay: ", round(totalPaid, 2), "\n")         cat("Saved: ", round(totalShouldPay - totalPaid, 2), "\n")         return(data.frame(capital=capital, payments=payments))     }          overPayments <- data.frame(         month   = c(20,    25,    30,    35,    40),          payment = c(20000, 20000, 20000, 20000, 20000))          cap <- mortgage(         over = overPayments     )          #print(cap)     par(mfcol=c(1, 2))     plot(cap$payments, type="s", main = "Payments", xlab = "Month", ylab = "Payment")         plot(cap$capital, type="l", main = "Capital", xlab = "Month", ylab = "Remaining loan")


title: How do I produce a set of predictions based on a new set of data using predict in R?

text: I'm struggling to understand how the predict function works and can be used with out of sample data. For instance the following code...               my <- data.frame(x=rnorm(1000))       my$y <- 0.5*my$x+0.5*rnorm(1000)     fit <- lm(my$y ~ my$x)        mySample <- my[sample(my, 100),]         predict(fit, mySample)      I would understand should return 100 y predictions based on the sample. But it returns 1,000 row with the warning message :               'newdata' had 100 rows but variables found have 1000 rows      How do I produce a set of predictions based on a new set of data using predict? Or am I using the wrong function? I am a noob so apologise in advance if I am asking stupid questions.


title: Relationship between the latent variable as a function of regressors and the logit model?

text: Can anyone give the intuition behind of the relationship between these two? I see a lot of proofs in books, but no real intuition. Thanks


title: Different group size average comparisons

text: I have 3800 stocks for which I calculate their return correlations. So I end up with a 3800 by 3800 matrix. Each stocks is related in some way to all the other 3799 stocks. For each stock, I create seven groups out of 3799 stocks based on their relatedness. So say for stock A, `group 1` (closest in relatedness - based on a pre-specified knowledge based linkage of like- industries) will have 60 stocks, `group 2` will have 300 stocks, and `group 3` 700 stocks and so on up to seven groups. I then average for each group their correlations with stock A. I want to see if there is a decay in the correlation when we move from `group 1` to `group 7`. But the group size varies... can this be a serious issue? I then take an average across all `group`s for each stock and see if the decay holds cross-sectionally? Now the issue is that for each stock, say `stock A` and `stock B` the size of each group is different. There may be 60 stocks in `group 1` for `stock A` while there is 40 in `group 1` for `stock B`. How problematic can this be?


title: Obtaining an estimate of slope error on a linear fit to data with y error

text: I have a problem similar to what was discussed here: I have a dataset (x_i,y_i) with y-error which is nonuniform over the interval of interest. I want to fit a line to this and extract a value of uncertainty of the slope of the fit which takes into account the y errors. I think it's safe to assume the errors are normal with mean zero. I am using Origin for my fits, and as far as I can tell it is only calculating the basic least squares standard error and ignoring my error bars. Unfortunately I am not familiar with Mathematica so I am unable to decipher the codes given in the linked question. They seem to do what I want but I can't make sense of them. Can someone suggest a different route?


title: Gamma vs. lognormal distributions

text: I have an experimentally observed distribution that looks very similar to a gamma or lognormal distribution. I've read that the lognormal distribution is the maximum entropy probability distribution for a random variate $X$ for which the mean and variance of $\ln(X)$ are fixed. Does the gamma distribution have any similar properties?


title: How to make use of known constants when modeling from data?

text: As a (perhaps contrived) example, let's say we want to discover from some empirical data Coulomb's law for an electric field: $$F = \frac{1}{4 \pi \epsilon_0} \cdot \frac{|q|}{r^2}$$ In this case we would have a matrix whose columns include force $F$, charge $q$, and radius $r$; as well as perhaps some columns containing irrelevant features. We also know in advance the values of $\epsilon_0$ and $\pi$. Supplying that matrix as input to some out-of-box machine learning tools may indeed have us find that $$F = c \cdot \frac{|q|}{r^2}$$ for some $c$ that is close to but probably not exactly $(4 \pi \epsilon_0)^{-1}$. My question is: how do we inform such systems to make use of the known constants so that it may more closely approximate the true form of the equation? Would it useful to include, for example, an entire column containing only $\pi$ and another entire column containing only $\epsilon_0$? Edit: It may be the case here that $\epsilon_0$ is the interesting quantity to be discovered (in addition to the form of the equation). However, the question still stands regarding how to make use of $\pi$.


title: simultaneous equations

text: I have the following relationships logY ~ logX1 + logX2 + logX3 + logX4 + logX5 and X1 ~ Z1 + Z2 + Z3 + Z4 + Z5 X2 ~ Z1 + Z2 + Z3 + Z4 + Z5 X3 ~ Z1 + Z2 + Z3 + Z4 + Z5 where Y and Z1, Z2, Z3, Z4, Z5 are endogenous (Say while the Z's play a role in determining Y, the values of Z's are fixed depending upont he values of Y - Kind of like advertising expense has an impact on sales revenue but at the same time managers determine the advertisement expense on the expected sales revenue). So all the variable are changing simultaneously. Can anyone help me on how I can estimate this relationship? I also have instruments for each of the Z's (lagged variables have been treated as instruments and I have the previous year data for the problem as well. Thank you for all your help and suggestions.


title: Find music tempo with a audio signal

text: I'd like to know if there is statistical techniques that aims to find the rhythm/tempo of a music based on audio signal : ![Audio signal](http://i.stack.imgur.com/KgbnR.png)


title: How to treat age-specific left-censored data

text: I would like to study how the death of a respondent's parent affects his personality (or whatever DV). I have an age of respondent, and also an age when respondent lost his parent. The problem is that these records are available only 23 years back. That is if somebody lost his parents 25 years ago this information is not available. That means that the two variables (age and age-at-parent's-death. APD) are highly correlated - because nobody can lose his parents at age that exceeds his own, and nobody can lose them earlier than 23 years back (because there is no this info in the dataset). If I want to study if the loss of a parent at early age affects him more than when he gets older, is it any reliable way to do it with existing dataset? Here is scatterplot of two variables (age, and APD): ![enter image description here](http://i.stack.imgur.com/rMETW.jpg)


title: Average of performance estimates

text: When I learn several classifiers for a task and want to compare their performance on it, I thought about calculating several performance measures, such as:   * sensitivity   * specificity   * positive and negative predictive value   * ... Then of course I can compare the single performance values, but does it statistically make sense to compare the average of all performance values also? If so, then I guess just under the assumption that all performance measures are worth the same? Thanks!


title: kurtosis, positive skewed and negative skewed for probability distribution

text: When discussing probability distribution, I always read something such as excess kurtosis, positive kurtosis, positive skewed and negative skewed. What exactly do these concepts indicate? In practical applications, such as market return or something else, what can these characteristics tell us?


title: What is the name for this distribution defined on a circle

text: One can define the probability distribution: $$ p(\theta; \alpha, \theta_0) = \frac{ e^{ \alpha \cos( \theta-\theta_0) }}{ 2 \pi I_0(\alpha)} $$ over an angular variable $\theta \in [0,2 \pi]$. By what name(s) is this distribution called? ($I_0$ is the modified Bessel function and serves to normalize the distribution)


title: What are the connections between: the normal, the $\chi^2$ & the F distributions?

text: I have often read that there is a huge connection between the normal distribution and several other distributions. But these were only mathematical explanations. What's the "real" connection between these 3 distributions (normal, $\chi^2$, and F)?


title: Are there some nice introductions to IDI, NRI, biomarker validation?

text: Are there some nice introductions to IDI (integrated discrimination improvement index), NRI (net reclassification index), multivariable likelihood ratios, Bayesian methods, biomarker validation, risk scores from multi- biomarker data? Where do they arise and how are they obtained or performed? I guess the topics are from survival analysis. But I can't find it in some survival analysis books.


title: Why bother with low rank approximations?

text: If you have a matrix with n rows and m columns, you can use SVD or other methods to calculate a low-rank approximation of the given matrix. However, the low rank approximation will still have n rows and m columns. How can low-rank-approximations be useful for machine learning and natural language processing, given that you are left with the same number of features?


title: Optimal matching distances as percentages?

text: I'm using R and the library TraMineR to do sequence analysis on social science data. I have computed an optimal matching (OM) distance matrix.                     [,1]     [,2]      [,3]      [,4]     [1,]   0.00000 94.80628 110.00000 107.46826     [2,]  94.80628  0.00000  53.71086  78.02520     [3,] 110.00000 53.71086   0.00000  47.88362     [4,] 107.46826 78.02520  47.88362   0.00000      Here is a dput() of the OM matrix:               structure(c(0, 94.8062777425491, 110, 107.468255963965, 94.8062777425491,      0, 53.7108606841158, 78.0251951544836, 110, 53.7108606841158,      0, 47.8836216007405, 107.468255963965, 78.0251951544836, 47.8836216007405,      0), .Dim = c(4L, 4L))      It is based on this data: https://github.com/aronlindberg/VOSS-Sequencing- Toolkit/blob/master/twitter_exploratory_analysis/cassandra.csv and I computed it using this code:               # Read and transpose data     twitter_sequences <- read.csv(file = "cassandra.csv", header = TRUE)     twitter_sequences_transposed <- t(twitter_sequences)          # Define the sequence object     twitter.seq <- seqdef(twitter_sequences_transposed, left="DEL", right="DEL", gaps="DEL", missing="")          # Define the substitution matrix     twitter_costs <- seqsubm(twitter.seq, method="TRATE")          # Calculate the OM matrix     twitter.om <- seqdist(twitter.seq, method="OM", indel=1, sm=twitter_costs, with.missing=FALSE)      Now I want to express the similarity between two sequences in %. How can I accomplish this?


title: Comparing the proportion of a subsample to the proportion of sample or another

text: I have one sample and three sub-samples and I would like to calculate the p-value for the hypothesis that the proportion in the sub-sample is not different from the proportion in another sub-sample or not different from the proportion in the sample. The sample consists of counting how often there is agreement between one individual and another (1 if they agree and 0 otherwise). The first sub-sample removes some observations from the sample. The second sub-sample removes some observations from the first sub-sample. The third sub-sample removes other observations from the first sub-sample. The sample sizes of the sample and sub-samples are 400, 340, 260, 170. I have computed the 95 % confidence interval and find that these largely overlap for the first three, but the fourth confidence interval overlaps little if at all with the other three. I have also computed the p-values for the hypothesis that the proportions are the same under the assumption that the samples are independent. The results here are of course the same. What I would like to know is how large the effect is when the assumption of independence is violated, does anyone know by how much it can skew the p-values? Otherwise I'd like to know if there is a better way to compute the p-values for this problem? Should I use the McNemar test (I encountered it in some of the answers around the web, though I am not sure if it applies)? Preferably I'd also like to know how to compute the power of such a test. Cheers, Martin


title: Ordinal Regression

text: In my data, the response variable is ordinal and there are 16 predicting variables which are also ordinal. I plan to perform an ordinal regression. But before doing that I want to graphically show the relationship between the response and predictor variables. In regular regression where the variables are numerical, we can produce a scatterplot matrix for that purpose. I am wondering whether for ordinal regression with ordinal predictors there is any graph that can be used to visually display the relationship before the formal analysis. Another possibility is to calculate the association between the response and the ordinal predictors before the analysis. Any recommendation on that part? Also, I plan to perform a proportional odds model using R. Can anyone recommend a R package or functions which can be used to test the proportionality and also used for the analysis?


title: Specific questions: Statistical tests, fair analysis

text: **Background:** My dataset is small (3 to 10 values). But each set comes from 18 acquisitions of different areas of the same sample. So they would not be considered independent. But 3-10 values are too small to analyze whether the distribution is normal. And the raw values differ a lot between independent experiments but we express them as percentage of the 'control' and then are very comparable. **Questions:**   1. We measure the background noise in our experiment and subtract it from each dataset. Does this skew the standard deviation of the samples? Should it not be subtracted?   2. Is it fair to test our sample against the control which will always be 100%?   3. Also, although I can calculate the standard deviation of the control within each independent experiment, how can I express the standard deviation between independent experiments? I'm sure you can't just average the standard deviations!   4. For comparing one sample to the control, is the Wilcoxon signed rank test reasonable for such a non-parametric low sample data set? Can we compare it to the hypothetical value of 100%?   5. For comparing a sample treated with increasing doses of a chemical (multiple datasets), is the Kruskall Walis test a reasonable one? Can it also be compared to the hypothetical value of 100%?


title: How do you calculate the effect size of one-sample Wilcoxon signed-rank test?

text: Can you / how do you calculate the effect of a one-sample Wilcoxon signed-rank test? Is there a way to perform the calculation in SPSS or another software package?


title: How to analyze full-factorial within-subject experiment with ordinal (likert) data?

text: How can I analyze an within-subject full-factorial experiment with ordinal data from likert scale (with repetitions)? I would greatly appreciate your help. **Study design** Here is a description of the study design, below you can find the data we collected for clarity. We conducted a _within-subject_ experiment where we tested 2 independent variables (A & B) with 3 levels each (0, 1 & 2). Here is a schematic of our conditions:                    | A(0) | A(1) | A(2) |     B(0) |   1  |   2  |   3  |     B(1) |   4  |   5  |   6  |     B(2) |   7  |   8  |   9  |      For each condition, participants rated the conditions on a 7-point likert scale (1=bad, 7=good) which is our dependent variable. We also collected multiple data points per condition and participant: we repeated a block of all 9 conditions (order of conditions was randomized) 3 times. Thus, we collected 3 data points per condition and participant (27 data points in total per participant). Here is our collected data for your reference. ![all results](http://i.stack.imgur.com/2q6kI.png) **Question** Now I am wondering how I should analyse the data. This would be a classical case for repeated-measures ANOVA, if the data would not be ordinal. So, do I have to use a non-parametric test due to the ordinal nature of my data? (A Friedman-test in this case, because it's a within-subject design?) If so, how would I correct for multiple comparisons? Or could I use a repeated-measures ANOVA even on ordinal data, if the data would be normally distributed? I've also read that I could (align) rank transform my data and then run an (two-way) ANOVA, but can I also do this for within-participant designs? I've only seen examples where this was used for between-subject designs. My stats knowledge is beginner level, so low-level explanations would be appreciated. I'm using SPSS (but could also use other packages if they help solve this :)


title: Paired t-test as a special case of linear mixed-effect modeling (still unresolved)

text: We know that a paired _t_ -test is just a special case of one-way repeated- measures (or within-subject) ANOVA as well as linear mixed-effect model, which can be demonstrated with lme() function the nlme package in R as shown below.               #response data from 10 subjects under two conditions     x1<-rnorm(10)     x2<-1+rnorm(10)          # Now create a dataframe for lme     myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), rep(paste("S", seq(1,10), sep=""), 2))     names(myDat) <- c("y", "x", "subj")      When I run the following paired t-test:               t.test(x1, x2, paired = TRUE)      I got this result (you will get a different result because of the random generator):               t = -2.3056, df = 9, p-value = 0.04657      With the ANOVA approach we can get the same result:               summary(aov(y ~ x + Error(subj/x), myDat))          # the F-value below is just the square of the t-value from paired t-test:               Df  F value Pr(>F)     x          1  5.3158  0.04657      Now I can obtain the same result in lme with the following model, assuming a positive-definite symmetrical correlation matrix for the two conditions:               summary(fm1 <- lme(y ~ x, random=list(subj=pdSymm(form=~x-1)), data=myDat))          # the 2nd row in the following agrees with the paired t-test     # (Intercept) -0.2488202 0.3142115  9 -0.7918878  0.4488     # xx2          1.3325786 0.5779727  9  2.3056084  0.0466      Or another model, assuming a compound symmetry for the correlation matrix of the two conditions:               summary(fm2 <- lme(y ~ x, random=list(subj=pdCompSymm(form=~x-1)), data=myDat))          # the 2nd row in the following agrees with the paired t-test     # (Intercept) -0.2488202 0.4023431  9 -0.618428  0.5516     # xx2          1.3325786 0.5779727  9  2.305608  0.0466      With the paired t-test and one-way repeated-measures ANOVA, I can write down the traditional cell mean model as               Yij = μ + αi + βj + εij, i = 1, 2; j = 1, ..., 10      where i indexes condition, j indexes subject, Yij is the response variable, μ is constant for the fixed effect for overall mean, αi is the fixed effect for condition, βj is the random effect for subject following N(0, σp2) (σp2 is population variance), and εij is residual following N(0, σ2) (σ2 is within- subject variance). I thought that the cell mean model above would not be appropriate for the lme models, but the trouble is that I can't come up with a reasonable model for the two lme() approaches with the correlation structure assumption. The reason is that the lme model seems to have more parameters for the random components than the cell mean model above offers. At least the lme model provides exactly the same F-value, degrees of freedom, and p-value as well, which gls cannot. More specifically gls gives incorrect DFs due to the fact that it does not account for the fact that each subject has two observations, leading to much inflated DFs. The lme model most likely is overparameterized in specifying the random effects, but I don't know what the model is and what the parameters are. So the issue is still unresolved for me.


title: Crossed random effects and unbalanced data

text: I am modeling some data where I think I have two crossed random effects. But the data set is not balanced, and I'm not sure what needs to be done to account for it. My data is a set of events. An event occurs when a client meets with a provider to perform a task, which is either successful or not. There are thousands of clients and providers, and each client & provider participates in varying numbers of events (roughly 5 to 500). Each client and provider has a level of skill, and the chance that the task is successful is a function of the skills of both participants. There is no overlap between clients and providers. I am interested in the respective variances of the population of clients and providers, so we can know which source has a bigger effect on the success rate. I also want to know the specific values of the skills among the client and providers we actually have data for, to identify best/worst clients or providers. Initially, I want to assume that the probability of success is driven solely by the combined skill levels of the client and provider, with no other fixed effects. So, assuming that x is a factor for the client and y is a factor for provider, then in R (using package lme4) I have a model specified as:                 glmer( success ~ (1 | x) + (1 | y), family=binomial(), data=events)      One problem is that clients are not evenly distributed across providers. Higher skill clients are more likely to be matched up with higher skill providers. My understanding is that a random effect has to be uncorrelated with any other predictors in the model, but I'm not sure how to account for it. Also, some clients and providers have very few events (less than 10), while others have many (up to 500), so there's a wide spread in the amount of data we have on each participant. Ideally this would be reflected in a "confidence interval" around each particpants skill estimate (although I think the term confidence interval isn't quite correct here). Are crossed random effects going to be problematic because of the unbalanced data? If so, what are some other approaches I should consider?


title: Homogeneity of variance for time variable

text: I have four groups of plant treated with different temperatures and I conducted repeated measurement on their growth rate. I use ANOVA with Mixed Model to analyze the data by specifying both temperature (TEMP), time (TIME) and their interaction (Temp*Time) as fixed factor, I also include TIME as repeated factor using plant as subject. I used Levene's test to test the assumption of homogeneity of variance. For the temperature factor, the result is not significant. However, for the Time factor, the result is significant (p=0.012). May I know do I need to have homogeneity of variance for Time as well? Further more, do I need to do Levene's test on the interaction (Temp*Time) as well? If I could not have homogeneity of variance for Time even after transformation, should I do One-Way ANOVA separately for each time point? Or should I use a non-parametric test? Or is there any other methods to analyze the data? Thank you so much for all the generous help!


I want to calculate: $\operatorname{var}[\frac{b}{a} B(a-b)-b B(b)]$ with $b\leq a$ and $b\geq 0$; $B=$brownian motion. I started like this: $(\frac{b}{a})^2 \operatorname{var}[ B(a-b)]+-b ^2 \operatorname{var}[ B(b)]+ 2 \operatorname{Cov}(B(a-b), B(b))$ Can someone help me here or comment on my approach?


title: Bayesian and Fisher's approaches to linear discriminant analysis

text: I know 2 approaches to do LDA, the _Bayesian approach_ and the _Fisher's approach_. Suppose we have the data $(x,y)$, where $x$ is the $p$-dimensional predictor and $y$ is the dependent variable of $K$ classes. By _Bayesian approach_ , we compute the posterior $$p(y_k|x)=\frac{p(x|y_k)p(y_k)}{p(x)}\propto p(x|y_k)p(y_k)$$, and as said in the books, assume $p(x|y_k)$ is Gaussian, we now have the discriminant function for the $k$th class as \begin{align*}f_k(x)&=\ln p(x|y_k)+\ln p(y_k)\\\&=\ln\left[\frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}}\exp\left(-\frac{1}{2}(x-\mu_k)^T\Sigma^{-1}(x-\mu_k)\right)\right]+\ln p(y_k)\\\&=x^T\Sigma^{-1}\mu_k-\frac{1}{2}\mu_k^T\Sigma^{-1}\mu_k+\ln p(y_k)\end{align*}, I can see $f_k(x)$ is a linear function of $x$, so for all the $K$ classes we have $K$ linear discriminant functions. However, by _Fisher's approach_ , we try to project $x$ to $(K-1)$ dimensional space to extract the new features which minimizes **within-class variance** and maximizes **between-class** variance, let's say the projection matrix is $W$ with each column being a projection direction. This approach is more like a **dimension reduction** technique. My questions are (1) Can we do dimension reduction using Bayesian approach? I mean, we can use the Bayesian approach to do classification by finding the discriminant functions $f_k(x)$ which gives the largest value for new $x^*$, but can these discriminant functions $f_k(x)$ be used to project $x$ to lower dimensional subspace? Just like _Fisher's approach_ does. (2) Do and how the two approaches relate to each other? I don't see any relation between them, because one seems just to be able to do classification with $f_k(x)$ value, and the other is primarily aimed at dimension reduction. **UPDATE** Thanks to @amoeba, according to ESL book, I found this: ![enter image description here](http://i.stack.imgur.com/CCMxo.png) and this is the linear discriminant function, derived via Bayes theorem plus assuming all classes having the same covariance matrix $\Sigma$. And this discriminant function is the SAME as the one $f_k(x)$ I wrote above. Can I use $\Sigma^{-1}\mu_k$ as the direction on which to project $x$, in order to do dimension reduction? I'm not sure about this, since AFAIK, the dimension reduction is achieved by do the **between-within variance** analysis. **UPDATE AGAIN** From section 4.3.3, this is how those projections derived: ![enter image description here](http://i.stack.imgur.com/vxjWl.png) , and of course it assumes a shared covariance among classes, that is the _common covariance matrix $W$ (for within-class covariance)_ , right? My problem is how do I compute this $W$ from the data? Since I would have $K$ different within-class covariance matrices if I try to compute $W$ from the data. So do I have to **pool** all class' covariance together to obtain a common one?


title: When measuring attitudes, does rank-ordering choices have any advantages over rating choices?

text: I can see an advantage for rating in that you potentially end up with interval data rather than ordinal data. However, I'm struggling to think of what potential advantages there might be to rank-ordering.


title: How to obtain the standard deviation from means or confidence interval of the mean difference?

text: I am trying to calculate the sample size for a new study based on a previous one. The previous study compared means using a _t_ test. **I have:** the means, the _t_ test result, the confidence interval of the mean difference, and the effect size. However, the calculator that I am using to get the sample size requires the standard deviation.   **How can I obtain the standard deviation from what I have?**


title: How is the poisson distribution a distribution? It seems more like a formula

text: I just watched this video: https://www.youtube.com/watch?v=Fk02TW6reiA It shows a formula to calculate an answer for the following problem:   * There are 2 customers expected every 3 minutes in a store   * Therefore there are 6 customers expected every 9 minutes   * What is the likelihood of there being 4 or less in the store in 9 minutes?   * Answer is: P(0;6)+...+P(4;6) which is about 0.28 This makes sense and is well described. However, numpy treats the poisson distribution basically like a random number generator: http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.poisson.html We can specify lambda as say 5 and how many numbers are desired (the second argument) and get a big list of integers:               >>> import numpy as np     >>> s = np.random.poisson(5, 10000)     >>> s     array([2, 4, 4, ..., 3, 4, 3])     >>> len(s)     10000      These seem like two totally different things. How do you get from using the Poisson formula to calculate the possibility of a certain number of events in a timeframe, to a list of seemingly random integers?


title: conducting multi-level regression on ordinal DVs with imputed data in R

text: Do you know of an approach/package that facilitates mixed model regression of ordinal dependent variables on multiply imputed datasets in R? Ideally, the function takes: a list of multiply imputed datasets a list of target variables (dependent variables, one or more) a list of factors (independent variables, one or more) a list of dummy coded conditions (for analysis of multifactor IVs) and returns a table similar to the result of other regressions in r After extensive searching, I had to create such a function using CLMM in the ordinal package. If you can't answer the first question perhaps you can advise me wrt code adaptations, statistical appropriateness of my approach, efficiency (it takes a LONG time with many imputed datasets), etc here's some data that mirrors mine               numimpdatasets = 3; N = 170; datalist = list()     for (datasetnum in 1:numimpdatasets){         dvone = sample(1:5, N, replace=T)            dvtwo = sample(1:5, N, replace=T)         teacher = c('Tom','Dick','Harry')[sample(1:3, N, replace=T)]         studentclass = c('Class1','Class2','Class3','Class4','Class5',                          'Class6')[sample(1:6, N, replace=T)]         aptitude = runif(N, -3.5, 3.5)         randord = sample(1:3, N, replace=T)         conddummycode1 = c('a_cond1', 'b_cond2', 'c_cond3')[randord]         conddummycode2 = c('c_cond1', 'b_cond2', 'a_cond3')[randord]         datalist[[datasetnum]] = data.frame(cbind(dvone, dvtwo,teacher,                 studentclass, aptitude, conddummycode1, conddummycode2))         }     dvs = colnames(datalist[[1]])[1:2]     conditions = colnames(datalist[[1]])[6:7]     ivs = c("+as.numeric(aptitude)","+(1|teacher/studentclass)")


title: Hierarchical Linear Modeling - Summation in Level 1

text: I'm currently trying to write a Level 1 model for the following research question. RQ: Does the school-level relationship between race and the proportion of reading proficient students (outcome) vary by the urbanicity of the school? RACE is grouped into 7 categories and URBAN is grouped into 3 categories. I'm just unsure what the limits for the double summation are for the interaction term between race and urbanicity. There are 21 interactions. [Note: A number of control variables are included in my model as well]. Level 1: \begin{align} READPCT_{ij} &= \beta_{0j} + \sum\limits_{q=1}^7 \beta_{qj}(RACE_{qij}) + \sum\limits_{p=8}^{10}\beta_{pj}(URBAN_{pij}) + \beta_{11}(PCTFEMALE_{ij}) + \beta_{12j}(FRL_{ij})\\\ &\+ \beta_{13j}(PTR_{ij}) + \beta_{14j}(MAGNET_{ij}) + \sum\limits_{q=?}^?\sum\limits_{p=?}^?\beta_{pq}(RACE_{qij} \times URBAN_{pij}) + e_{ij} \end{align}


title: Interpreting results from a mixed model with interactions

text: I need your help interpreting my model's results. To make story short, I have a continuous dependent variable $Y$, and two factors $A$ and $B$, each with $2$ levels $(A_1, A_2, B_1, B_2)$. My main interest is the differences between $A_1$ and $A_2$, however I am taking into account the affect of $B$. The observations are not independent, in other words, each "subject" contributed more than $1$ observation, therefore, I used a mixed model with $A$, $B$ and the interaction $AB$ as fixed effects, and the "subject" as a random effect. All $3$ fixed effects were statistically significant, A with a $\text{p -value}$ of $0.0001$, $B$ with a $\text{p -value}$ of $0.0038$ and $AB$ with a $\text{p -value}$ of $0.00032$. I asked the computer to make "slicing" and to test a pair, $A_1$ vs $A_2$ for $B_1$ and for $B_2$. For $B_2$ it was statistically significant, however, for $B_1$ not despite seeing a mean difference, with a $\text{p -value}$ of $0.057$... So for $B_1$ the differences are not statistically significant (could be due to lack of power, I won't use bad phrases like "almost significant"), but overall, the factor $A$ IS statistically significant ! Now what should be my conclusion from this analysis ? What should I report and how ? For your convenience I attach the two plots that the computer gave me, will help you see my problem. ![Interatcion plot](http://i.stack.imgur.com/gtpVo.jpg) ![Means plot](http://i.stack.imgur.com/9VUuY.jpg) All the means you see are adjusted means (what SAS, NCSS and others call least square means). Again, my "goal" is to look if $A_1$ is superior over $A_2$, where lower values of $Y$ are better. Thank you in advance.


title: Why does the least square solution give poor results in this case?

text: There is an image in page 204, chapter 4 of "pattern recognition and machine learning" by Bishop where I don't understand why the Least square solution gives poor results here: ![enter image description here](http://i.stack.imgur.com/fHyoq.png) The previous paragraph was about the fact that least-squares solutions lack robustness to outliers as you see in the following image, but I don't get what's going on in the other image and why LS gives poor results there as well. ![enter image description here](http://i.stack.imgur.com/HwN5w.png)


title: Is multilevel analysis with aggregation terms the best approach for this situation?

text: My research model has two levels of dependent variable, at the organizational level and a predictor variable at individual level. I have 500 answers from 5 middle managers (in each firm) in 100 firms at an individual level and 200 from 2 top managers (in each firm) from the same 100 firms. My question is whether aggregation of each data set and using mean score in regression is the best option? If yes, what are the aggregation terms? If not what other option or software is the best to fit this model?


title: Is it valid to include several interaction terms but exclude main effects because only interaction effects are significant?

text: I have a statistical problem in connection to interaction effects and main effects in the social sciences. I am doing a regression where I have two interaction terms with one common variable. The interaction terms are used in order to prove or falsify two hypotheses. If I include all main effects and interaction effects, all effects become insignificant. My approach is to leave the main effects out. This will leave with only significant interaction effects. Is this a valid approach?


title: Pearson correlation output in SPSS

text: In my research, I am trying to show the relationship between an IV and a DV. As I compute the Pearson correlation through SPSS, it outputs a correlation matrix. This is not what I want. I need an $r$ and a $p$ value only.


title: How to calculate Zipf's law coefficient from a set of top frequencies?

text: I have several query frequencies, and I need to estimate the coefficient of Zipf's law. These are the top frequencies:               26486     12053     5052     3033     2536     2391     1444     1220     1152     1039


title: Failure detection method

text: I receive 1000 _points_ per day from installations who produces electricity. Every installation must proportionally produce the same amount of energy. I have to spot failures in those data. The actual method :   1\. I make an average.   2\. I make a second average where I remove the installation who are 30% lower than the first average   3\. I compare every installation with the filtered average. If an installation is lower than x % from the filtered average, there is a problem.   Is there some better model for that kind of purpose ? Do you have books reference about that ? Thanks in advance


title: Assessing conditional independence of genes in Trans-eQTL cluster

text: I've identified a Single Nucleotide Polymorphism (SNP) that is associated with a variable of interest (lets just call it height) using a Genome Wide Association Study (GWAS). The SNP can be of the values A/A, A/G or G/G in each of the 69 samples. I also have gene expression data for the 69 samples and it appears that the SNP is affecting the expression of a large number (~500) genes. I'm wondering if there is a convenient statistical package/technique available to assess if the variable (height) becomes independent of genotype (A/A, A/G or G/G) conditioning on the expression of each of the ~500 genes who's expression are associated with the SNP (i.e. assessing conditional independence of each of the ~500 genes)?


title: Automatic detection of level changes in series of prices

text: I have a large number of time series which consist of pricing data of consumer goods. As expected the prices show trend and seasonality. However my main problem is to detect large level changes in the price series. In these data sets this has occurred mainly due to the price being recorded for the wrong packaging size. For example instead of the price for a can of beer being recorded, the price for a six-pack has been recorded. For any given time series it is likely the case that most prices are correct, but a series of observations may be wrong. Most likely the errors are at the beginning of the series but not necessarily. This is easy to pick up when the time series is plotted and can be modeled with a dummy variable. However I would like to automatically detect the location and magnitude of the level change. How can I go about this ? Unfortunately my stats education ended at the undergrad level so I'm not sure where to begin ? The trend and seasonality that exists in my time series is not my main concern. Do I need to worry about autocorrelation etc ? Or should I just be worried about level changes. Although I have access to R, ultimately the algorithm may have to be implemented in Java. I'm not sure if this is an appropriate place for the question, but I hope someone can help me !!


title: Convergence of neural network weights

text: I came to a situation where the weights of my Neural Network are not converging even after 500 iterations. My neural network contains 1 Input layer, 1 Hidden layer and 1 Output Layer. They are around 230 nodes in the input layer, 9 nodes in the hidden layer and 1 output node in the output layer. I wanted to know, if i do early stopping condition (say stopping my neural network training after 100 iteration) What effect the model will have on this. Also wanted to know what is the industry standard of work around if the weights in the neural network are not converging?


title: What is the practical significance of a linear regression, lasso regression, and ridge regression outputting the same coefficients?

text: The lasso and ridge regression are tuned to the same alpha parameter. No matter what I tune the parameter to [0,1], the results of all three regressions are always the same (linear, ridge, lasso), which is a bit befuddling because I had previously suspected that there was some multicollinearity in my data due to high variance inflation factors. I thought that either the lasso or ridge would address the issue of multicollinearity, but this does not seem to be the case. If anyone can provide some insight that would be great. Thanks.


title: Confusion related to EM algorithm

text: I was reading this tutorial related to EM algorithm at http://aass.oru.se/~tdt/ml/extra-readings/EM_algorithm.pdf. As given in the tutorial ![enter image description here](http://i.stack.imgur.com/v5bqe.png) we can see that at each E step we calculate the expectation of the likelihood over the posterior distribution of the hidden variables and them maximize it.The expected log likelihood over the posterior distribution hidden variables is upper bounded by the likelihood over the observed data only. My questions are:   * Why is $L(\theta)$ non convex? Isn't it possible for $L(\theta)$ to be convex?   * Further why is $l(\theta | \theta_n)$ a convex function as shown in the figure?


title: Linear Regression with $L_2$: Different penalty strengths yield the same parameters?

text: Suppose I have model 1: $$Y=aX$$ where $X$ is $n\times1$, consisting of a single feature. suppose I fit this model with $L_2$ penalty with coefficient $\lambda=1$ : $$ \underset{a}{\min} \sum_{i=1}^{n}(y_i-ax_i)^2 +a^2$$ and I find that $a=b$. Now consider model 2 : $Y=a'X'+a''X''$ where $X',X''$ are just duplicates of $X$ so $X'=X''=X$. If we fit this model with $L_2$ penalty using $\lambda=1$, $$ \underset{a',a'' }{\min} \sum_{i=1}^{n}(y_i-a'x'_i-a''x''_i)^2 +a'^2 +a''^2$$ We find that $a'=a''=\frac b2$. If we set $\lambda=2$, we get $$\underset{a',a'' }{\min} \sum_{i=1}^{n}(y_i-a'x'_i-a''x''_i)^2 +2(a'^2 +a''^2) = \min_{a' } \sum_{i=1}^{n}(y_i-2a'x_i)^2 +4a'^2$$ since $X'=X''$ and $a'=a''$. Letting $c=2a'$, this reduces to $$\underset{c}{\min} \sum_{i=1}^{n}(y_i-cx_i)^2 +c^2$$ which is the same optimization as in model 1, which yields $c=b$ and thus $a'=a''=\frac b2$ - same parameters even though the penalty strength is different. This seems to conflict with the analytic solution for Ridge regression: $$a_{ridge} = (X^TX+\lambda I)^{-1}X^TY$$ where different $\lambda$ lead to different parameters. How do you reconcile this?


title: Point Biserial correlation coefficient in Java

text: I've recently installed Apache Commons Math library to help me with some basic statistics in my project. Basically, what I need is to calculate the correlation coefficient for 2 arrays, which I can do using the class "PearsonsCorrelation". However, I realized that I should be using a special variant of Pearson when one of the variables is nominal - Point Biserial correlation coefficient. Could you tell me if there is a similar Java library which can calculate the Point Biserial coefficient, since I had no luck in my search. Thanks for your help! Best regards,


title: model noise with Gaussian(0, Gamma)

text: Can anyone help me in understanding what kind of noise is it? Noise = Normal(0, v) v = GammaFromShapeAndRate(alpha, beta) I mean what is the advantage of making a normal noise with a Gamma distribution? And if it shows any specific kind of noise, then how should I determine alpha and beta? should I set them in my graphical model? or I should make them other hyperparameters comming from normal distribution and try to learn them. Thanks.


title: Learning Conditional Random Fields using EM Algorithm (from unaligned data)

text: I am trying to learn CRF from the unaligned data in Natural Language Understanding application. There is one paper in this field which does exactly the same, Learning conditional random fields from unaligned data for natural language understanding. In this paper, the authors mention that first log-likelihood $L(\theta)$ is calculated with expectation over the abstract annotation using this equation: $ L(\theta; \theta^t) = \sum_{i}^{M}\sum_{Y_{i}^{u}}P(Y_i^u|X_i;\theta)\log P(Y_i^u|X_i;\theta) =\sum_{i}^{M}\sum_{Y_{i}^{u}}P(Y_i^u|X_i;\theta)\sum_t\sum_k\theta_kf_k(y',y,X_i)-\sum_i^k\log Z(X_i) $ where $Y_i^u$ is the unknown semantic tag sequence of $i$-th word sequence, and $Z$ is the normalization factor. This equation is somehow like auxiliary $Q(\theta; \theta^t)$ function in the EM algorithm; however, there are some differences. On the other hand I cannot find any relation between this equation and HCRF objective function on the hidden variables. Can anyone give any further explanations or comments on this paper and specifically the mentioned equation? I am going to implement this paper and need to know exact details on this equation?