title: Overall standard deviation

text: I want to determine the standard deviation for a set of measurements of a cable profile. However, the cable follows the terrain, and thus using a mean value for the entire cable length does not make any sense. My current thought is to compute the mean and standard deviation at every point based on the point itself and the two previous and the two next points. Then I can construct a confidence interval for the standard deviation, via the chi-squared distribution, at each of these points. How can I get a single meaningful overall standard deviation and a confidence interval for it for the entire cable?


title: Can your manipulation condition be a predictor in a multiple regression?

text: I recently ran an experiment which manipulated participants in one of two conditions, by either promoting similarities between two groups, or emphasizing differences between two groups. I then ran a multiple regression to find out what predicts the amount of tax that a participant would be willing to pay for another group. As the analysis was exploratory, I added all of the variables including my nominal manipulation condition (group 1 = promote similarities, group 2 = emphasize differences). The analysis determined that Symbolic threat (threat felt from cultural values of others) and condition (which manipulation was assigned to participant, i.e which group individual belonged to) predicted the DV (amount will pay for other group in tax). I was just wondering, is it weird and unusual to add your manipulation condition as a predictor in a multiple regression? I made a scale which was supposed to measure the manipulation and it differed significantly between groups. However, when doing a MR this is not included as a predictor... When using the manipulation condition (group 1 or 2) it is included. I don't really understand why this is. Any advice would be really appreciated. Thank you for taking the time to read this, and for any helpful comments. Richard


title: Null hypothesis of the runs test

text: From Wikipedia > The runs test (also called Wald–Wolfowitz test after Abraham Wald and Jacob > Wolfowitz) is a non-parametric statistical test that checks a randomness > hypothesis for a two-valued data sequence. More precisely, it can be used to > test the hypothesis that the elements of the sequence are mutually > independent.   1. I was wondering if its null hypothesis is that the elements of the sequence are iid? (note that the quote only mentions mutual independence.)   2. What if the sequence is three-valued? Thanks and regards!


title: Test for finite variance?

text: Is it possible to test for finiteness (or existence) of the variance of a random variable given a sample? As a null, either {the variance exists and is finite} or {the variance does not exist/is infinite} would be acceptable. Philosophically (and computationally), this seems very strange because there should be no difference between a population without finite variance, and one with a very very large variance (say > $10^{400}$), so I am not hopeful this problem can be solved. One approach that had been suggested to me was via the Central Limit Theorem: assuming the samples are i.i.d., and the population has finite mean, one could check, somehow, whether the sample mean has the right standard error with increasing sample size. I'm not sure I believe this method would work, though. (In particular, I don't see how to make it into a proper test.)


title: Misunderstanding of use of standard deviation across groups

text: I'm trying to understand standard deviation. I understand how to determine standard deviation across a group, but I am not understanding using standard deviation to compare groups. Am I attempting to use the correct tool? For example, group 1, 2011, for 10 people, these are the results:               99 98 91 89 88 74 74 72 67 56     average = 81     Variance = 18 17 10 8 -7 -7 -9 -14 -25          324 289 100 64 49 49 81 196 625          1777/10 = 178 so stdev = 13      Group 2, 2012, for 10 people, these are the results:               99, 91, 90, 88, 87, 86, 84, 83, 74, 65     average = 85     variance = 80     stdev = 9      1) If we assume that the tests were similarly hard, but may have slightly different questions, would it be correct to assume both students with 99 have the same ability? Or should I be comparing apples to apples by adjusting both sets of scores to a set average, say, 80? If I need to adjust, what do I do to compare fairly, and what is it called?


title: Finding correlations in longitudinal data analysis

text: We are doing research on video lecture watching. We offer a course which lasts four 5 weeks, and there are 9 or 10 videos in each week. We organize group- watching activities. A group usually is composed of 4-5 people, and we have 3 groups. In each group, each participant watched all videos of that week on his/her own device with our software (so we could collect data). We collected data on how they navigated the videos. Right after each week's video watching session, we asked each person to rate the difficulty level of the whole weeks' lecture with 5 point Likert scale. After 5 weeks, we get a dataset which includes the number of pauses per user per week and the difficulty rating per user per week. The data set is organised as follows: For the Likert-scale questionnaire, the data is like                 group, person, video.difficultyOfTheWeek, week       apricot,   A,     5,                       1       apricot,   B,     3,                       1       apricot,   C,     4,                       1       apple,     A,     3,                       1       apple,     B,     2,                       1       apple,     C,     2,                       1       orange,    A,     4,                       4       orange,    B,     3,                       4       orange,    C,     4,                       4      We also have pause data similar as follows:                 group, person, numOfPauses, week, totalLengthOfVideoInTheWeekInMinute       apricot,   A,    15,         1,             125       apricot,   B,    23,         1,             125       apricot,   C,    24,         1,             125       apple,     A,    13,         1,             125       apple,     B,    12,         1,             125       apple,     C,     8,         1,             125       orange,    A,    11,         4,             156       orange,    B,     4,         4,             156       orange,    C,     9,         4,             156      What I want to do for the next step is to answer the following question. **Does the number of pauses correlate with the perceived difficulty?** I tried to normalize the number of pauses. I composed a new data set by introducing a new variable "pauseFrequency", which is computed by dividing the number of pauses with the total length of the video of that week (same for all users in that week).                 group, person,  difficulty, week,      pauseFrequency       apricot,   A,     5,         1,             15/125       apricot,   B,     3,         1,             23/125       apricot,   C,     4,         1,             24/125       apple,     A,     3,         1,             13/125       apple,     B,     2,         1,             12/125       apple,     C,     2,         1,             8/125       orange,    A,     4,         4,             11/156       orange,    B,     3,         4,             4/156       orange,    C,     4,         4,             9/156      Then the problem seems to be easy. It seems that I just need to make correlation test with the difficulty column and the pauseFrequency column. I did it. I actually treat all difficulty/frequency pairs in the same way, no matter they are from the same group or from in the same week or whatever. I treat them as individual observation. I am using R for analysis, then I did a Kendall's Rank Correlation test like the following:               cor.test(pauseFrequency,video.difficulty,method="kendall")      This of course has generated a result, but I am not confident this is correct.I have two concerns:   1. My measures are actually longitudinal and repeated, each user in each group is measured 5 times.   2. The "task" for each week is actually different. Although they watched videos of the same course for 5 weeks, but each week they watch a different series of videos (with different content and lengths) The observations are actually not independent. How can I compensate for this non-independence in my data analysis?


title: p-value adjustment for Local Moran's I statistic (LISA)

text: I'm working with some exploratory spatial analysis in R using spdep package. I came across an option to adjust _p_ -values of local indicators of spatial association (LISA) calculated using `localmoran` function. According to the docs it is aimed at: > ... probability value adjustment for multiple tests. Further in the docs of `p.adjustSP` I read that the options available are: > The adjustment methods include the Bonferroni correction ('"bonferroni"') in > which the p-values are multiplied by the number of comparisons. Four less > conservative corrections are also included by Holm (1979) ('"holm"'), > Hochberg (1988) ('"hochberg"'), Hommel (1988) ('"hommel"') and Benjamini & > Hochberg (1995) ('"fdr"'), respectively. A pass-through option ('"none"') is > also included. > > The first four methods are designed to give strong control of the family- > wise error rate. There seems no reason to use the unmodified Bonferroni > correction because it is dominated by Holm's method, which is also valid > under arbitrary assumptions. > > Hochberg's and Hommel's methods are valid when the hypothesis tests are > independent or when they are non-negatively associated (Sarkar, 1998; Sarkar > and Chang, 1997). Hommel's method is more powerful than Hochberg's, but the > difference is usually small and the Hochberg p-values are faster to compute. > > The "BH" (aka "fdr") and "BY" method of Benjamini, Hochberg, and Yekutieli > control the false discovery rate, the expected proportion of false > discoveries amongst the rejected hypotheses. The false discovery rate is a > less stringent condition than the family-wise error rate, so these methods > are more powerful than the others. Couple of questions that appeared:   1. In plain words - what is the purpose of this adjustment?    2. Is it necessary to use such corrections?   3. If yes - how to choose from available options?


title: Machine learning predicted value

text: When we fit a generalized linear regression (e.g., logistic regression, gamma regression) we are estimating the population average Y given the predictors $X$ ( i.e., $E(Y | X)$ ). When we fit a machine learning model such as an ANN, SVM, or a decision tree, does this notion still apply? In other words, are we estimating the population average value of $Y$ or isn't that idea applicable and we are just predicting "Y"? **ADD After Dikran response:** I. What aspect of the theory of a predictive modeling algorithm tells us that we are modeling E(Y|X) versus just Y|X ? Is it the use of an error term that follows a certain distribution? For instance, what is it about ANN versus a decision tree which tells us the former models E(Y|X) while the latter is modeling Y|X? II. Is there any connection between these and say a confidence interval versus prediction interval in linear regression?


title: Using Anselin Local Moran's I Values in Regression

text: I am doing a multiple linear regression of factors related to poverty and would like to include some spatial-statistical data. I have come up with Anselin Local Moran's I values (cluster/outlier) data for census tracts in a metropolitan area, which I would like to include in the regression. However, the Anselin Local Moran's I values are indexed with either positive or negative values. A positive value for I indicates that a feature has neighboring features with similarly high or low attribute values; this feature is part of a cluster. A negative value for I indicates that a feature has neighboring features with dissimilar values; this feature is an outlier. These values are, however, only considered significant if indicated so by their corresponding z or p scores. Is it possible to use a negative/positive index like this as an independent variable in a regression? If so, how could this problem be approached considering only some of the observations (Moran's I values for census tracts) are considered significant? More details about Anselin Local Moran's I can be found here: http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//005p00000012000000


title: Why is the ROC curve of a random classifier the line $x=y$?

text: The title is my whole question.


title: Analysis of a dataset and tool to get the incidence of an event

text: I have a dataset and I want to do statistical analysis:               session    day     type   entity   name           RunningTime  PLACE     ----------------------------------------------------------------------      1      29/09/10   1      Tera     NikolaT        24:33        unknown      2      14/12/10   1      Vera     EmmanuelF      22:13        Aldebaran      3      29/01/11   1      Eon      Eon            null         all          4      01/03/10   1      NA       DSG            12:25        earth      5      02/11/10   1      Ardan    alvaroR        43:56        Aldebaran      for instance I want to get the probability distribution of the incidence of an event, ocurrence of sessions...., I was planing first to sort by day and then getting $\mu$ and $\sigma$ for each month and then by year, but what is the best approach to get valuable information from this dataset?   what statistical analysis can be done in a dataset like this?   what tool instead of matlab or R can I use, that is online?


title: Dealing with dependant data when estimating probability of an event happening

text: I have 10 years worth of data from 1970 to 1980 (40 quarters). For each quarter I have five measurements M1, M2, M3, M4 and M5. **TWIST:** Although the data I have is on individual patient level, the outcome variable (T) is common to all patients during epidemic times. For example, all patients tend to have their target flags set to 1 (positive) in quarters Q19, Q20 and Q21 as well as in quarters Q33, Q34, Q35, Q36, Q37, Q38, Q39 and Q40. The reason for this is because patients are flagged with 1 even if they have mild conditions of a particular disease and they all tend to have it during epidemic times. The common target variable makes me think that I should aggregate my measurements on from individual to population level and then estimate probability of widespread disease occurring. However, that means that I’m “losing” information i.e. I’m ignoring potentially relevant information on an individual level. My “gut feeling” is that I shouldn’t be aggregating and averaging my data. On the other side if I use a simple logistic regression model and treat my data as if it was independent (which it clearly isn’t) I’ll probably get biased estimates. My ultimate goal is to be able to (given individual measurements) calculate the probability for that particular individual of being ill.


title: Relative advantages of multiple imputation and expectation maximization (EM)

text: I've got a problem where $$y = a + b $$ I observe y, but neither $a$ nor $b$. I want to estimate $$b = f(x) + \epsilon$$ I can estimate $a$, using some sort of regression model. This gives me $\hat b$. I could then estimate $$\hat b = f(x) + \epsilon$$ First problem: a regression model to predict $a$ could lead to $\hat b$ being negative, which wouldn't make any sense. Not sure how to get around this (not the sort of problem I've dealt with a lot) but seems like the kind of thing that others deal with routinely. Some sort of non-gaussian GLM? The main problem is how to account for the uncertainty in the main model that comes from estimating $\hat b$. I've used multiple imputation before for missing covariates. But this is a missing "latent parameter." Alternatively, it is outcome data, which seems OK to impute. However I often hear of EM used for "latent" parameters. I am not sure why, nor do i know whether EM is any better in these contexts. MI is intuitive both to understand, implement, and communicate. EM is intuitive to understand, but seems more difficult to implement (and I haven't done it). Is EM superior for the sort of problem I've got above? If so, why? Second, how does one implement it in R for a linear model, or for a semiparametric (GAM) model?


title: Significance test for classification systems

text: I have several classification systems which I want to compare. The way I evaluate each one of them is using precision, recall and f-score. I want to compare the performance of every two of them using some significance test to find out if the difference in the performance is due to noise or not. My question is - what would be a good choice for significance test, given the fact that I want to use the f-score values as an input. So far, most of the tests I have seen compare means w.r.t. standard deviation, but I am not sure this is what I need in my case.


title: Plotting two densities on the same graph in R

text: I generated two density functions in R and trying to plot them on the same graph, but for some reason I cannot see the full plots. Here is the R code:               v3 <- rt(100000, 1)/sqrt(3-2)      w3 <- rchisq(100000,2)      z3 <- rnorm(n=100000, m=0, sd=1)      z_eff_3 <- v3 + z3 * sqrt((3*(1+v3*v3))/w3)     plot(density(z_eff_3,from=0,to=5),xlim=c(0,5))          v4 <- rt(100000, 1)/sqrt(4-2)      w4 <- rchisq(100000,3)      z4 <- rnorm(n=100000, m=0, sd=1)      z_eff_4 <- v4 + z4 * sqrt((4*(1+v4*v4))/w4)          lines(density(z_eff_4,from=0,to=5),xlim=c(0,5),col ="blue")      Do you guys know how to fix this?


title: How many correlations are too many?

text: Following my earlier posts here, here and here, I understand the dangers of 'fishing' for a significant correlation. I also understand the 'problem' of running several correlations on the same data set. However, I am testing the relationship between variables A and B by several demographic characteristics, so I cannot help but run several correlation tests. For example, on gender I have 2 correlation tests (one for male and one for female), on age I have 3 correlation tests (one for less than 18 years, one for 18 to 44 years and one for over 45 years), on marital status, I have 3 correlation tests (one for married, one for never married and one for other). **With all the concerns expressed in the earlier posts, is there anything I am missing on the use of Pearson r?** Remember, I am just focusing on the relationship between A and B by demographic characteristics of the respondents. **Is there an alternative way to test this relationship?**


title: How can I predict the odds that a dodgeball team is going to win based on the winning history of its players?

text: Imagine there are 80 dodgeball players in the world. Each of them has played thousands of dodgeball games with the other 79 players in more-or-less random order. This is a world without teams (e.g., every player has a chance of being drafted in either team each game). I know the previous win rate of each player (e.g., one has won 46% of all previous games, another has won 56% of all his previous games). Lets say there is a match coming up and I know who is playing on each team. I also know their previous win rate. What is the best way to calculate the probability of each team winning based on the composition of the team? If it requires relatively advanced calculation (e.g., logistic regression) let me know some of the specifics. I am pretty familiar with SPSS, but I rather not need to ask a follow-up question. Moreover, how would I explore the accuracy of my method using archival data? I know it won't be clear cut since most players hover around 40-60%, but still. To be specific, what are the odds that team A is going to win? A - comprised of individuals with previous win rate of 52%, 54%, 56%, 58%, 60% B - comprised of individuals with previous win rate of 48%, 55%, 56%, 58%, 60% (this is just a random example for illustrative purposes. Two pretty good teams.) Edit: Is there a way to start with a very simple algorithm and then see how it works? Maybe we could simply sum the percentages of each team and predict that the one with the highest percentage is going to win. Of course our classification would not be accurate, but over thousands of archived games we could see if we can predict better than chance.


title: Problem printing greyscale or B&W ggplot2 images

text: I have an issue with printing/creating black and white and greyscale images with ggplot2. I have been trying to make some black and white, and greyscale graphics for publication, they look perfect on mine and everyone's monitors, but came out in colour when printed. (An awful photo of the result is here: http://i.imgur.com/UuB6A.jpg) For example:               p <- qplot(displ, hwy, data = mpg)     p <- p + facet_wrap(~ cyl) + theme_bw()           pdf("testplot.pdf",width=7.6,height=3.6)          p      dev.off()      The generated image and output pdf looks perfect onscreen, but the facets borders and titles are in a shade of red when printed. Try it! I get the same results when outputting to a .png, and also when I print directly from the R graphics window. According to my friend with some experience with image manipulation, the images I've sent him have bits of colour in them according to the colour histogram. Is there any way to generate true B&W and greyscale images using ggplot? I tried forcing greyscale (perhaps very naively!) by using 'colormodel="gray"' in the pdf() command, but it does not make a difference. Any help would be greatly appreciated.


title: Do Julian day number handle timezone offset?

text: I am currently working on some FX rate modeling. I decided to use the dataset from UBC: http://fx.sauder.ubc.ca/data.html But I find that I have a problem about the timezone. The dataset is daily basis, with a Julian Day Number, does this number handle timezone shift? Or does it use GMT+0 or something else?


title: Wavelets and machine learning

text: I am trying to learn features from a signal using Wavelet transform and then apply ML techniques on it to classify a signal. The problem I am facing is that, at each of level of decomposition, my signal has a vector, and I have more than 10,000 signals. Or let me rephrase: Suppose I want to classify if a house will sell or not. Each house has a lot of features that let me decide. In my feature matrix, each row represents a single house. And in each row, each column represents a random variable pertaining to size, lawn, floorplan etc.. But in the case of Wavelets, I have a signal which I have decomposed. Every level of decomposition has a vector associated with it. And there are many signals. Now in my feature matrix, every row represents a single signal. And each rows' each column is a random variable which is supposed to be a scalar, but in this case, it is a vector. Does anyone know how to tackle this problem? I can't do averaging or mean or suppose it belongs to any distribution.


title: Can you compare p-values of Kolmogorov Smirnov tests of normality of two variables to say which is more normal?

text: I have applied the one sample Kolmogorov Smirnov test of normality to two variables and one has a larger p value but both are greater than .05. e.g.,   * $x_1$ (p-value) = 0.09    * $x_2$ (p-value) = 0.06  Does this mean that $x_1$ is better or more normal than $x_2$?


title: Using linear discriminant analysis to validate the cluster groups resulting from kmeans

text: I'm currently working on a cluster analysis project and ran kmeans on the data for k=2. I was reading similar articles on similar experiments, and the investigators used discriminant analysis to verify their cluster groups. I was wondering how to use linear discriminant analysis is used to verify the cluster groups from kmeans, is LDA done separately and then done compared to the kmeans result or if a result from the kmeans is used as a parameter in the LDA call in R? I'm sorry if this question is vague, but I have never been exposed to discriminant analysis before. Also the data gathered is repeated measures, and unbalanced per individual, not sure if that valids some assumptions for LDA.


title: Algorithm for a 5 star rating system

text: I'm creating a reporting tool that produces a data set and calculates a metric for each item in the data set. I have two algorithms that I'm working with: an algorithm that produces the metric and an algorithm that produces a rating(1-5 star rating) for visual representation of the calculated metric. The algorithm to calculate the Metric will be fine-tuned on an ongoing basis. It's a prediction tool, so I will constantly compare the metric to real data and modify the algorithm accordingly. On the other hand, I'd like the Rating algorithm to remain somewhat consistent. My initial approach to the Rating algorithm was to use a percentile function similar to Excel's to find the value at: 98th percentile(5 stars), 90th(4 stars), 75%(3 stars), etc. Here is the function I used. (This is the Ruby programming language, but it should be readable enough to get the jist of it)               def percentile(values, percentile)         values_sorted = values.sort         k = (percentile*(values_sorted.length-1)+1).floor - 1         f = (percentile*(values_sorted.length-1)+1).modulo(1)              return values_sorted[k] + (f * (values_sorted[k+1] - values_sorted[k]))       end      If I understand this equation correctly, it's producing a percentile value based on the number of items rather than the actual values of the items. Typically, I will have a data set with 150 items. This is a poor method of calculating a Rating because it will always produce exactly 3 five star items. I need a better approach that will calculate a Rating based on the actual Metric values. Any ideas? Here are a few sample data sets: Sample 1 - 150 items               2,466,505     1,118,729     935,150     431,367     271,558     243,603     168,266     166,507     164,543     137,216     132,801     120,329     109,000     103,088     101,341     89,685     80,173     74,703     70,748     66,970     63,408     60,368     54,433     54,184     51,760     46,812     45,659     43,140     42,682     41,864     40,151     39,989     35,766     34,020     33,379     32,803     29,393     28,791     28,570     27,724     27,341     26,282     24,431     23,537     22,804     22,634     22,217     21,812     21,648     21,230     20,657     20,312     19,602     18,157     17,924     17,398     16,771     15,359     15,238     15,156     15,156     15,126     15,068     14,654     14,300     14,109     13,924     13,893     12,813     12,612     12,419     12,236     11,736     11,586     11,367     10,349     10,297     9,988     9,900     9,112     8,453     8,420     7,912     7,893     7,525     7,484     7,440     7,429     6,736     6,518     6,497     6,174     5,613     5,473     5,242     5,067     4,657     4,262     4,152     4,104     4,069     4,052     4,052     4,045     3,990     3,952     3,814     3,621     3,387     3,280     3,280     3,189     2,806     2,804     2,789     2,701     2,701     2,540     2,476     2,320     2,320     2,227     2,129     2,064     2,020     1,989     1,951     1,878     1,783     1,713     1,641     1,590     1,341     1,341     1,341     1,326     1,248     1,105     1,030     912     723     609     482     368     314     237     74     65     34     28      Sample 2 - 150 items                3,553,644       723,572       356,706       321,258       184,440       173,914       133,995       117,048       101,442       87,061       80,692       62,248       39,193       32,277       31,656       29,794       27,932       27,666       26,792       26,604       25,438       23,459       22,612       22,346       22,249       21,104       20,949       19,508       18,533       17,380       16,582       16,582       16,138       14,897       14,631       14,631       13,996       13,656       12,985       12,104       12,104       11,731       11,627       10,729       10,182       9,822       9,154       9,062       8,876       8,415       8,415       8,069       7,819       7,645       7,200       7,162       6,840       6,172       6,172       5,657       5,574       5,544       5,303       5,143       5,089       4,946       4,946       4,895       4,895       4,895       4,895       4,685       4,564       4,357       4,256       4,180       3,752       3,671       3,628       3,482       3,386       3,339       3,104       3,009       2,979       2,979       2,913       2,902       2,554       2,451       2,314       2,217       2,147       1,951       1,951       1,874       1,862       1,756       1,664       1,643       1,561       1,561       1,561       1,561       1,561       1,561       1,483       1,426       1,419       1,419       1,345       1,341       1,330       1,330       1,277       1,234       1,229       1,194       1,170       1,170       1,158       1,149       1,149       1,146       1,135       1,117       1,117       1,064       958       859       851       851       766       721       681       638       585       575       575       479       383       383       319       284       276       234       211       195       44       44       Sample 3 - 58 items               60,349     43,949     31,562     24,505     17,767     16,928     14,801     11,564     10,695     9,567     8,883     8,712     7,572     7,468     7,340     7,302     6,638     6,472     6,223     5,974     5,808     5,744     5,476     5,393     3,717     3,566     3,556     3,368     3,242     3,126     2,863     2,626     2,613     2,555     2,409     2,391     2,351     2,309     2,125     2,042     1,936     1,825     1,797     1,659     1,521     1,475     1,294     1,171     1,120     1,045     878     757     631     543     461     431     426     237


title: Collinearity diagnostics problematic only when the interaction term is included

text: I've run a regression on U.S. counties, and am checking for collinearity in my 'independent' variables. Belsley, Kuh, and Welsch's _Regression Diagnostics_ suggests looking at the Condition Index and Variance Decomposition Proportions:               ## colldiag(, scale=TRUE) for model with interaction     Condition     Index   Variance Decomposition Proportions                (Intercept) inc09_10k unins09 sqmi_log pop10_perSqmi_log phys_per100k nppa_per100k black10_pct hisp10_pct elderly09_pct inc09_10k:unins09     1    1.000 0.000       0.000     0.000   0.000    0.001             0.002        0.003        0.002       0.002      0.001         0.000                 2    3.130 0.000       0.000     0.000   0.000    0.002             0.053        0.011        0.148       0.231      0.000         0.000                 3    3.305 0.000       0.000     0.000   0.000    0.000             0.095        0.072        0.351       0.003      0.000         0.000                 4    3.839 0.000       0.000     0.000   0.001    0.000             0.143        0.002        0.105       0.280      0.009         0.000                 5    5.547 0.000       0.002     0.000   0.000    0.050             0.093        0.592        0.084       0.005      0.002         0.000                 6    7.981 0.000       0.005     0.006   0.001    0.150             0.560        0.256        0.002       0.040      0.026         0.001                 7   11.170 0.000       0.009     0.003   0.000    0.046             0.000        0.018        0.003       0.250      0.272         0.035                 8   12.766 0.000       0.050     0.029   0.015    0.309             0.023        0.043        0.220       0.094      0.005         0.002                 9   18.800 0.009       0.017     0.003   0.209    0.001             0.002        0.001        0.047       0.006      0.430         0.041                 10  40.827 0.134       0.159     0.163   0.555    0.283             0.015        0.001        0.035       0.008      0.186         0.238                 11  76.709 0.855       0.759     0.796   0.219    0.157             0.013        0.002        0.004       0.080      0.069         0.683                      ## colldiag(, scale=TRUE) for model without interaction     Condition     Index   Variance Decomposition Proportions                (Intercept) inc09_10k unins09 sqmi_log pop10_perSqmi_log phys_per100k nppa_per100k black10_pct hisp10_pct elderly09_pct     1    1.000 0.000       0.001     0.001   0.000    0.001             0.003        0.004        0.003       0.003      0.001             2    2.988 0.000       0.000     0.001   0.000    0.002             0.030        0.003        0.216       0.253      0.000             3    3.128 0.000       0.000     0.002   0.000    0.000             0.112        0.076        0.294       0.027      0.000             4    3.630 0.000       0.002     0.001   0.001    0.000             0.160        0.003        0.105       0.248      0.009             5    5.234 0.000       0.008     0.002   0.000    0.053             0.087        0.594        0.086       0.004      0.001             6    7.556 0.000       0.024     0.039   0.001    0.143             0.557        0.275        0.002       0.025      0.035             7   11.898 0.000       0.278     0.080   0.017    0.371             0.026        0.023        0.147       0.005      0.038             8   13.242 0.000       0.001     0.343   0.006    0.000             0.000        0.017        0.129       0.328      0.553             9   21.558 0.010       0.540     0.332   0.355    0.037             0.000        0.003        0.003       0.020      0.083             10  50.506 0.989       0.148     0.199   0.620    0.393             0.026        0.004        0.016       0.087      0.279              `?HH::vif` suggests that VIFs >5 are problematic:               ## vif() for model with interaction             inc09_10k           unins09          sqmi_log pop10_perSqmi_log      phys_per100k      nppa_per100k       black10_pct        hisp10_pct               8.378646         16.329881          1.653584          2.744314          1.885095          1.471123          1.436229          1.789454          elderly09_pct inc09_10k:unins09               1.547234         11.590162           ## vif() for model without interaction             inc09_10k           unins09          sqmi_log pop10_perSqmi_log      phys_per100k      nppa_per100k       black10_pct        hisp10_pct               1.859426          2.378138          1.628817          2.716702          1.882828          1.471102          1.404482          1.772352          elderly09_pct               1.545867       Whereas John Fox's _Regression Diagnostics_ suggests looking at the square root of the VIF:               ## sqrt(vif) for model with interaction             inc09_10k           unins09          sqmi_log pop10_perSqmi_log      phys_per100k      nppa_per100k       black10_pct        hisp10_pct               2.894589          4.041025          1.285917          1.656597          1.372987          1.212898          1.198428          1.337705          elderly09_pct inc09_10k:unins09               1.243879          3.404433      ## sqrt(vif) for model without interaction             inc09_10k           unins09          sqmi_log pop10_perSqmi_log      phys_per100k      nppa_per100k       black10_pct        hisp10_pct               1.363608          1.542121          1.276251          1.648242          1.372162          1.212890          1.185108          1.331297          elderly09_pct               1.243329       In the first two cases (where a clear cutoff is suggested), the model is problematic only when the interaction term is included. The model with the interaction term has until this point been my preferred specification. I have two questions given this quirk of the data:   1. Does an interaction term always worsen the collinearity of the data?   2. Since the two variables without the interaction term are not above the threshold, am I ok using the model with the interaction term. Specifically, the reason I think this might be ok is that I'm using the King, Tomz, and Wittenberg (2000) method to interpret the coefficients (negative binomial model), where I generally hold the other coefficients at the mean, and then interpret what happens to predictions of my dependent variable when I move `inc09_10k` and `unins09` around independently and jointly.


title: Code binary categorical (participates: yes/no) as a continuous variable (proportion of cells within participating country) in linear regression model?

text: I have a question regarding the coding of variables in r for linear regression models: I did regression analyses to find out the relative importance of 12 variables in predicting different continuous response variables. I want to use different types of linear regression models (OLS, GLMs with binomial distribution, simultaneous autoregressive models). All predictor variables are continuous, except for one which describes whether or not a country is participating in an international programme. It is essentially a factor with binary coding (country participates: 1/0). I ran the models at different spatial resolutions (i.e. differently sized grid cells). Because the coarser grids have quite a few cells that overlay with participating an non-participating countries, instead of coding the variable as binary, I calculated the proportion of every grid cell that lies within a participating country. I arcsine-square-root-transformed the resulting values (between zero and one) and treated them as continuous in my regression models. I would now interpret them as the degree to which a grid cell is participating. The problem is that while the scatterplots look OK at coarse resolutions, at fine resolutions only a small proportion of grid cells really lies at the boarders between countries, so the vast majority of grid cells are still coded as zeros and ones. Accordingly, the scatterplot looks quite awkward. Can anyone explain to me which (if any) assumptions I violate for OLS, binomial GLM, SAR, or in general? And of course whether it is a major problem or whether I can still get away with it (or would reviewers tear me apart) ? I've already run all my models and it took quite some time, so I'd be a bit hesitant to redo everything if that's not absolutely necessary. Thanks a lot for your help!! Carsten


title: How does centering the data get rid of the intercept?

text: I keep reading about instances where we center the data (e.g., with regularization or PCA) in order to remove the intercept (as mentioned in this question). I know it's simple, but I'm having a hard time intuitively understanding this. Could someone provide the intuition or a reference I can read?


title: Multiple Linear Regression Simulation

text: I'm new to the R language. I would like to know how to simulate from a multiple linear regression model that fulfills all four assumptions of the regression.


title: What does changing the time-frame of the dependent variable tell me?

text: I have a question concerning a cross sectional regression model. I want to know how to interpret the change in the explanatory variables gathered from a base year (say 2005) when I replace my observed dependent variable form 2005 to values observed in another year (2010). So, if my initial model looks like this: Y (2005)= β1 + β2X(2005) + β3X(2005)…. With the observed variables all constructed with data from the year 2005. What would happen if I only change my dependent variable Y(2005) to Y(2010) (for example observed house prices year 2005 and 2010)? Leaving the others as they are. Can I interpret it as change in the explanatory variables effect over time? Is this a common approach?   I appreciate any help I can get, Thank you


title: Meaning of Qqnorm plot in R

text: I am testing the normality of a sample with R using qqnorm. I obtain this: ![Normal Q-Q plot](http://i.stack.imgur.com/lv2Di.png) I understand that the meaning of this plot is that the sample has fat tails. But what is the meaning of the values on the $x$ axis? Are they standard deviations? Are they quantiles? If they are standard deviations, I understand that more that a 99% of points follows quite nicely a normal distribution. If they are quantiles, I can only say that for 2/3 of the points. However, it looks really nice here: ![Histogram vs Normal distribution](http://i.stack.imgur.com/OW19H.png)


title: Determining the confidence interval of Monte Carlo data

text: I want to determine the confidence interval for my set of data. I have obtained the data by sampling from several Normal Distributions and running a Monte Carlo Simulation. I was wondering how I could determine the confidence interval interval for the $\mu$ of my data? I don't know the $\mu$ or $\sigma$ of the population data. Also, is there a way by which I can determine the % confidence associated with my sampled data being a good representation of the overall data? _Is this possible_? **Edit** :               % Each of the parameters input to the Nielsenneworiginal() are samples from Normal  Distributions.      % A single value from the output of Nielsenneworiginal() is then assigned to thresh_strain and a plot of this will provide a Normal Distribution itself.          function [final_matrix] = MCsolutionsimilar()     tic           no_iterations = input('No. of iterations?:');          thresh_strain = zeros(no_iterations, 1);     %      casechoice =input('Enter 1 for 1st Layup and 2 for 2nd layup:');          diamchoice = input('Enter 1 for Christmas tree-shaped and 2 for Cylindrical:');          [layup, overall_p0] = fibreanglerng(no_iterations, casechoice);          % [E11 E22] = elasticmodulusrng(); % N/m^2     Ef = 235e9; Em = 3.5e9 ; %GPa     [vol_f, overall_p1] = vol_fraction(no_iterations);     E11 = (Ef.*vol_f) + (Em.*(1-vol_f));     E22 = 1./((vol_f./Ef) + ((1-vol_f)./Em));          % Poisson Ratio     nuf = 0.340; num = 0.33;     v12 = (nuf*vol_f) + num*(1-vol_f);          % % Shear Modulus     Gf = 96.52; Gm = 1.8; %GPa     G12 = (Gf*Gm)./(vol_f .* Gm + (Gf.*(1-vol_f)));     [GIC, overall_p2] = shearmodulusstrainenergyrng(no_iterations);          [diamfromimpact_energy, dentsize, overall_p3] = impactenergyrng(no_iterations);          [diam, overall_p4] = diameterng(no_iterations,diamchoice, diamfromimpact_energy);          [time] = detectiontime(no_iterations);          [new_dentsize] = dent_size(dentsize,time);          for i=1:no_iterations         for j=1:16             [J] = Nielsenneworiginal(layup,E11(1,i),E22(1,i),v12(1,i),G12(1,i),GIC(1,i));             if (isreal(J(1,j)))==0                 J(1,j) = sqrt(imag(J(1,j))^2 + (real(J(1,j)))^2);             end                 [thresh_strain(i,1), I] = min(J,[],2);          end     end          %%%%%%% Plot 1 %%%%%%%%%%%%%     figure(1); clf(1)          [mu_j,sigma_j] = normfit(thresh_strain);     x=linspace(mu_j-4*sigma_j,mu_j+4*sigma_j,200);     pdf_x = 1/sqrt(2*pi)/sigma_j*exp(-(x-mu_j).^2/(2*sigma_j^2));     plot(x,pdf_x/10000);     title(sprintf('Mu = %g, sigma = %g',mu_j,sigma_j));     xlabel('Threshold Strains','FontSize',12);     ylabel('Probabilities of occurrence','FontSize',12);     title('\it{Threshold Strains versus Probabilities of occurrence(Christmas tree diameter configuration)}','FontSize',16);          figure(2); clf(2)     X_j = min(thresh_strain) : (max(thresh_strain) - min(thresh_strain))/100 : max(thresh_strain);     P_j = normcdf(abs(X_j),real(mu_j),real(sigma_j));     title('Normal cdf')     plot(X_j,P_j,'blue.-')


title: Boxplots with lognormally distributed data

text: ### Context Environmental data (e.g., pollutant concentrations in water, soil, air) are often lognormally distributed. Even when they are not, we tend to assume that they are (for better or worse). Because of this, 99.5% of the time that I create boxplots, they are presented with a log-scaled concentration axis. Here's an example: ![enter image description here](http://i.stack.imgur.com/jcUyE.png) It occurs to me that computing the lower fence as $Q_1 - (1.5 \times \mathrm{IQR})$ may not be the best way when the data are presumed to be lognormal. ### (Restated) Question Given that, for me, the main value of defining the fences as with $1.5 \times \mathrm{IQR}$ is detecting potential outliers, should I define those fences and outliers in log-space, or keep everything in arithmetic space? I'm currently leaning towards log-transforming the data, but have concerns that this may cause undue confusion or even not be an acceptable practice. ### Similar questions The accepted answer to this question: Is there a boxplot variant for Poisson distributed data? suggests simply transforming the data -- in that case by taking the square root. I'm specifically curious if the fences should be computed in log-space and then converted _back_ to arithmetic space.


title: Whether to use mean or median to summarise the central tendency of length of time to perform a task

text: I work for a commission-based company that gives accounts to employees for about a month and the employees try and resolve the account, if they are successful they get a commission, otherwise the account goes to another employee to try and resolve. We are running some reports to measure some performance characteristics for the company. We have years of data that we can analyze. My question is: **What summary(median or mean) can best describe the central tendency of length of time to resolve an account per employee?** Almost every employee's median is greater than their mean, some are further apart than others. One thing that we are trying to find out is if we give an employee an account, how long should it take for him or her to resolve it. It seems to me that median is a better summary.


title: Good form to remove outliers?

text: I'm working on statistics for software builds. I have data for each build on pass/fail and elapsed time and we generate ~200 of these/week. The success rate is easy to aggregate, I can say that 45% passed any given week. But I'd like to aggregate elapsed time as well, and I want to make sure I don't misrepresent the data too badly. Figured I'd better ask the pros :-) Say I have 10 durations. They represent both pass and fail cases. Some builds fail immediately, which makes duration unusually short. Some hang during testing and eventually time out, causing very long durations. We build different products, so even successful builds vary between 90 seconds and 4 hours. I might get a set like this:               [50, 7812, 3014, 13400, 21011, 155, 60, 8993, 8378, 9100]      My first approach was to get the median time by sorting the set and picking the mid-value, in this case 7812 (I didn't bother with the arithmetic mean for even-numbered sets.) Unfortunately, this seems to generate a lot of variation, since I only pick out one given value. So if I were to trend this value it would bounce around between 5000-10000 seconds depending on which build was at the median. So to smooth this out, I tried another approach -- remove outliers and then calculate a mean over the remaining values. I decided to split it into tertiles and work only on the middle one:               [50, 60, 155, 3014, 7812, 8378, 8993, 9100, 13400, 21011] ->     [50, 60, 155], [3014, 7812, 8378, 8993], [9100, 13400, 21011] ->     [3014, 7812, 8378, 8993]      The reason this seems better to me is two-fold:   * We don't want any action on the faster builds, they're already fine   * The longest builds are likely timeout-induced, and will always be there. We have other mechanisms to detect those So it seems to me that this is the data I'm looking for, but I'm worried that I've achieved smoothness by removing, well, truth. Is this controversial? Is the method sane? Thanks!


title: Sum of Bernoulli variables with different success probabilities

text: Let $x_i$ be independent Bernoulli random variables with success probabilities $p_i$. That is, $x_i=1$ with probability $p_i$ and $x_i=0$ with probability $1-p_i$. Is there a closed expression or an approximate formula for the distribution of the sum $\sum_i x_i$?


title: How to update ARIMA forecast in R?

text: I have a time series data of 30 years and found that ARIMA(0,1,1) has best model among others. I have used the simulate.Arima (forecast package) function to simulate the series into the future.               library(forecast)          series <- ts(seq(25,55), start=c(1976,1))          arima_s <- Arima(series, c(0,1,1))          simulate(arima_s, nsim=50, future=TRUE)      Later on, i have found the updated value of first forecasted year (i.e. series[31] <\- 65). Now i want to simulate the series with this updated value. I am wondering how to do this in R.


title: Is there literature on staggered within subjects designs? What are the consequences of such a design?

text: I am creating a term to describe a design that may not be the proper term ("Staggered within subjects design"). What I mean by this is one can imagine a within subjects design with 3 levels. Each participant provides data for two levels, and between participant the order in which the data is collected is counter-balanced. For example (condition order is in the cells of the table below; sorry for the table formatting, I did my best, but I think you'll get the gist):               Participant Condition A     Condition B     Condition C     1                     1               2      2                     2               1      3                                     1               2     4                                     2               1     5                     1                               2     6                     2                               1     7                     1               2      8                     2               1      9                                     1               2     10                                    2               1     11                    1                               2     12                    2                               1      Does this sort of design have a name? What are the consequences of such a design? How would one analyze it? I understand that I could probably use a mixed model to analyze this data but I'm curious about the application of ANOVA/t-test style analyses to this sort of data as these statistical approaches are more familiar to many in my field. The rationale for such a design might be that you have limited time with each participant and so you can't run them through all three conditions however each participant must be rewarded the same whether they participate in 1 condition or 2 conditions. It seems like, with special consideration, that this data could be run as a between subjects ANOVA (is the same as between subjects t-tests) or a within- subjects ANOVA. I will assert points, and I would appreciate it if the person to answer this question would accept or reject these points.   1. I could conduct a within Subjects anova on this data comparing any two conditions crossing condition with trial order also as a within-subjects factor (I expect that the trial order effects will be non-significant).    2. If there were main effects of order it would influence the validity of my estimate of the exact amount of effect of being demonstrated by each condition but would not influence subsequent analyses.    3. If there were interactions between a condition and order I should not conduct any between subjects analyses on data coming from order 2.   4. I could run two between subjects t-tests for each pair of conditions. One for the data from order 1 and the other from order 2 (as using data from both in a single between subjects analysis would violate the assumption of independence). If I could run all of these statistics, using the same data I would have 3 statistical tests. Each of the between subjects tests would be a full replication of the comparison between condition A and B. The within subjects tests would validate my hypothesis that trial order does not matter and allow the between subjects tests to be applied. If I fail to find significance in my between subjects analysis my within subjects analysis may provide additional power which will allow me to evaluate whether my hypothesis is sufficiently supported to run an additional study or whether I'm barking up the wrong tree.


title: Normalize histogram with different bin width?

text: I know how to normalize histogram (so that the area =1) with the same bin width, but how to do it when the histogram has different bin width? Any idea?


title: Finding a suitable distribution for data in R?

text: I have a set of data with summary as below: ![data summary](http://i.stack.imgur.com/K3fO6.png) The data are non-normal and transformations did not obtain normality. Given the data are continuous I have narrowed the range of distribution down to 'Cauchy','Weibull','Gamma','f','Laplace'. Below is the density plot for the data: ![density plot](http://i.stack.imgur.com/NQ1A4.png) How do I determine the right distribution for the data?


title: Intutive difference between hidden Markov models and conditional random fields

text: I understand that HMM are generative models, and CRF are discriminative models. I also understand how CRFs' are designed and used. What I do not understand is how they are different from HMMs'? I read that in the case of HMM, we can only model our next state on the previous node, current node, and transition probability, but in the case of CRFs' we can do this, and can connect an arbitrary number of nodes together to form dependencies or contexts? Am I correct here?


title: Estimating a distribution based on three percentiles and a mean

text: I have a set of percentiles the 10-th, 50-th and 90-th. Furthermore I have the mean value. I am trying to reconstruct the underlying distribution. The question is similar to Estimating a distribution based on three percentiles except I have a mean value as well. My present approach was to fit a lognormal distribution based on the 50th percentile in which the variance was an unknown. Subsequently I determined the variance_1 based on the 10th percentile keeping the mean as found from the 50th. I did the same for the 90th percentile giving me variance_2. Then I searched for the variance somewhere between variance_1 and variance_2 which gave me the mean I was given in the first place. 10-th percentile: 0.004   50-th percentile: 0.007   90-th percentile: 0.02   arithmetic mean : 0.009 I came up with a log-normal (\mu = log(0.007), \sigma = 0.70825) distribution What do you think of the approach, and how would you tackle this problem?


title: What does the parameter $\alpha$ do in the Jaccard method for binaryRatingsMatrix in R recommenderlab?

text: What is the role of the parameter 'alpha' in the `recommenderlab` `R` package's use of Jaccard method in the recommender model for boolean user- preferences matrix? i.e. `method="Jaccard",param=list(...,alpha=0.5)` ? I saw the code for IBCF ("Item-Based Collaborative Filtering") and they used a dissimilarity function. But the function is not defined in the official CRAN PDF for the `recommenderlab` package. Could someone please help?


title: How to model a multiplicative effect of a parameter

text: I am having difficulty in fitting a model on data. Basically, I have data about the evaluation of phenotypic property (i.e. hard) of 65 palm trees by 5 judges. As an evaluation scheme, each judge provides score to each sample. For 3 judges sample data look like this:               Judge       Product                  Hard          aa             1                      5          ab             1                      6          ac             1                      3          aa             1                      7           ab             1                      5          ac             1                      4          aa             2                      5          ab             2                      8          ac             2                      6          aa             2                      7          ab             2                      4          ac             2                      4       Main objective here is to get product coefficients with less judge errors, for which I want to fit this kind of model: $$Y_{ij} = α_i + β_iθ_j + ε_{ij}$$ _i_ = judge, _j_ = product Here, $α_i$ is judge main coefficients, $_i$ is judge coefficients due to difference in their scoring pattern and $θ_j$ is product coefficients and $ε_i$ is assessor dependent. I was trying to fit this model using `lme` function in R, but difficulty I am facing to fit the interaction term because model here fitted for parameters rather than co-variates. This model looks quite accurate for my kind of data. I have seen Bayesian version (http://www.r-bloggers.com/extending-the-sensory-profiling-data- model/) of it and I don't know how to do using mixed-modelling approach or in a frequentist way. My queries here are: a) What can be an appropriate method to fit this kind of model? I had referred so much literature where description about iterative generalized least squares, multi-level model, separate regression model, weighted least-square model are given. But still I am not getting how to use and fit estimated value of parameters in interaction terms and get separate coefficients for both interaction parameters? b) How can I get heterogeneous error in this form? c) which R package can I use?


title: What factor (trials or blocks) to use in a repeated measures ANOVA

text: I am conducting a psycholinguistic experiment. Each trial consists of the subject responding to a word by pressing a button.The design of my experiment is as follows: 5 blocks of a 100 trials each (each trial is a response). In each block, 50 trials are `Regular` and 50 are `Random`. Response time (RT) is the dependent variable. This experiment is partly a learning experiment. I expect RTs to be faster in the `Regular` condition because there is learning which allows for faster responses. I therefore want to test the difference in mean response time (RT) between the two conditions of my experiment.I know that the model for the ANOVA would be something like               RTs ~ Type of Trial (`Regular` and `Random`) + ...      I have two related questions   1. Should I use `Trial` or `Block` as the other factor? (or both?) Normally I would compare mean RT for each condition across blocks, but it is possible that within blocks, there is already changes in RT between trials. So, is `Block`, `Trial` or both my factor? (besides the condition factor).    2. Since I expect variation across participants, should I include `Subject` as an Error variable in the model? I am doing my analysis in R, in case anyone wants/can provide some advice in that format, but also general statistical advise is much appreciated.


title: What is the definition of "source population"?

text: **What is the definition of source population?** Suppose we screen $100,000$ people and then include $50,000$ people in three studies $A$, $B$ and $C$. For study $A$, is the source population the $50,000$ people or the $100,000$ people screened?


title: Sufficient statistic and hypothesis testing

text: Suppose I have a family of (continuous) distributions $\mathcal{P}=\\{P_\theta(x),\theta\in\mathbb{R}^+\\}$. I also have a statistic $T(x)$ that is sufficient for $\theta$. The value of the parameter $\theta$ is unknown and my interest in it is limited to testing $H_0:\theta=0$ vs. $H_1:\theta>0$. Clearly, the alternate hypothesis $H_1$ is composite, thus Neyman-Pearson lemma doesn't apply. Also, unfortunately, the "easy" uniformly most powerful (UMP) test formulation doesn't apply as I don't think that the family of densities $p_\theta(x)$ has a monotone likelihood ratio (see Lehmann and Romano p. 65). As the matter of fact, while, for any $\theta <\theta'$, the distributions $P_\theta$ and $P_{\theta'}$ are distinct, the ratio $p_{\theta'}(x)/p_\theta(x)$ is a non- _increasing_ (as opposed to non- decreasing) function of $T(x)$. I understand that the maximum likelihood estimator (MLE) $\hat{\theta}_{MLE}$ of $\theta$ is a function of $T(x)$. But what is $T(x)$'s relationship to hypothesis testing? I ask because a trivial test that solves my problem can be constructed as follows: first, pick threshold $\eta$, then compute estimate $\hat{\theta}_{MLE}$, and accept $H_0$ if $\hat{\theta}_{MLE}<\eta$, reject otherwise. The error probabilities then depend on the choice of $\eta$ and the value of $\theta$. However, can anything be said about optimality (in some way, shape or form) of the said test? Are there any conditions on $\mathcal{P}$ other than monotone likelihood function described above that lead to optimality (maybe something weaker than UMP)? The test using $\hat{\theta}_{MLE}$ seems like the best I can do, however, is there a way to prove this statement?


title: What are the major philosophical, methodological, and terminological differences between econometrics and other statistical fields?

text: Econometrics has substantial overlap with traditional statistics, but often uses its own jargon about a variety of topics ("identification," "exogenous," etc.). I once heard an applied statistics professor in another field comment that frequently the terminology is different but the concepts are the same. Yet it also has its own methods and philosophical distinctions (Heckman's famous essay comes to mind). What terminology differences exist between econometrics and mainstream statistics, and where do the fields diverge to become different in more than just terminology?


title: Determining parameters (p, d, q) for ARIMA modeling

text: I am fairly new to statistics and R. I would like to know the process to determine the ARIMA parameters for my dataset. Can you help me figure out the same using R and theoretically (if possible)? The data range from Jan-12 to Mar-14 and portrays the monthly sales. Here's the data set: > 99 58 52 83 94 73 97 83 86 63 77 70 87 84 60 105 87 93 110 71 158 52 33 68 > 82 88 84 And here's the trend: > ![enter image description here](http://i.stack.imgur.com/sFNrZ.png) The data don't exhibit a trend, seasonal behavior or cyclicity.


title: Best metric for evaluation of mixture-of-Gaussian clusters on big-data

text: I have made a new algorithm that is specifically crafted for clustering very large datasets. In order to document it as a research paper, I have to choose one or two **internal (no-label) cluster evaluation measures** to evaluate my algorithm. Which algorithm do you think is generally the best choice for **big datasets**? And **why?** **EDIT:** My algorithm is a modified version of Expectation-Maximization Gaussian mixture models. In other words, the whole data is described by $$P_(x) = \sum_{k=1}^{K}\pi(c_k)\mathcal{N}(\mu_k,\Sigma_k)$$ where $\pi(c_1),\ldots,\pi(c_K)$ are mixture weights. The main difference between my algorithm and regular EM is that it uses some sampling and approximation tricks that accelerate EM. The objective function is the same (the log-likelihood which is to be maximized). Should I use log-likelihood as the evaluation metric? or use other (which?) internal measures for such task? Is it rational?


title: Laplace smoothing parameter choice for Markov chain transitions

text: Let $Y_{t}$ be the state of the process at time $t$, ${\bf P}$ be the transition matrix then: $$ {\bf P}_{ij} = P(Y_{t} = j | Y_{t-1} = i) $$ Since this is a Markov chain, this probability depends only on $Y_{t-1}$, so it can be estimated by the sample proportion. Let $n_{ik}$ be the number of times that the process moved from state $i$ to $k$. Then, $$ \displaystyle\hat{{\bf P}}_{ij} = \dfrac{ n_{ij} }{ \sum_{k=1}^{m} n_{ik} } $$ where $m$ is the number of possible states ($m=5$ in this case). The denominator, $\sum_{k=1}^{m} n_{ik}$, is the total number of movements out of state $i$. Estimating the entries in this way actually corresponds to the maximum likelihood estimator of the transition matrix, viewing the outcomes as multinomial, conditioned on $Y_{t-1}$. However if the observed Markov chain was too small a sparse matrix may arise, suggesting some transitions are impossible: \begin{align} \widehat{P}_{\text{personal care}}=&\begin{pmatrix} 0& 0 & 0& 1.00& 0\\\ 0& 0 & 0.22& 0.44& 0.33\\\ 0& 0.19& 0.57& 0.10& 0.14\\\ 0.03& 0.19& 0.09& 0.50& 0.19\\\ 0& 0.06& 0.13& 0.56& 0.25\\\ \end{pmatrix} \end{align} This is unlikely and so Laplace smoothing could be used to replace mle: $$ \hat{{\bf P}}_{ij} = \frac{ n_{ij}+ \alpha}{ \sum_{k=1}^{m} n_{ik} +n_{ij}\alpha}, $$ where $\alpha$ is some smoothing parameter. What criteria can I use to choose $\alpha$?


title: Linear kernel and non-linear kernel for support vector machine?

text: When using support vector machine, are there any guidelines on choosing linear kernel vs. nonlinear kernel, like RBF? I once heard that non-linear kernel tends not to perform well once the number of features is large. Are there any references on this issue?


title: Dispersion of choices in a survey

text: I have asked a couple of hundred people to pick 5 tasks from a list of some 50 options, with no ordering. > Is there a standard measure of the dispersion of the choices as an unordered > 5-tuple? I am not interested in the frequency with which each option was picked, but on a measure on how much correlation there is between the five things each person picked. I suspect each person makes his choices from a very small subset of the available options, and I'd love to quantify this. In fact, I also asked them to attach to each of the picked options a weight, from 0 to 100, with zero indicating more preference. > Is there a standard measure of dispersion taking into account this extra > information?


title: Best way to handle autocorrelation in residuals from panel data regression (FE)

text: What would be the best way to handle autocorrelated residuals in my panel data analysis (Fixed Effects). What I've tried: -Adding different regressors (no luck there), and besides the current model is theoretically and intuitively justified. -Adding AR(1) term on the right hand side, will take away most of the autocorrelation, but I'm not fond of the new model. I'm working in Eviews but if necessary I could use Matlab or R too. I read a lot about Autocorrelation robust standard errors, applying this to my residuals would solve the problem right? Possibly making coefficients insignificant, but the current result (with autocorrelations) is not that valid. So how can I come up with these autocorrelation robust errors, the easiest way (software-wise)?


title: Distribution/expected length of the shortest path in infinite random geometric graphs

text: Consider an infinite random geometric graph $G(\rho,d)$ in which vertices are uniformly and independently scattered over the 2D plane with density $\rho$ and edges connect the vertices that are closer than $d$. What is the distribution/expected value of the length of the shortest path from a vertex $v_0$ to another vertex $v_1$ which is $k$ hops faraway? _**Note:_** We know that the length of the edges follow the following PDF: $$ f(l)= \begin{cases} \frac{2 l}{d^2} \;\quad l \le d \\\ 0 \qquad\; l > d \end{cases} $$ However, I guess the expected length of the shortest path is not simply $k \times$ the expected length of the edges because in a shortest path longer edges are more likely to be chosen, right?


title: Multilevel model specification for count data (with overdispersion) in R

text: I would like to specify a multilevel model including the following variables:   * DV: count data, i.e. a score value between 0 and 13 (resulting from an additive index) for each individual   * Individual level IVs: 4 socioeconomic variables   * Group level IVs: 2 neighborhood characteristics Data comprise approximately 2000 observations nested in 150 neighborhoods. Of interest are the context effects of the neighborhoods, i.e., their mediating role on the individual level variables. I would thus like to specify a varying slope, varying intercept multilevel model including cross-level interaction terms. My question however relates to the DV. Due to the property of the DV as being count data, I thought of testing for overdispersion and, in case of overdispersion, I would have to take it into consideration when choosing an appropriate distribution. However, I am quite confused regarding the problem of overdispersion in multilevel models and the implementation in `R`, even after having read a lot in this forum and relevant literature. I am currently using the `lme4` package and have two rather basic questions:   1. How do I test for overdispersion?   2. How do I account for it in the model?


title: How to interpret the overlap or non-overlap of two related confidence intervals

text: Say I want to determine if there is a relation/association between birth weights in male and female populations. Hence I calculate the mean for birth weights in males and females and a calculate 95% CI for each using this formula:               CI = mean +/-  1.96 x SE      I then see if the female group birth weight mean lies within the male CI and vice versa. If, in the typical scenario, both genders lie within their opposite genders CI then gender is associated/related with birth weight. However, if not, then gender and birth weight are not associated but are significantly different. My question is what happens if, e.g., the female population lies in the male CI, but the male population does not lie in the female CI, is it possible? Is it due to a statistical error?


title: Visualize before and after score

text: Any suggestions on how to visualize the change in a score before and after an event? I don't want to plot two bars next to each other, though I realize that is the simplest. Is there something more interesting?


title: Expected Value of Quadratic Form

text: On page 9 of Linear Regression Analysis 2nd Edition of Seber and Lee there is a proof for the expected value of a quadratic form that I don't understand. Let $X = (X_i)$ be an $n \times 1 $ random vector and let $A$ be an $ n \times n$ symmetric matrix. If $\mathbb{E}[X] = \mu$ and $Var[X] = \Sigma = (\sigma_{ij})$ then $\mathbb{E}[X^{T}AX] = tr(A\Sigma) + \mu^{T}A\mu$ The problem I have is almost right out the gate, I can't see how $\mathbb{E}[X^{T}AX] = tr(\mathbb{E}[X^{T}AX])$ I think I get the rest of the proof, but if someone here can't point me in the right direction on this part, I'd be forever grateful! Thanks


title: A Conditional Probability Problem

text: This is a biological problem, however the problem requires a statistical solution, so therefore I'll keep it as abstract as possible in order for statisticians to be also able to contemplate it - I'm also only partially sure it requires a conditional probability solution. Basically, I have a bunch of 'miRNAs' and these target certain proteins:               miRNA1 ----> protein1     miRNA1 ----> protein2     miRNA1 ----> protein3     miRNA1 ----> protein4     miRNA1 ----> protein5          miRNA2 ----> protein3     miRNA2 ----> protein4     miRNA2 ----> protein5      There are many different miRNAs that target many different proteins. However, up to two different miRNAs come from the same 'area'. Therefore I want to determine whether given that these miRNAs come from the same area: are their gene targets more similar? Any help on this matter would be extremely useful!


title: Estimating a single proportion from a marginal or conditional model

text: I had a single sample of a binary outcomes (success / failure), and I wanted to estimate the population proportion with a point estimate and a confidence interval. The problem was that some subjects contributed 2 samples, while others just 1, so my observations, or at least some of them, weren't independent. To solve this problem, I ran both a GEE model and a GLMM model. I used the binary outcome as a dependent variable, and a vector of 1's as an independent variable, making the intercept my parameter of interest. I thought that it would give me the proportion I look for. The vector of 1's was declared in SAS as a classification variable. And now for the problem. When I used the `logit` link function, I did not get the correct results. When I used the log link function, I did. I thought that the `log` link function suits the Poisson model, so I am confused. I suspect it might have something to do with the difference between proportion and predicted probabilities, but I'm not sure. I wanted to ask if someone can explain to me (and maybe please attach a formula) why the log link function is what I need here? To be more specific, if $b$ is the estimate of the intercept, I did at first $e^b/(1+e^b)$, which was wrong. When I did $e^b$, it was correct. I need to know why... My SAS code was:               PROC GENMOD DATA=data descending;         CLASS ID D * D is vector of 1's;         MODEL Y  = D / dist = bin  link=log cl;         REPEATED subject= ID;         LSMEANS D/ cl exp;     RUN;


title: Frailty conundrum: how can area effects be large when shared frailty variance for area is non-sig?

text: I’m running Cox proportional hazards models (in stata) to predict risk of teenage birth, where individuals (c. 50 000, c. 3000 births) are clustered within wards (>500). In null models with shared frailty at ward level there is large variance in the outcome between wards. When I control for individual characterisitcs I find that the variance term for frailty becomes very low, close to 0, and non- significant. AIC implies models without shared frailty are a better fit. However, when I then add ward-level variables to a model without shared frailty, ward level variables have large effects and have lower AIC values than the ind variable only model, despite the fact that the previous comparison between non-shared and shared frailty models showed that there was no significant variation in the outcome between wards. I am confused as to why this is. Can anyone explain why this happens, or direct me to some relevant literature? I don’t know if it matters that my outcome is quite rare, and there isn’t a success/birth event in every ward. My ward level variables are continuous (mortality rate, crime rate etc) but made into quartiles, and are significant only at the highest levels. Many thanks, Caroline


title: What's the minimum number of measured variables required when fitting a one latent variable model?

text: I am envisaging a situation in which I have designed a small questionnaire, and want to assess whether the measured variables are measuring the same thing. I'm not really interested in practical tips so much as an statement of what the minimum number is to obtain an identified model and an explanation of why that is the minimum.


title: Compare means - heterogeneous variance, non-normal

text: I would like to perform some basic means comparisons (3 conditions, balanced, n = 21 for a total of 63 observations). The main dependent variable is total task completion time, i.e., the experimental question of interest is: Does task completion time differ across the 3 experimental conditions? The data are non-normal and the variances are not homogeneous across the groups (smallest variance is about 10 times the largest). Due to the large difference in variance, I understand that a non-parametric alternative (like Kruskal-Wallis) would likely be inappropriate. A log transformation of task completion time does improve the situation--normalizes the residuals in 2 of the 3 conditions and results in homogeneous variances. My main question is: Is a log transformation justified here? If so, could I still use ANOVA after the transformation, even though 1 of the conditions is still non-normal?


title: hazard rate and cumulative incidence function

text: Assume an event of interest $E$ and the competing risk of death $D$. Let $h(t)$ the hazard rate $h(t)=\frac{Pr\left\\{(t < T \leq t + \Delta t)| (T > t)\right\\}}{dt}\\\$ and $S(t)$ the survival function If we calculate the cumulative incidence function $CIF$ for $E$, as $CIF=Pr(T<t)=\int_{0}^{t}hSdu$ What is the relationship between the $h(t)$ and the $CIF´=\frac{CIF}{dt}\\\$ ? Is then $CIF´={Pr(t<T<T+dt)}\\\$?


title: Dealing with outliers when comparing variances with Bartlett's test

text: I have four different groups (with unequal sample sizes of 100 to 120) and want to test if the variance differs. For ANOVAs I used the winsorized mean to get a more robust estimate and I am wondering if it is recommended to use the winsorized data (all values extremer than the 5th or 95th percentile are set to the nearest percentile) when testing whether the variances are equal. (e.g. with a Levene test). Or would you work with a different kind of outlier removal? I checked the normality assumptions of the groups and saw that it stays approximately the same (some slight violations of normality) when I compare the tests of the original data with the tests of the winsorized one. Any thoughts on that would be highly appreciated.


title: calculation of 12 week PFS rate

text: what type of setup is needed to calculate 12 weeks PFS rate using Kaplan Meier method in SAS? Do we calculate PFS at 84 days and also the patients status? alive, dead, lost to follow-up. then what? Thanks!


title: Examining the relative efficiency of the trimmed mean

text: Using the following reference A survey of sampling from contaminated distributions, I am trying to investigate the relative efficiency (RE) for the mean vs the trimmed mean, given the following contaminated exponential distribution \begin{align} F(x) = F_\lambda(x) + F_{k\lambda} (x) \end{align} When trying to simulate the RE for the this distribution in R for trim = 0, 0.01, 0.03 and 0.05, I am getting the following plots ![Plots for the Asymptotic Relative Efficiency of the trimmed mean](http://i.stack.imgur.com/rSw4Y.png) Unfortunately, like the reference above, I suspect that the plots should crossover each other at some point, but this is not happening. So far I have come up with the following code               number.of.simulations <- 100     n <- 1000     error <- 0          #using random uniform variables to create exponential distribution     uniform.variables <- runif(n)          no.trim.mean <- trim.mean <- rep(0,number.of.simulations)          uncensored <- rep(0,n)     fraction.of.trim <- seq(from = 0 , to = 0.1, by = 0.001)          amount.trim <- c(0, 0.01, 0.03, 0.05)     relative.efficiency <- matrix(         NA,nrow= length(fraction.of.trim) , ncol = length(amount.trim)     )          for(i in 1:length(amount.trim)){       for(j in 1:length(fraction.of.trim)){          uncensored <- matrix(             (uniform.variables <= 1 - error) * rexp(n = n * number.of.simulations) +              (uniform.variables > 1  - error) * rexp(n = n * number.of.simulations),             ncol = number.of.simulations          )             no.trim.mean <- apply(uncensored, 2, mean, trim = amount.trim[i])         trim.mean <- apply(uncensored, 2, mean, trim = fraction.of.trim[j])            relative.efficiency[j,i] <- mean((no.trim.mean - 1)^2)/mean((trim.mean - 1)^2)       }       }      What is it that I am missing or doing wrong with this simulation?


title: Why would you use logistic regression for count data?

text: Wouldn't it be better to use poisson regression for count data? Also in logistic regression what is the advantage of using the log link versus the logit link? I know you can get the log relative risk with the log link. But why use relative risks as opposed to odds ratios?


title: Expected number I will be on after drawing cards until I get an ace, 2, 3, and so forth

text: I am having some trouble solving the following. You draw cards from a standard 52-card deck without replacement until you get an ace. You draw from what is remaining until you get a 2. You continue on with 3. What is the expected number you will be on after the entire deck runs out? It was natural to let   * $T_i = \text{first position of card whose value is }i$   * $U_i = \text{last position of card whose value is }i$ So the problem essentially amounts to figuring out the probability that you will be on $k$ when the deck runs out, namely: $$Pr(T_1<\cdots<T_k \cap U_{k+1} < T_k)$$ I can see that $$Pr(T_1<\cdots<T_k) = 1/k! \\\ \text{and} \\\ Pr(U_{k+1} < T_k) = 1/70$$ but could not get any further...


title: Transforming nonlinearity in genearlized linear models

text: Are the buldging rules applicable in generalized linear models? Specifically, to transform the independent variables? I've only seen it disscused/used in OLS regression. Thanks in advance


title: Mixed design anova assumptions

text: What are the assumptions that I need to consider when running a mixed design ANVOA? Let's say I have the following data set:                                |       Scenario 1       |     Scenario 2         |                      |Trial 1|Trial 2| Trial 3|Trial 1|Trial 2| Trial 3|      -------------------------------------------------------------------                   S1 | ...      Condition 1  S2 | ...                   S3 |      -------------------------------------------------------------------                   S5 |      Condition 2  S6 |                   S7 |      Where `Condition x` is the between-subject variable, `S1, S2, ...` are the subjects, and `Scenario`s are within-subject variables, where each `Scenario` has some nested `Trial`s. Using `lmer` I can describe this model by: `lmer(value~Condition*Scenario + (1+Scenario|Player) + (1|Scenario/Trial)`. Taking a more traditional approach I can collapse over `Trials` so I have the a mean value for each `Scenario`. I think I should do the following:   1. I should check Homoskedasticity per `Scenario/Trial`   2. Check for the absence of influential data points Is there something I am missing?


title: Two studies, opposite results: Am I allowed to compute an overall mean using random-effect model?

text: I am performing a meta-analysis. I have a subgroup made only of 2 studies and they show opposite results. The 1st has excellent outcomes the 2nd very bad ones. Is it true that I am not allowed to compute the data from these two (different) studies under a random-effect model? I have been told that it is methodologically not correct to perform a meta-analysis of 2 studies with such opposite results. They suggested me to report the raw data from these studies separately (without an overall outcome calculated by a random-effect model). Any suggestion? Or any place where I can find the explanation about this point?


title: Clustering a noisy data or with outliers

text: I have a noisy data of two variables like this.               x1 <- rep(seq(0,1, 0.1), each = 3000)     set.seed(123)     y1 <- rep (c(0.2, 0.8, 0.3, 0.9, 0.65, 0.35,0.7,0.1,0.25, 0.3, 0.95), each = 3000)     set.seed(1234)     e1 = rnorm(length(x1), 0.07,0.07)     set.seed(1223)     e2 = rnorm(length(x1), 0.07,0.07)     set.seed(1334)     yn <- rnorm(20000, 0.5,0.9)     set.seed(2344)     xn <- rnorm(20000, 0.5,0.9)     y <- c(y1 + e1,yn)      x <- c(x1 + e2, xn)      plot(x,y,  xlim=c(0,1.2), ylim = c(0,1.2), pch = ".", col = "gray40")       I can visually see there are potential 10 clusters in closer look. ![enter image description here](http://i.stack.imgur.com/fkHBd.jpg) However the whole data has much points spread:               plot(x,y,   pch = ".", col = "gray40")       ![enter image description here](http://i.stack.imgur.com/uKAYl.jpg) I would like to make 10 clusters. I tried K-means cluster analysis.               xm1 <- cbind(x,y)     cl1 <- kmeans(xm1, 10)     colrs <- c("red", "green", "blue1", "pink", "green4","tan",       "gray40", "yellow", "black", "purple")      plot(xm1, col = colrs[cl1$cluster], pch = ".", xlim=c(0,1.2), ylim = c(0,1.2))      ![enter image description here](http://i.stack.imgur.com/gJj8C.jpg)               plot(xm1, col = colrs[cl1$cluster], pch = ".")      ![enter image description here](http://i.stack.imgur.com/F0x1w.jpg) Is there anyway (may be kernel k-means, nearest neighbors) that can do more justice to this type of data. If so how can I do this ?


title: Single time-series difference before and after treatment

text: I am doing panel data analysis. I have N= 103 firms in T= 5 years (that is around 507 observations with an “unbalance data” fixed effect). The study period is from 2008-2012. The period of the financial crisis is 2008. I am not sure whether I can use a differences in differences model to get the effect of the crisis. I can’t see a control group because all companies were affected by the crisis. Michael R. Roberts talks in his book "Endogeneity in Empirical Corporate Finance" about single time-series difference before and after treatment. Does he mean to interact my independent and moderating variables with a crisis dummy which is the year 2008? And is it right to interact the crisis dummy with all the observations to get the crisis effect?


title: Wrong values in Dunnett post-hoc test

text: I can't get the same results in R as in GraphPad Prism for repeated measures anova. The experiment was a stimulation time course, so I have as DV=response and as factor "time" within groups, also I add a factor sample for each experiment               data <- read.csv("http://dl.dropbox.com/u/4828275/datos.csv")     options(contrasts=c("contr.sum","contr.poly"))          ## Convert variables to factor     data <- within(data, {     sample <- factor(sample)     time <- factor(time)     })     aov <- aov(response~time+sample, data=data)     summary(glht(aov, linfct=mcp(time="Dunnett")))               Simultaneous Tests for General Linear Hypotheses          Multiple Comparisons of Means: Dunnett Contrasts               Fit: lme.formula(fixed = response ~ time, data = data, random = ~1 |      sample)          Linear Hypotheses:                 Estimate Std. Error z value Pr(>|z|)     2 - 0 == 0    1.1789     2.0800   0.567        1     5 - 0 == 0    1.2966     2.0800   0.623        1     10 - 0 == 0   1.0555     2.0800   0.507        1     15 - 0 == 0   0.4317     2.0800   0.208        1     30 - 0 == 0   0.2148     2.0800   0.103        1     (Adjusted p values reported -- bonferroni method)      For repeated measures I have this code               aov.repeated <- ezANOVA(       data       , dv = .(response)       , wid = .(time)       , within = .(sample)       , type = 1       , return_aov = TRUE     )$aov      The GraphPad Prism results for the same data was               Table Analyzed  Data 1                        Repeated Measures ANOVA                        P value   0.0415                     P value summary   *                      Are means signif. different? (P < 0.05)   Yes                    Number of groups  6                      F 2.863                      R square  0.4172                         Was the pairing significantly effective?                          R square  0.1980                     F 2.119                      P value   0.1162                     P value summary   ns                     Is there significant matching? (P < 0.05) No                        ANOVA Table SS  df  MS             Treatment (between columns)   130.6   5   26.12              Individual (between rows) 77.30   4   19.32              Residual (random) 182.4   20  9.121              Total 390.3   29                    Dunnett's Multiple Comparison Test  Mean Diff.  q   Significant? P < 0.05?  Summary 95% CI of diff       0 vs 2    -2.861  1.498   No  ns  -8.085 to 2.362       0 vs 5    -5.777  3.024   Yes *   -11.00 to -0.5531       0 vs 10   -6.009  3.146   Yes *   -11.23 to -0.7855       0 vs 15   -4.621  2.419   No  ns  -9.844 to 0.6029       0 vs 30   -2.581  1.351   No  ns  -7.805 to 2.642      How can I get the same results as above in R? Is there a way to get Dunnett's Multiple Comparison Test in aov.repeated?


title: Approximation in a equality

text: Let $X_1,X_2$ be independent distributed with cdfs $F_1(x),F_2(x)$, so that $$\overline F_i(x) := 1 - F_i(x) = x^{-\alpha}L_i(x),\ \alpha \geq 0 \>,$$ where $L_i(x)$ is a slowly varying function, that is, it satisfies $L(tx)/L(x) \to 1$ as $x \to \infty$ for all $t > 0$. For $0<\delta<\frac{1}{2}$ I have the following line: $\overline{F_{1}}((1-\delta)x)+\overline{F_{2}}((1-\delta)x)+\overline{F_{1}}(\delta x)\overline{F_{2}}(\delta x) =\left(\overline{F_{1}}((1-\delta)x)+\overline{F_{2}}((1-\delta)x)\right)\left(1+o(1)\right).$ Can someone explain why we have $o(1)$ as $x \rightarrow \infty$ in the second line of the equation?


title: Does the samples mean and the histogram mean have to be the same?

text: I'm reading this book. Here, he gives me a method for calculating the mean: ![enter image description here](http://i.stack.imgur.com/ipzQY.gif) Here he introduces a new method of calculating the mean through histogram: ![enter image description here](http://i.stack.imgur.com/vSuRl.gif) He says that now we will count and index the samples in the histogram, this list: > {7, 8, 0, 9, 2, 10, 4, 1, 6, 0, 8, 7, 6, 0, 0} Is going to be counted and it will result in: > {{7, 2}, {8, 2}, {0, 4}, {9, 1}, {2, 1}, {10, 1}, {4, 1}, {1, 1}, {6, 2}} It means there are two 7's in the list, there are two 8's in the list and so on. Then i've calculated the mean of this list: > {4, 1, 1, 0, 1, 0, 2, 2, 2, 1, 1} Which is almost as the ordered version of the last list, four 0's, one 1 and so on The problem is that the mean of the first one is different of the second, while the mean of the first list is 4.53333, the second list is 1.36364, do these means have to be identical?


title: Maximum likelihood estimate: Is this possible to solve?

text: I have the following problem: Formulate the likelihood function, the log-likelihood function, and the maximum-likelihood estimate as well as the Fisher information and the observed Fisher information for each of the following problems. c) $X_i,\dots,X_n \overset{\mathrm{i.i.d}}{\sim} \mathcal{N}(\theta,\theta^2)$ d) $X_i,\dots,X_n \overset{\mathrm{i.i.d}}{\sim} \mathcal{N}(\theta,\theta)$ I have the following so far: ![enter image description here](http://i.stack.imgur.com/O66Jj.png) First, did I do the right things so far? And more important: How can I calculate the maximum likelihood estimate? Or is this even possible? * * * **EDIT** After CoolSerdash's hint I got the following: ![enter image description here](http://i.stack.imgur.com/RlKZQ.png) Is that correct? But I have some doubts to formulate the observed Fisher information. Wouldn't this term get a beast inserting the maximum-likelihood estimate? * * * **Similar question** I have to do the same for a Gamma distribution: $X_i,\dots,X_n \overset{\mathrm{i.i.d}}{\sim} \Gamma(\alpha,\beta)$ where $\alpha$ is known. I have the following so far: ![enter image description here](http://i.stack.imgur.com/5Jatl.png) Does this make sense?


title: Why does my proof for showing that the Kaplan-Meier estimate is unbiased not work?

text: I know that the Kaplan-Meier estimator is biased because my textbook says so. However, I don't understand why the following proof doesn't work: Let $\hat{S}(t)$ be the Kaplan-Meier estimate for the survival function $S(t)\equiv P(T_i > t)$ where $T_i$ are iid failure times. Let $\hat{\Lambda}(u)$ be the Nelson-Aalen estimator for the cumulative hazard function $\Lambda(u)$. It is known that $\frac{\hat{S}(t)}{S(t)}-1 = -\int\limits_{0}^{t}\frac{\hat{S}(u^-)}{S(u)}d\\{\hat{\Lambda}(u)-\Lambda(u)\\}$. Now, $\int\limits_{0}^{t}\frac{\hat{S}(u^-)}{S(u)}d\\{\hat{\Lambda}(u)-\Lambda(u)\\}$ is a martingale because $\hat{\Lambda}(u)-\Lambda(u)$ is a martingale and because $\frac{\hat{S}(u^-)}{S(u)}$ is a predictable process. So, $\mathbb{E}[\frac{\hat{S}(t)}{S(t)}-1]=0 \Rightarrow \mathbb{E}[\frac{\hat{S}(t)}{S(t)}]=1 \Rightarrow \mathbb{E}[\hat{S}(t)]={S(t)}$ since $S(t)$ is a non-stochastic function.


title: PDF restoration from samples — what is the correct name of the method?

text: I have an independent sample $x_1 \ldots x_N$, identically distributed. I make an empirical CDF as $P_{\mathrm{emp}}(x)=\sum\limits_{i=1}^N H(x-x_i)$, where $H(x)$ is a Heaviside step function. Then I make an interpolation of the raw CDF as $$P_{\mathrm{interpolated}}(x) = \frac{P_{\mathrm{emp}}(x_{i+1}) - P_{\mathrm{emp}}(x_i)}{x_{i+1}-x_i}\cdot(x_{i+1}-x),$$ $$x \in [x_i, x_{i+1}]$$ Then I take a known window function $w(x)$, convolve $P_{\mathrm{interpolated}}(x)$ with it: $$P_{\mathrm{smoothed}}=[P_{\mathrm{interpolated}} \ast w](x)$$ Finally, I take the PDF as the derivative of $P_{\mathrm{smoothed}}(x)$ $$f_{\mathrm{PDF}}(x)=P_{\mathrm{smoothed}}'(x)$$ Can you tell me the proper name of such a method of PDF restoration? Or, maybe, the name of the wide class of PDF restoration methods?


title: Research on debt recovery

text: My final year project is on debt recovery data for a debt collection firm. Data such as original/current balances,payments made,DOB, number of contacts made,whether or not a debtor has made insuarance payments are given We need to model the given problem in order to predict the collectability of debts of a given debtor. Only approach I started on was a logit model. Any other ideas or hint on possible ways to approach the problem would be appreciated.


title: R-squared vs. t-test confusion

text: I am using regression tree for this educational exercise. I have 500 observation points of independent variables $X$ and dependent variable $Y$. I use 400 points to generate my model and keep 100 points for out of sample (OOS) testing. I then, use the model to predict new Y (Y_predicted) for OOS data. I then, obtain R-squared using original Y from OOS data and Y_predicted. I also perform a t-test on (Y - Y_predicted). The R-squared comes out to be negative and t-test fails to reject null at 5%. The null being the mean difference of Y and Y_predicted is zero. I am confused about how to interpret this result. R-squared is saying my model is basically bogus but t-test is saying that since mean difference of Y and Y_predicted is 0 at 5% level, so the model is appropriate. R-squared is calculated on OOS original Y and predicted Y using model. t-test was run on OOS (Y - Y_predicted)


title: Distribution for percentage data

text: I have a question about the correct distribution to use for creating a model with my data. I conducted a forest inventory with 50 plots, each plot measures 20m × 50m. For each plot, I estimated the percentage of tree canopy that shades the ground. Each plot has one value, in percent, for canopy cover. Percentages range from 0 to 0.95. I am making a model of percent tree canopy cover ( _Y_ variable), with a matrix of independent _X_ variables based on satellite imagery and environmental data. I am not sure if I should use a binomial distribution, since a binomial random variable is the sum of _n_ independent trials (i.e., Bernoulli random variables). The percentage values are not the sum of trials; they are the actual percentages. Should I use gamma, even though it does not have an upper limit? Should I convert percentages to integer and use Poisson as counts? Should I just stick with Gaussian? I have not found many examples in the literature or in textbooks that try to model percentages in this way. Any hints or insights are appreciated. * * * Thank you for your answers. In fact, the beta distribution is exactly what I need and is thoroughly discussed in this article:   * Eskelson, B. N., Madsen, L., Hagar, J. C., & Temesgen, H. (2011). Estimating Riparian understory vegetation cover with Beta regression and copula models. Forest Science, 57(3), 212-221.  These authors use the betareg package in R by Cribari-Neto and Zeileis. The following article discusses a good way to transform a beta-distributed response variable when it includes true 0's and/or 1's in the range of percentages:   * Smithson , M., and J. Verkuilen, 2006. A better lemon squeezer? Maximum-likelihood regression with beta-distributed dependent variables, Psychological Methods, 11(1):54–71.


title: Conditional vs unconditional expectation

text: I'm having trouble understanding the calculation of conditional versus unconditional expectations in this case: \begin{array}{c c c c} & \quad &X=1\quad &X=-1\quad &X=2 \\\ &Y=1\quad &0.25\quad &0.25\quad &0 \\\ &Y=2\quad &0\quad &0\quad &0.5 \end{array} \begin{align} ~ \\\ ~ \\\ ~ \\\ E(Y)&=\sum_Y y f(y) \\\ ~ \\\ ~ \\\ ~ \\\ E(Y|X)&=\sum_Y y f(y|x) \end{align} To me, both calculations are $1*0.25 + 1*0.25 + 2*0.5 = 1.5$. What am I doing wrong?